{"paper_id":"1f13149f-0f22-4f42-8596-7fc01b297bfc","body_text":"Towards high-resolution modeling of small molecule – ion channel \ninteractions \nBrandon J. Harris1,2, Phuong T. Nguyen1, Guangfeng Zhou3, Heike Wulff4, Frank DiMaio3, 1 \nVladimir Yarov-Yarovoy1,2,5* 2 \n1Department of Physiology and Membrane Biology, University of California Davis, Davis, 3 \nCalifornia, United States 4 \n2Biophysics Graduate Group, University of California Davis, Davis, California, United States 5 \n3Department of Biochemistry, University of Washington, Seattle, Washington, USA; Institute for 6 \nProtein Design, University of Washington, Seattle, Washington, USA. 7 \n4Department of Pharmacology, School of Medicine, University of California Davis, Davis, 8 \nCalifornia, United States 9 \n5Department of Anesthesiology and Pain Medicine, University of California Davis, Davis, California, 10 \nUnited States 11 \n* Correspondence:  12 \nVladimir Yarov-Yarovoy 13 \nyarovoy@ucdavis.edu 14 \nKeywords: Rosetta, ligand docking, ion channel, computational modeling, ligand – protein 15 \ninteractions, computer-aided drug discovery.  16 \n 17 \nAbstract 18 \nIon channels are critical drug discovery targets for a wider range of pathologies, such as epilepsy, 19 \nchronic pain, itch, autoimmunity, and cardiac arrhythmias. To develop effective and safe 20 \ntherapeutics, it is necessary to design small molecules with high potency and selectivity for specific 21 \nion channel subtypes. There has been increasing implementation of structure-guided drug design for 22 \nthe development of small molecules targeting ion channels. We evaluated the performance of two 23 \nRosetta ligand docking methods, RosettaLigand and GALigandDock, on structures of known ligand - 24 \ncation channel complexes. Ligands were docked to voltage-gated sodium (NaV), voltage-gated 25 \ncalcium (CaV), and transient receptor potential vanilloid (TRPV) channel families. For each test case, 26 \nRosettaLigand and GALigandDock methods were able to frequently sample a ligand binding pose 27 \nwithin 1-2 Å root mean square deviation (RMSD) relative to the native ligand – channel structure. 28 \nHowever, RosettaLigand and GALigandDock scoring functions cannot consistently identify native 29 \ndrug binding poses as top-scoring models. Our study reveals that the proper scoring criteria for 30 \nRosettaLigand and GALigandDock modeling of ligand - ion channel complexes should be assessed 31 \non a case-by-case basis using sufficient ligand and receptor interface sampling, knowledge about 32 \nstate specific interactions of the ion channel and inherent receptor site flexibility that could influence 33 \nligand binding. 34 \n 35 \nIntroduction 36 \nThe voltage-gated cation channel families consist of pore-forming transmembrane proteins that 37 \nselectively conduct ions across lipid bilayers, and mediate physiological processes such as signal 38 \ntransduction, gene expression, synaptic transmission, and the activation and proliferation of cells in 39 \nthe immune system (Catterall, 1995; Hille, 2001; Catterall, 2011; Feske, Wulff, and Skolnik, 2015; 40 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\nNanou and Catterall, 2018). Cation channels function in a finely regulated manner across spatial and 41 \ntemporal domains to complete these cellular functions. Current drug discovery efforts aim to 42 \nmodulate channel activity by targeting specific channel domains. For instance, a voltage-gated 43 \nsodium channel's therapeutically relevant structural domains are the selectivity filter (otherwise 44 \nknown as the outer pore vestibule), the central pore cavity (otherwise known as the inner pore 45 \nvestibule), and the voltage-sensing domain (Nguyen and Yarov-Yarovoy, 2022). There have been 46 \nconsiderable academic and industry efforts to identify therapeutically relevant small molecules that 47 \nselectively target ion channels (Bagal et al., 2013; Zhang et al., 2020). However, developing effective 48 \nand safe therapeutics targeting ion channels has been challenging (Wulff et al., 2019).  49 \nTo address these challenges, drug discovery pipelines are trending towards incorporation of 50 \nvirtual drug screening and computer-aided drug design processes, for their ability to minimize drug 51 \ndevelopment time and cost (Maia et al, 2020). Among these processes, molecular docking has 52 \ndemonstrated its usefulness for structure-based drug discovery. Molecular docking involves 53 \npredicting conformations of the small molecule’s orientation with the protein (known as the pose) 54 \nand scoring the poses to rank the likely protein-ligand interaction (Meng et al, 2011). 55 \nAmong the numerous molecular docking software packages, Rosetta is a protein modeling 56 \nsoftware and design suite with two established small molecule docking methods: RosettaLigand 57 \n(Meiler and Baker, 2006; Davis and Baker, 2009; Smith and Meiler, 2020) and GALigandDock (Park 58 \net al, 2021). RosettaLigand uses a Monte Carlo minimization procedure using the Rosetta energy 59 \nfunction (Alford et al, 2017) to dock a pre-generated set of ligand conformers, while allowing side 60 \nchain flexibility within a protein-receptor site. GALigandDock utilizes a different approach with two 61 \ndistinct features. First, the scoring function, RosetteGenFF, is a new generalized energy function 62 \ntailored for small molecules. RosettaGenFF was trained from the Cambridge Structural Database 63 \n(Groom et al, 2016), which at the time contained 1,386 small molecule crystal lattice arrangements to 64 \ncreate a balanced force field that discriminates true lattice packing arrangements of the ligand from 65 \ndecoy (alternative lattice packing and conformational) arrangements. During docking, an orientation-66 \ndependent water-bridging energy term is incorporated within RosettaGenFF to further discriminate 67 \nthe protein-ligand orientation (Pavlovicz, Park, and DiMaio, 2020). Second, GALigandDock samples 68 \nconformational space using a genetic algorithm. The ligand rigid body degrees of freedom and 69 \nrotatable torsions are encoded as ‘genes’ to generate new ligand inputs for successive docking 70 \niterations. This allows efficient sampling of the protein-ligand energetic landscape when paired with 71 \nthe RosettaGenFF score function, canonical Monte Carlo optimization, and quasi-Newtonian 72 \nminimization procedures within the Rosetta framework (Park et al, 2021).  73 \nWhile Rosetta protein-ligand docking methods are able to perform well with soluble protein – 74 \nligand benchmarks, the application of these methods to membrane-embedded ion channels has not 75 \nbeen explored. Since there is a need to better assess and screen small molecules targeting different 76 \nion channel domains, we selected a diverse set of ten known cation channel-ligand structures for 77 \nevaluation. From this set of high-resolution ion channel – small-molecule complexes, we assessed the 78 \naccuracy of the RosettaLigand and GALigandDock methods in sampling ligand poses near the native 79 \nstructural coordinates and predicting the closest matching pose by energy ranking. 80 \nOur case studies include four voltage-gated sodium (NaV) channel structures, five voltage-gated 81 \ncalcium (CaV) channel structures, and one transient receptor potential vanilloid (TRPV) channel 82 \nstructure. The ion channel ligand binding sites include the voltage-sensing domain, the selectivity 83 \nfilter, and the central pore cavity. Our results demonstrate that RosettaLigand and GALigandDock 84 \nmethods can frequently sample ligand binding poses within 1-2 Å root-mean-square deviation 85 \n(RMSD) from the native channel – ligand structure. However, the ability to identify a near-native 86 \npose from energy ranking remains a challenge. When considering factors like the targeted ion 87 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n \n \nchannel domain, the ligand library features, and the sampling of ligand and receptor site 88 \nconformations, our work demonstrates that high-resolution structures paired with RosettaLigand or 89 \nGALigandDock can support drug discovery pipelines on a case-by-case basis. 90 \n 91 \nMaterial and Methods 92 \nLigand generation 93 \nLigands were extracted as Structure Data Files (.sdf) from PubChem (Kim et al, 2023). Using the 94 \nAvogadro software (Hanwell et al, 2012), each ligand structure underwent bond-correction, 95 \nprotonation at pH 7.4, and energy minimization using the Merck Molecular Force Field (Halgren, 96 \n1996a; b; c; Halgren and Nachbar, 1999; Halgren, 1999a; b). The resulting models were saved as 97 \nTripos Mol2 (.mol2) files. The protonation and bond order of saxitoxin and tetrodotoxin was matched 98 \nto experimentally reported work (Hinman and DuBois, 2003; Thomas-Tran and Du Bois, 2016). Both 99 \nexperimentally resolved structures of verapamil docked to rabbit CaV1.1 were tested (Zhao et al, 100 \n2019). 101 \nNext, using Amber Tools’ Antechamber protocol, the partial atomic charge, atom, and bond type 102 \nassignments for each ligand were AM1-BCC corrected (Case et al, 2021; Salomon-Ferrer, Case, and 103 \nWalker, 2012; Appendix S1). AM1-BCC correction is commonly used in Rosetta-based ligand 104 \ndocking protocols (Smith and Meiler, 2020; Park et al, 2021) and has demonstrated a similar 105 \nperformance correlation with other RosettaLigand input preparation protocols (Smith and Meiler, 106 \n2020). 107 \nThe AM1-BCC corrected ligands were used for subsequent steps specific to each method. For 108 \nRosettaLigand, an in-house script (Appendix S2) using the OpenEye Omega toolkit (Hawkins et al, 109 \n2010) was used to generate the conformer library, followed by Rosetta to generate the associated 110 \nligand parameters file. For GALigandDock, the input conformer was generated using the 111 \nRosettaGenFF crystal structure prediction protocol (Park et al, 2021; Appendix S3), taking the 112 \nlowest energy packing arrangement as input. 113 \nIon channel preparation 114 \nIon channel structures were downloaded from the Protein Data Bank (Berman et al, 2000). Prior 115 \nto RosettaLigand docking, structures were relaxed with backbone constraints using the RosettaRelax 116 \nprotocol (Nivón, Moretti, and Baker, 2013). This protocol allows the repacking of protein sidechains 117 \nand minimization of the structure into the Rosetta score function for comparison between poses. The 118 \nlowest energy pose from 100 relaxed poses was used for docking. 119 \nRosettaLigand docking 120 \nRosettaLigand docking was performed using previously described RosettaScripts protocols 121 \n(Davis and Baker, 2009; DeLuca, Khar, and Meiler, 2015; Appendix S4-S7). Briefly, the initial 122 \nplacement of all ligands into their respective ion channels was performed by superimposing the initial 123 \nligand poses onto the experimental ligand coordinates. RosettaLigand uses grid-based sampling to 124 \nscore the ion channel – ligand interface, whereby all ligand atoms and ion channel atoms a set 125 \ndistance from the ligand are scored with the defined scoring function. The scoring grid width was 126 \ncalculated uniquely for each ligand to ensure all ligand atoms were bound to the scoring grid. As 127 \nused previously, the scoring grid width was calculated as the maximum conformer atom-atom 128 \ndistance plus twice the box size value used in the Transform mover (Moretti, Bender, Allison, and 129 \nMeiler, 2016; Appendix S4). 130 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\nRosettaLigand has a low-resolution and high-resolution sampling phase. For an unbiased 131 \nsampling of the domain, an initial transformation during the low-resolution sampling phase is 132 \nperformed on the initial ligand position using Rosetta’s Transform mover with a box size of 7-8 Å. 133 \nThe low-resolution Transform step is a grid-based Monte Carlo simulation, where the ligand is 134 \ntranslated up to 0.2 Å and rotated up to 20 degrees per iteration, for a total of 500 iterations. The pose 135 \nwith the lowest score is then used as the starting pose for high-resolution docking. For high-136 \nresolution docking using Rosetta’s HighResDocker mover, six docking cycles of rotamer 137 \nsampling are performed with repacking of sidechains every third iteration. Lastly, the ion channel – 138 \nligand complex is minimized using Rosetta’s FinalMinimizer mover, and the interface scores 139 \nare reported using the InterfaceScoreCalculator.  140 \nFor all RosettaLigand docking runs, the ligand_soft_rep and hard_rep scoring functions 141 \nwere reweighted based on previous work assessing Rosetta score functions with the Comparative 142 \nAssessment of Score Function 2016 (CASF-2016) dataset (Su et al., 2019; Smith and Meiler, 2020). 143 \nSpecifically, reweights of Coulombic electrostatic potential (fa_elec), Lennard-Jones repulsive 144 \nenergy between atoms in the same residue (fa_intra_rep), sidechain-backbone hydrogen bond 145 \nenergy (hbond_bb_sc), sidechain-sidechain hydrogen bond energy (hbond_sc), and 146 \nRamachandran preferences (rama) were applied for the soft-repulsive and hard-repulsive docking 147 \nphases (Appendix S6). 148 \nFor each docking run, either 20,000 poses or 100,000 poses were generated to assess if there are 149 \nstatistically significant differences in nearest native pose RMSD. In RosettaLigand, a ligand interface 150 \nis defined either by a representative ligand atom (a “neighbor atom”, defined as the geometric center 151 \nof mass by default), or all ligand atoms, relative to all ion channel Cβ atoms within a specified radius 152 \nfrom the ligand (commonly default to 6 or 7 Å). For RosettaLigand, two mutually exclusive ligand 153 \narea interface modes are available for scoring the pose interface: the ligand neighbor atom cutoffs 154 \nmode (add_nbr_radius=”True”), and the all ligand atom cutoffs mode 155 \n(all_atom_mode=”True”). In previous work, both modes were used (Moretti, Bender, Allison, 156 \nand Meiler, 2016; Smith and Meiler, 2020), thus, both modes were evaluated for any differences in 157 \nperformance. Four individual RosettaLigand docking sets were performed for each PDB structure, by 158 \ncombining different pose totals and ligand area interface modes, for performance comparison.  159 \nGALigandDock docking 160 \nAs described in the original study, docking was performed using the RosettaScripts’ 161 \nGALigandDock mover in DockFlex mode (Park et al, 2021; Appendix S8-S10). Replicate runs of 162 \nGALigandDock were performed in parallel for each structure evaluated. Each run consisted of 20 163 \ngenerations with a pool of 100 poses, where each generation updates the pool by total energy. By 164 \ndefault, GALigandDock outputs the top 20 structures from the final generation, however, for this 165 \nstudy, the entire pool of poses were used. For each docking run, 20,000 poses (1,000 runs) or 166 \n100,000 poses (5,000 runs) were output with a padding value of 2 Å, 4 Å, or 7 Å to test for statistical 167 \ndifferences in the lowest RMSD pose (RMSDMin). In total, six individual GALigandDock docking 168 \nsets were performed for each PDB structure, by combining different pose totals and padding sizes, 169 \nfor performance comparison.  170 \nPNear calculation 171 \nPNear is a quantitative metric to evaluate an energy funnel’s quality from a population of poses, 172 \nwith the lowest-scoring pose as the converged minima, or native state. PNear was calculated as 173 \npreviously defined (Bhardwaj, Mulligan, Bahl, et al., 2016), and applied as previously described for 174 \nligand docking (Smith and Meiler, 2020), with the native state being the ligand structure coordinates. 175 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n \n \nPNear ranges from 0 (protein will not converge to native state) to 1 (protein will always converge to 176 \nthe native state).  177 \n𝑃!\"#$ = ∑ 𝑒%$&'(!\n\"\n)\"!\n*+, 𝑒\n% -!\n.#/\n∑ 𝑒\n% -$\n.#/!\n0+,\n\t 185 \nThe numerator is the Boltzmann probability of an individual pose being near the native state, 178 \ngoverned by the “native-ness” parameter (λ), and the thermal energy, the product of the Boltzmann 179 \nconstant and absolute temperature (kBT). The denominator is the partition function of a canonical 180 \nensemble. For small-molecule docking, the root mean square deviation (RMSD) of the pose’s ligand 181 \ncoordinates relative to the structure’s ligand coordinates when the protein pose is superimposed to the 182 \nentire protein structure, or the receptor site is used. The energy scoring metric (Ei) is the interface 183 \nenergy: the energy solely composed of protein-ligand interactions.  184 \nA previous Rosetta small-molecule docking study (Park et al., 2021) categorized native-like 186 \nmodels at 1 Å or 2 Å RMSD without calculating PNear, while another study evaluating RosettaLigand 187 \nperformance calculated PNear using λ =1.5 Å, and kBT =0.62 (Meiler and Smith, 2020). Therefore, we 188 \ncalculated PNear with native-like models defined by λ=2.0 Å and kBT =0.62 Rosetta energy units 189 \n(REU), however we calculated PNear using all previously reported values for the parameters λ and kBT 190 \n(Tables S4.1-4.10). A specific PNear cutoff indicative for drug discovery pipelines has not been 191 \nestablished, hence we will refer to PNear ≥ 0.5 as a ‘first pass’ cutoff for this study when evaluating 192 \nenergy funnel convergence to the ligand structure coordinates. 193 \nStatistical tests 194 \nTests for normality, heteroscedasticity, and Pearson’s correlation between covariates were 195 \nperformed in Python using NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), and Pingouin 196 \n(Vallat, 2018) with a significance level (α) of 0.05. Population data consisted of the RMSDMin for 197 \neach docking set. Not all RMSDMin data for each docking set fit a normal distribution. Since 198 \nRMSDMin data should skew towards zero Å, a lognormal base 10 transformation was applied to all 199 \nRMSDMin population data when comparing across docking sets. Shapiro Wilk tests for normality and 200 \nLevene heteroscedasticity tests were performed to ensure transformed data were normal and of equal 201 \nvariance, respectively.  202 \nCovariates for both methods were the number of rotatable bonds of the ligand, the number of 203 \nligand heavy atoms, and the resolution of the structure. Statistical tests for RosettaLigand also 204 \nincluded the number of conformers generated, transformed with a lognormal base 10, as a covariate. 205 \nThe ligand molecular weight was initially included as a covariate but was discarded after identifying 206 \nstrong Pearson’s correlation with the number of ligand heavy atoms (Table S3, Figure S1.). 207 \nTo assess RMSDMin across all docking sets for a given method, a repeated measures factorial 208 \nanalysis of variance (ANOVA) with covariates were performed using IBM SPSS (version 29). For 209 \nRosettaLigand, the two factors were sample size and ligand area interface mode, while for 210 \nGALigandDock the two factors were sample size and padding value. Mauchly’s test for sphericity 211 \nwas performed for levels of padding (p = 0.63). 212 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\nResults 213 \nTable 1. Summary of ligand – ion channel structures docked, domain, and covariates included in the study. 214 \nPDB \nID Description \nIon channel \nmodulation / \ntherapeutic \nclass \nDomain \nTotal \nligand \nrotatable \nbonds \nTotal \nligand \nheavy \natoms \nLigand \nmolecular \nweight (Da) \nStructure \nresolution \n(Å) \nTotal \nRosettaLigand \nConformers \n5EK0 \nhuman NaV1.7-\nVSD4-NaVAb / \nGX-936 \nInhibitor / NA \nVoltage-\nsensing \ndomain \n8 39 575.6 3.53 2179 \n6J8G human NaV1.7 / \nSaxitoxin Blocker / NA Selectivity \nfilter 3 21 299.3 3.20 28 \n6J8I human NaV1.7 / \nTetrodotoxin Blocker / NA Selectivity \nfilter 1 22 319.3 3.20 2 \n6JP5 rabbit CaV1.1 / \nNifedipine \nBlocker / \nvasodilator, \nantihypertensive, \nantianginal \nCPC† 5 25 346.3 2.90 22 \n6JP8 \nrabbit CaV1.1 / \n(S)-(−)-Bay K \n8644 \nActivator / NA CPC 3 25 356.3 2.70 3 \n6JPA rabbit CaV1.1 / \n(S)-Verapamil \nBlocker / \nantiarrhythmic, \nvasodilator \nCPC 13 33 454.6 2.60 785 \n6JPB rabbit CaV1.1 / \nDiltiazem \nBlocker / \nantihypertensive, \nvasodilator \nCPC 7 29 414.5 2.90 45 \n6KZP human CaV3.1 / \nZ944 Blocker / NA CPC 6 26 383.9 3.10 1992 \n6U88 rat TRPV2 / \nCannabidiol \nActivator / \nanticonvulsants CPC 6 23 314.5 3.20 71 \n6UZ0 rat NaV1.5 / \nFlecainide \nBlocker / \nantiarrhythmic \nclass 1c \nCPC 7 28 414.3 3.24 1169 \n† CPC = Central pore cavity.215 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  \n 216 \nIn this study, we evaluated ligands that are suggested to modulate ion channel activity for a given 217 \nclinical effect, such as antiarrhythmic, anticonvulsant, and antihypertensive (Table 1). Our study also 218 \nincludes canonical ion channel blockers, such as tetrodotoxin and saxitoxin, small molecules whose 219 \nstructures have been resolved for drug discovery campaigns, such as GX-936, Z944, and Bay K 220 \n8544, and drugs currently approved for therapy, such as nifedipine, diltiazem, and flecainide (Figure 221 \n1).222 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  \n 223 \nFigure 1. Two-dimensional structures of ligands docked in this study. The depicted 224 \nstereochemistry reflects that resolved from structures. The protonation and bond order of saxitoxin 225 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n9 \nand tetrodotoxin from prior reported work was used (Hinman and DuBois, 2003; Thomas-Tran and 226 \nDu Bois, 2016). 227 \n 228 \nPerformance criteria for Rosetta ligand docking 229 \nA previous study has emphasized that two criteria must be satisfied for accurate modeling of 230 \nprotein-ligand interactions (Kaufmann and Meiler, 2012). First, the method must produce native-like 231 \nposes through sufficient sampling. For small molecule docking, a native-like ligand pose has a root 232 \nmean square deviation (RMSD) below a context-dependent predetermined value. Previously reported 233 \nRMSD ranges when assessing a native-like pose are 1.0 Å or 2.0 Å, with 2.0 Å being a common 234 \ncutoff for drug discovery pipelines (Maia et al., 2020; Park et al., 2021). When calculating RMSD, 235 \nthe protein pose is superimposed to the entire protein structure or the receptor site, and then the 236 \nligand coordinates of the pose are evaluated relative to the ligand coordinates of the structure.  237 \nSecond, the pose population must produce an energy funnel converging to the lowest-energy, 238 \nnative-like pose (Meiler and Baker, 2006). Rather than the lowest-energy pose, we use empirical 239 \nstructural evidence in this study to define the native pose as the ligand’s structural coordinates. An 240 \nenergy funnel is qualitatively evaluated using the interface energy between ligand and protein with 241 \nrespect to ligand RMSD (Meiler and Baker, 2006, Smith and Meiler, 2020). A quantitative metric of 242 \nfunnel convergence, PNear, has been utilized in Rosetta protein design (Guarav, Mulligan, Bahl, et al., 243 \n2016; Mulligan et al., 2021), and has been adopted for small molecule docking (Smith and Meiler, 244 \n2020). Due to the large number of poses, funnel quality is often assessed with a population subset 245 \nconsisting of the best-scoring poses (Combs et al., 2010; Lemmon, Kaufmann, and Meiler, 2012; 246 \nShim et al., 2019). 247 \nWith these criteria in mind, we evaluated the performance of RosettaLigand and GALigandDock 248 \nfor the following ion channel receptor sites: the voltage-sensing domain, the pore-forming domain, 249 \nand the central pore cavity (Table 1). We evaluated RosettaLigand using combinations of sample 250 \nsize (20,000 poses vs 100,000 poses) and ligand area interface mode (ligand neighbor atom vs all 251 \nligand atom). We evaluated GALigandDock using combinations of sample size (20,000 poses vs 252 \n100,000 poses) and padding of the sampling grid (2 Å, 4 Å, and 7 Å). Increasing the padding of the 253 \nsampling grid enables additional rotamer sampling around the receptor and increases translational 254 \nsampling of the ligand around the receptor site.255 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  \nTable 2. The minimum RMSD (Å) pose for each ligand-channel docked that is detailed in the 256 \nprimary text. 257 \nPDB ID Description \nRosettaLigand, all ligand \natom interface cutoff, \n100,000 poses \n GALigandDock, padding \n7 Å, 100,000 poses \n5EK0 human NaV1.7-VSD4-\nNaVAb / GX-936 0.91  0.65 \n6J8G human NaV1.7 / Saxitoxin 0.70  0.76 \n6J8I human NaV1.7 / \nTetrodotoxin 0.54  0.81 \n6JP5 rabbit CaV1.1 / Nifedipine 0.77  1.3 \n6JP8 rabbit CaV1.1 / (S)-(−)-\nBay K 8644 0.93  0.97 \n6JPA-1† rabbit CaV1.1 /  \n(S)-Verapamil 1.4  2.0 \n6JPA-2 rabbit CaV1.1 /  \n(S)-Verapamil 2.2  2.1 \n6JPB rabbit CaV1.1 / Diltiazem 2.5  1.1 \n6KZP human CaV3.1 / Z944 0.73  0.84 \n6U88 rat TRPV2 / Cannabidiol 1.0  0.90 \n6UZ0 rat NaV1.5 / Flecainide 1.2  0.94 \nAverage  1.2 ± 0.6  1.1 ± 0.5 \nMedian  0.93  0.94 \n† Verapamil was resolved with two binding orientations to the rabbit CaV1.1 central pore cavity in 258 \n6JPA and are named 6JPA-1 and 6JPA-2 in this study respectively. 259 \n 260 \n 261 \nCovariates potentially influencing RMSD 262 \nWe controlled for the following covariates that are dependent upon the ligand or structure used 263 \nfor docking. We speculated that the number of ligand rotatable bonds and the number of heavy atoms 264 \ncould influence RMSDMin by increasing the amount of internal sampling required during the docking 265 \nrun. We also speculated that a poorer, higher reported structural resolution could result in an 266 \nincreased overall RMSDMin. The total number of ligand heavy atoms was also used as a covariate for 267 \nRMSDMin, which is strongly correlated to the molecular weight of each ligand (Figure S1). Lastly, 268 \nfor RosettaLigand, the total number of conformers provided as input could affect the RMSDMin since 269 \nmore conformers could inherently require more sampling.  270 \n 271 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n11 \nRosettaLigand and GALigandDock generate native-like ion channel-ligand poses 272 \nThis study evaluated docking sets using population data consisting of the lowest RMSD pose 273 \n(RMSDMin) from each channel – ligand test case. We provide four docking sets for RosettaLigand 274 \nand six docking sets for GALigandDock, using combinations of the sample size with either the 275 \nRosettaLigand-specific ligand area interface calculation, or GALigandDock-specific padding value. 276 \nFor brevity, we will discuss the RosettaLigand docking set using a sample size of 100,000 poses and 277 \nthe all-ligand atoms interface mode, and the GALigandDock docking set using a sample size of 278 \n100,000 poses and padding value of 7 Å. These specific docking sets provide the greatest breadth of 279 \nreceptor sampling and scoring among our docking sets by appropriately encompassing the entirety of 280 \neach ion channel domain tested (Table 2). The results for other docking sets are provided in the 281 \nSupplemental Material (Tables S1.1-1.2). 282 \nFor RosettaLigand, the average RMSDMin across the docking set was 1.2 ± 0.6 Å (Table 2). 283 \nWhen comparing RosettaLigand docking sets and using repeated measures ANOVA with covariates, 284 \nan interaction from both factors (sample size and ligand area interface mode) did not result in 285 \nrejection of the null hypothesis of equivalent means from common logarithm transformed RMSDMin 286 \ndata (p=0.71). Likewise, individual interactions from sample size (p=0.38) and ligand area interface 287 \nmode (p=0.72) did not result in rejection of the null hypothesis of equivalent means. Each 288 \nRosettaLigand docking set produced comparable RMSDMin averages and standard deviations within a 289 \nfew sub-Angstroms (Table S1.1). Therefore, our results suggest that using either ligand area 290 \ninterface mode with 20,000 or 100,000 total poses could generate similar near-native models.  291 \nFor GALigandDock, the average RMSDMin across the docking set was 1.1 ± 0.5 Å (Table 2), 292 \nwhile all docking sets produced similar RMSDMin averages and standard deviations within a few sub-293 \nAngstroms (Table S1.2). Again, using repeated measures ANOVA with covariates, an interaction 294 \nfrom both factors (sample size and padding size) failed to result in a rejection of the null hypothesis 295 \nof equivalent means (p=0.91). Likewise, individual interactions from sample size (p=0.59) and 296 \npadding size (p=0.34) did not result in rejection of the null hypothesis of equivalent means. While 297 \nthere was no statistical advantage to using a specific padding value to achieve a lowered average 298 \nRMSDMin, we note that ion channel structures within the central pore cavity had greater sampling 299 \nwith a padding value of 7 Å, since a padding value of 5 Å did not encompass the entire pore. 300 \nTherefore, the appropriate padding size is context dependent and should be adjusted when using 301 \nGALigandDock to provide sufficient sampling grid space.302 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  \nTable 3. PNear of RMSDMin and interface energy for each ligand-channel docked that is detailed in the 303 \nprimary text, calculated with kBT=0.62 and λ=2.0. 304 \n  RosettaLigand, all ligand atom \ninterface cutoff, 100,000 poses  GALigandDock, padding 7 Å, \n100,000 poses \nPDB \nID Description Full \npopulation  Top 10% \ntotal_score  Full \npopulation  Top 10% \ntotal_score \n5EK0 \nhuman NaV1.7-\nVSD4-NaVAb / \nGX-936 \n0.22  0.26  0.83  0.83 \n6J8G human NaV1.7 \n/ Saxitoxin 1.2•10-2  9.3•10-3  1.6•10-2  1.6•10-2 \n6J8I human NaV1.7 \n/ Tetrodotoxin 0.12  8.8•10-2  0.27  0.27 \n6JP5 rabbit CaV1.1 / \nNifedipine 0.21  0.28  0.24  0.29 \n6JP8 \nrabbit CaV1.1 / \n(S)-(−)-Bay K \n8644 \n0.62  0.66  0.60  0.61 \n6JPA-\n1 \nrabbit CaV1.1 / \n(S)-Verapamil 4.7•10-3  4.3•10-3  2.4•10-3  2.8•10-3 \n6JPA-\n2 \nrabbit CaV1.1 / \n(S)-Verapamil 5.0•10-4  1.2•10-4  1.8•10-5  1.9•10-5 \n6JPB rabbit CaV1.1 / \nDiltiazem 6.8•10-4  6.4•10-4  2.8•10-6  7.8•10-6 \n6KZP human CaV3.1 \n/ Z944 0.50  0.62  1.2•10-4  1.1•10-4 \n6U88 rat TRPV2 / \nCannabidiol 8.3•10-2  4.3•10-2  5.0•10-3  3.8•10-3 \n6UZ0 rat NaV1.5 / \nFlecainide 0.14  0.12  5.9•10-6  2.9•10-5 \n 305 \n 306 \nRosettaLigand and GALigandDock energy funnels for ion channel-ligand docking 307 \nWe evaluated whether the entire population of generated poses and the top 10% of the lowest 308 \ntotal energy-scoring poses would achieve PNear values indicative of an energy funnel converging onto 309 \nthe ligand structural coordinates (Tables S4.1-4.10). Since interface energy describes how a ligand is 310 \ninteracting with the ion channel, PNear was calculated with the interface energy rather than the total 311 \nenergy. A specific PNear cutoff indicative for drug discovery pipelines has not been established. 312 \nHence, we will refer to PNear ≥ 0.5 as a ‘first pass’ cutoff for this study when evaluating energy funnel 313 \nconvergence to the ligand structure coordinates. For brevity, we will discuss PNear with a “native-314 \nness” parameter (λ) of 2.0, and a thermal energy parameter (kBT) of 0.62. Other PNear values 315 \nmatching parameter values reported in other work are provided in the Supplemental Material 316 \n(Bhardwaj, Mulligan, Bahl, et al., 2016; Smith and Meiler, 2020; Tables S4.1-4.10).  317 \nThe PNear of the full population and the top 10% of total energy-scoring poses were similar within 318 \nmethods, with the bulk of PNear values being in the thousandths or lower (Table 3). However, 319 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n13 \nRosettaLigand and GALigandDock were able to identify energy funnels for unique cases. For 320 \nRosettaLigand, the human CaV3.1- Z944 complex (PDB: 6KZP) achieved PNear ≥ 0.5, while for 321 \nGALigandDock, human NaV1.7-VSD4-NaVAb - GX-936 (PDB: 5EK0) achieved PNear ≥ 0.5. For 322 \nboth methods, rabbit CaV1.1-Bay K8644 (PDB: 6JP8) achieved PNear ≥ 0.5. 323 \nFollowing the standard of reporting a percentage of poses (Combs et al., 2010; Lemmon, 324 \nKaufmann, and Meiler, 2012; Shim et al., 2019), specific PNear values reported herein refer to the top 325 \n10% of total energy-scoring poses (Table 3).  326 \nNotably, there are visual distinctions in interface energy funnel plots between RosettaLigand and 327 \nGALigandDock (Figures S2-S12). Generally, RosettaLigand can sample more poses up to 2 Å 328 \nRMSD compared to GALigandDock, but with a less distinguishable energy funnel since some poses 329 \nscore greater than zero, indicating an unfavorable energy score. Since GALigandDock uses the 330 \nlowest total energy poses as input for future conformer generations, poses with interface energy 331 \ngreater than zero are infrequent.  332 \nLastly, when using the lowest interface energy pose from the full docking population as ranking 333 \ncriteria, the RMSDMin increases to a range between 1.2 Å and 10.8 Å (mean 4.5 ± 2.8 Å) with 334 \nRosettaLigand (Table S2.1), and a range between 0.83 Å and 14.5 Å (mean 6.4 ± 5.0 Å) with 335 \nGALigandDock (Table S2.5). This suggests that extracting the lowest energy pose is not a reliable 336 \nindicator of a near-native pose and does not necessarily reflect the native binding mode for ion 337 \nchannel-ligand docking. 338 \nLigand Docking into the voltage-sensing domain 339 \nThe only ion channel-ligand structure evaluated for ligand docking into the voltage-sensing 340 \ndomain (VSD) was GX-936 in complex with the VSD4 of the hNaV1.7-NaVAb chimera (PDB: 341 \n5EK0). The hNaV1.7 channel has been validated as a drug target for pain signaling, and aryl 342 \nsulfonamides have been reported as selective inhibitors. Specifically, GX-936 exhibits selectivity for 343 \nhNaV1.7 compared to other hNaV isoforms (McCormack et al., 2013; Ahuja et al., 2015; Nguyen and 344 \nYarov-Yarovoy, 2022).  345 \nAfter docking GX-936 to the hNaV1.7 VSD4, the RMSDMin poses were 0.91 Å using 346 \nRosettaLigand and 0.65 Å using GALigandDock (Figure 2, Table 2). The largest deviations from 347 \nthe native structure were observed at the peripheral pyrazole ring and the ethyl azetidine (Figure 2). 348 \nNotably, GX-936 has eight rotatable bonds, making it the second most flexible ligand used in this 349 \nstudy (Table 1; Figure 1). However, the sampled ligand poses were consistently below 2 Å RMSD 350 \nfor each docking run (Table S1.1-1.2). RosettaLigand was unable to identify the RMSDMin pose as 351 \nthe lowest interaction energy pose (Table S2.1), scoring the lowest energy pose 2.7 Å from the 352 \nstructure coordinates. This pose has the same binding orientation, but the pyrazole ring clashes with 353 \nthe structurally resolved lipid, while the ethyl azetidine withing the VSD is orientated up towards the 354 \nextracellular space, rather than down. Indeed, the lowest ten interface energy poses all possess these 355 \nfeatures (Figure S13.1). GALigandDock scored the lowest energy pose 0.83 Å from the structure 356 \ncoordinates (Table S2.5). Compared to RosettaLigand, the lowest ten interface energy poses 357 \nconverge well with the RMSDMin pose, with only one of the ten poses clashing with the lipid (Figure 358 \nS13.2). Further, GALigandDock was able to discriminate with good confidence an energy funnel 359 \n(PNear = 0.83; Table 3), whereas RosettaLigand was unable to distinguish a clear energy funnel (PNear 360 \n= 0.26; Table 3). 361 \n 362 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 14 \n363 \nFigure 2. GX-936 docking to hNa V1.7-NaVAb voltage sensing domain 4 (PDB 5EK0).  (A) 364 \nRosettaLigand (RMSDMin = 0.91 Å, PNear = 0.26) and (B) GALigandDock (RMSDMin = 0.65 Å, PNear 365 \n= 0.83). Left: transmembrane view . Middle: extracellular view . Right: ligand RMSD vs interface 366 \nenergy distribution of the top 10% of total energy poses. Grey dots: < 1 Å, blue dots: 1 -2 Å, yellow 367 \ndots: 2 -4 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown  in cyan, 368 \nGALigandDock in dark pink, native structure coordinates in dark grey and phospholipid head group 369 \nresolved from the structure in light grey. Non-carbon atoms match the Corey-Pauling-Koltun coloring 370 \nscheme. Black star indicates RMSDMin pose from the entire population. 371 \n 372 \nLigand Docking into the Pore-Forming Domain 373 \nTetrodotoxin (TTX) and saxitoxin (STX) are both guanidinium-based small molecules derived 374 \nfrom puffer fish and shellfish, and act as selective blockers of sodium channels by binding to the 375 \npore-forming domain (Hille, 2018; Shen et al., 2019). Being potent pore blockers for only some NaV 376 \nchannel isoforms, hNaV channel isoforms are classified in physiology as TTX-resistant (hNaV1.5, 377 \nhNaV1.8-1.9) and TTX-sensitive (hNaV1.1-1.4, hNaV1.6-1.7) (Stevens, Peigneur, and Tytgat, 2011). 378 \nLike aryl sulfonamides for VSDs, the discovery of TTX and STX binding to the pore-forming 379 \ndomain of NaV channels has spurred the design of similar blockers with greater selectivity for a 380 \nspecific hNaV isoform, usually hNaV1.7 for pain therapy (Hagen et al., 2017; Pajouhesh et al., 2020). 381 \nWe chose STX and TTX as test cases to evaluate ligand docking into the pore-forming domain of 382 \nhNaV1.7 (PDB: 6J8G and 6J8I, respectively; Shen et al., 2019). The protonation and bond order of 383 \nSTX and TTX from a previously reported work was used (Thomas-Tran and Du Bois, 2016). Based 384 \non their cage-like structures, STX has only 3 rotatable bonds and TTX has only 1 rotatable bond: 385 \nmaking them the most rigid ligands in this study (Table 1; Figure 1). In both cases, RosettaLigand 386 \ndocking resulted in RMSDMin values of 0.70 Å for STX and 0.54 Å for TTX, while GALigandDock 387 \nreported RMSDMin values of 0.76 Å for STX and 0.81 Å for TTX (Table 2; Figures 3 and 4). The 388 \nPNear values suggested little to no energy funnel convergence with STX (RosettaLigand PNear = 389 \n9.3•10-3 and GALigandDock PNear = 1.6•10-2), and TTX (RosettaLigand PNear = 7.8•10-2  and 390 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n15 \nGALigandDock PNear = 0.27) (Table 3). This lack of convergence is due to other energy minima 391 \noccurring within 3-6 Å RMSD (Figures 3 and 4).  392 \nFor STX, the lowest ten interface energy poses with RosettaLigand converged at 4.4 Å RMSD, 393 \nwith the carbamoyl group positioned deeper into the selectivity filter, while GALigandDock rendered 394 \na 3.8-6.4 Å RMSD, with the center of STX further away from the selectivity filter in various 395 \nrotations (Figure S14.1-2). For TTX, the lowest ten interface energy poses with RosettaLigand 396 \ncontained poses between 2.9-4.2 Å RMSD; five of ten poses converged to the structural coordinates 397 \nwith the amino group pointing away from the selectivity filter, while the other five poses converged 398 \nto the structural coordinates with the amino group pointing perpendicular to the selectivity filter path, 399 \n(Figure S15.1). With GALigandDock the poses were between 2.2-5.9 Å RMSD, with five of ten 400 \nposes having the TTX amino group pointing towards the pore and one set of P1/P2 helices, three of 401 \nten poses with the TTX amino group pointing toward the selectivity filter path, and one pose with the 402 \nTTX amino group pointing toward the extracellular space. (Figure S15.2). 403 \n404 \nFigure 3. Saxitoxin docking to hNaV1.7 selectivity filter (PDB 6J8G). (A) RosettaLigand 405 \n(RMSDMin = 0.70 Å, PNear = 9.3•10-3) and (B) GALigandDock (RMSDMin = 0.76 Å, PNear = 1.6•10-2). 406 \nLeft: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 407 \ndistribution of the top 10% of total energy poses. Grey dots: < 1 Å, blue dots: 1-2 Å, yellow dots: 2-4 408 \nÅ, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 409 \nin dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-410 \nPauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 411 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 16 \n412 \nFigure 4. Tetrodotoxin docking to hNaV1.7 selectivity filter (PDB 6J8I). (A) (RMSDMin = 0.54 Å, 413 \nPNear = 7.8•10-2) and (B) GALigandDock (RMSDMin = 0.81 Å, PNear = 0.27). Left: transmembrane view. 414 \nMiddle: extracellular view. Right: ligand RMSD vs interface energy distribution of the top 10% of total 415 \nenergy poses. Grey dots: < 1 Å, blue dots: 1-2 Å, yellow dots: 2-4 Å, green dots: > 4 Å. Carbon atoms 416 \nof RosettaLigand molecule are shown in cyan, GALigandDock in dark pink, and native structure 417 \ncoordinates in dark grey. Non-carbon atoms match the Corey-Pauling-Koltun coloring scheme. Black 418 \nstar indicates RMSDMin pose from the entire population. 419 \n 420 \nLigand docking into the central pore cavity 421 \nSeven test cases were used to evaluate ligand docking into the central pore cavity involving the 422 \nfollowing channels: four cases from rabbit CaV1.1, one from human CaV3.1, one from rat NaV1.5, 423 \nand one from TRPV2. The central pore cavity is a broad target, with small molecules reported to 424 \ntraverse through the gate or fenestration sites to reach their binding site and modulate channel activity 425 \nvia pore blockade or allosteric mechanism (Hille, 1977; Hockerman et al., 1997; Jiang et al., 2020). 426 \nDrugs targeting this region can act as vasodilators (dihydropyridines), antiarrhythmics 427 \n(benzothiazepines, phenylalkylamines, flecainide), antiepileptics (Z944, cannabidiol), or local 428 \nanesthetics (flecainide) (Pumroy et al., 2019; Zhao et al., 2019a; Zhao et al., 2019b; Jiang et al., 429 \n2020).  430 \nWe docked nifedipine (dihydropyridine channel blocker), Bay K 8644 (dihydropyridine channel 431 \nactivator), diltiazem (benzothiazepine channel blocker), and two conformations of verapamil 432 \n(phenylalkylamine channel blocker) into rabbit CaV1.1. The number of rotatable bonds was five, 433 \nthree, seven, and thirteen, respectively (Table 1; Figure 1). 434 \nNifedipine docking resulted in a RosettaLigand RMSDMin of 0.77 Å and a GALigandDock 435 \nRMSDMin of 1.3 Å (Figure 5, Table 2). The calculated PNear suggested little energy funnel 436 \nconvergence (RosettaLigand PNear = 0.28 and GALigandDock PNear = 0.29) (Table 3). This lack of 437 \nenergy convergence is exemplified by the interaction energy distributions containing multiple low 438 \nenergy minima 1-2 Å and 4-6 Å RMSD away from the ligand’s structural coordinates (Figure 5). 439 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n17 \nThe lowest ten interface energy poses with RosettaLigand contained poses between 1.1-4.4 Å 440 \nRMSD, with most poses being 1.7 Å or less; eight of ten poses converged to the structural 441 \ncoordinates with slight variation in rotamers, while two of the ten poses flipped the position of the 442 \ndihydropyridine and carbon ring relative to the structural coordinates (Figure S16.1). With 443 \nGALigandDock, the lowest ten interface energy poses were at 1.3 Å RMSD, in a similar position and 444 \norientation to the low RMSD conformations from RosettaLigand (Figure S16.2). 445 \n 446 \n447 \nFigure 5. Nifedipine docking to rabbit CaV1.1 central pore cavity (PDB 6JP5). (A) 448 \nRosettaLigand (RMSDMin = 0.77 Å, PNear = 0.28) and (B) GALigandDock (RMSDMin = 1.3 Å, PNear = 449 \n0.29). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs Interface 450 \nenergy distribution of the top 10% of total energy poses. Grey dots: < 1 Å, blue dots: 1-2 Å, yellow 451 \ndots: 2-4 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 452 \nGALigandDock in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match 453 \nthe Corey-Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire 454 \npopulation. 455 \n 456 \nBay K 8644 docking resulted in a RosettaLigand RMSDMin of 0.93 Å and GALigandDock 457 \nRMSDMin of 0.97 Å (Figure 6, Table 2). The calculated RosettaLigand PNear was 0.66 and the 458 \nGALigandDock PNear was 0.61 (Table 3). The PNear values, paired with the interaction energy 459 \ndistribution data, indicate well-converged energy funnels. Further, the lowest interaction energy 460 \nposes with RosettaLigand was within 0.3 Å of the RMSDMin pose (Table S2.1) and for 461 \nGALigandDock was within 0.2 Å of the RMSDMin pose (Table S2.5). With RosettaLigand, the 462 \nlowest ten interface energy poses were between 1.2-1.3 Å RMSD, with all ten poses converged to the 463 \nstructural coordinates with slight variation in rotamers (Figure S17.1). With GALigandDock, the 464 \nlowest ten interface energy poses were at 1.1 Å or 4.6 Å RMSD. Eight of ten poses had 1.1 Å RMSD 465 \nwith slight deviation in position to structural coordinates, while two of ten poses had 4.6 Å RMSD 466 \nwith the dihydropyridine in the correct position, but the carbon ring flipped 180 degrees such that the 467 \ntrifluoromethyl group was pointed in the opposite direction (Figure S17.2). 468 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 18 \n469 \nFigure 6. Bay K 8644 docking to rabbit CaV1.1 central pore cavity (PDB 6JP8). (A) 470 \nRosettaLigand (RMSDMin = 0.93 Å, PNear = 0.66) and (B) GALigandDock (RMSDMin = 0.97 Å, PNear 471 \n= 0.61). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface 472 \nenergy distribution of the top 10% of total energy poses. Grey dots: < 1 Å, blue dots: 1-2 Å, yellow 473 \ndots: 2-4 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 474 \nGALigandDock in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match 475 \nthe Corey-Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire 476 \npopulation. 477 \nDiltiazem docking resulted in a RosettaLigand RMSDMin of 2.5 Å and a GALigandDock 478 \nRMSDMin of 1.1 Å (Figure 7, Table 2). The calculated PNear suggested no energy funnel convergence 479 \nfor either method (RosettaLigand PNear = 6.4•10-4 and GALigandDock PNear= 7.8•10-6) (Table 3). 480 \nThis lack of energy funnel convergence is exemplified by the interaction energy distributions 481 \ncontaining local minima 5-14 Å RMSD away from the ligand’s structural coordinates (Figure 7). 482 \nWith RosettaLigand, the lowest ten interface energy poses were between 5.2-7.6 Å RMSD, with all 483 \nten poses at a similar channel depth at the pore center, but with no convergence in local position or 484 \nrotamer placement. (Figure S20.1). With GALigandDock, the lowest ten interface energy poses were 485 \nbetween 10.5-14.5 Å RMSD. Rather than being positioned in the pore center, all ten poses were in a 486 \nsimilar channel depth to the structural coordinates but positioned in the fenestration with no 487 \nconvergence in local position or rotamer placement (Figure S20.2). 488 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n19 \n489 \nFigure 7. Diltiazem docking to rabbit CaV1.1 central pore cavity (PDB 6JPB). (A) RosettaLigand 490 \n(RMSDMin = 2.5 Å, PNear = 6.4•10-4) and (B) GALigandDock (RMSDMin = 1.1 Å, PNear = 7.8•10-6). 491 \nLeft: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 492 \ndistribution of the top 10% of total energy poses. Yellow dots: 2-4 Å, green dots: > 4 Å. Carbon 493 \natoms of RosettaLigand molecule are shown in cyan, GALigandDock in dark pink, and native 494 \nstructure coordinates in dark grey. Non-carbon atoms match the Corey-Pauling-Koltun coloring 495 \nscheme. Black star indicates RMSDMin pose from the entire population. 496 \n 497 \nVerapamil was previously resolved in a complex with rabbit CaV1.1 in two binding modes with 498 \nflipped orientations (Zhao et al., 2019a). We evaluated if RosettaLigand and GALigandDock 499 \nsampling favored a particular orientation. For the first orientation, docking resulted in a 500 \nRosettaLigand RMSDMin of 1.4 Å and a GALigandDock RMSDMin of 2.0 Å (Figure 8, Table 2). The 501 \ncalculated PNear suggested no energy funnel convergence (RosettaLigand PNear = 4.3•10-3, 502 \nGALigandDock PNear = 2.8•10-3) (Table 3). This nonexistent energy funnel is evident by local 503 \ninteraction energy minima beyond 4 Å RMSD from the ligand structure coordinates (Figure 8). With 504 \nRosettaLigand, the lowest ten interface energy poses were between 4.1-9.7 Å RMSD, with all ten 505 \nposes at a similar channel depth, and some poses converging in local position and rotamer placement 506 \nat 5.7 Å and 9.2 Å RMSD (Figure S18.1). With GALigandDock, the lowest ten interface energy 507 \nposes were between 3.8-11.9 Å RMSD. All ten poses were at different channel depths, positioned in 508 \nthe pore center or at the pore center – fenestration interface, and did not converge in local position or 509 \nrotamer placement (Figure S18.2). 510 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 20 \n511 \nFigure 8. Verapamil docking in its first orientation to rabbit CaV1.1 central pore cavity (PDB 512 \n6JPA). (A) RosettaLigand (RMSDMin = 1.4 Å, PNear = 4.3•10-3) and (B) GALigandDock (RMSDMin = 513 \n2.0 Å, PNear = 2.8•10-3). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD 514 \nvs interface energy distribution of the top 10% of total energy poses. Blue dots: 1-2 Å, yellow dots: 515 \n2-4 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 516 \nGALigandDock in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match 517 \nthe Corey-Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire 518 \npopulation. 519 \n 520 \nFor the second orientation, docking resulted in a RosettaLigand RMSDMin of 2.2 Å and a 521 \nGALigandDock RMSDMin of 2.1 Å (Figure 9, Table 2). The calculated PNear values suggest no 522 \nenergy funnel convergence (RosettaLigand PNear = 1.2•10-4, GALigandDock PNear = 1.9•10-5) (Table 523 \n3). Similar to the first orientation, the interaction energy distribution for the second orientation yields 524 \nlocal energy minima greater than 4 Å RMSD from the ligand structure coordinates (Figure 9). With 525 \nRosettaLigand, the lowest ten interface energy poses were between 8.5-10.8 Å RMSD, with all ten 526 \nposes at a similar channel depth, and some poses converging in local position, similar to the docking 527 \nset for the first structural orientation of verapamil (Figure S19.1). With GALigandDock, the lowest 528 \nten interface energy poses were between 6.5-10.8 Å RMSD. All ten poses were at different channel 529 \ndepths, positioned in the pore center or at the pore center – fenestration interface, and did not 530 \nconverge in local position or rotamer placement, similar to the docking set for the first structural 531 \norientation of verapamil (Figure S19.2). 532 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n21 \n533 \nFigure 9. Verapamil docking in its second orientation to rabbit CaV1.1 central pore cavity 534 \n(PDB 6JPA). (A) RosettaLigand (RMSDMin = 2.2 Å, PNear = 1.2•10-4) and (B) GALigandDock 535 \n(RMSDMin = 2.1 Å, PNear = 1.9•10-5). Left: transmembrane view. Middle: extracellular view. Right: 536 \nligand RMSD vs interface energy distribution of the top 10% of total energy poses. Yellow dots: 2-4 537 \nÅ, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 538 \nin dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-539 \nPauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 540 \n 541 \nThe small molecule Z944 (channel blocker) has 6 rotatable bonds and is resolved bound to 542 \nhuman CaV3.1. Docking with RosettaLigand and GALigandDock both resulted in a RMSDMin of 0.73 543 \nÅ and 0.84 Å, respectively (Figure 10, Table 2). Interestingly, the calculated PNear suggests good 544 \nenergy funnel convergence for RosettaLigand (PNear = 0.62), but no energy funnel convergence for 545 \nGALigandDock (PNear = 1.1•10-4) (Table 3). For RosettaLigand, the interaction energy produced an 546 \nenergy funnel between 1-2 Å, while for GALigandDock produced an energy funnel near 10-15 Å 547 \n(Figure 10). With RosettaLigand, the lowest ten interface energy poses were between 0.90-1.7 Å 548 \nRMSD, with all ten poses converged to the structural coordinates with slight variation in rotamers. 549 \n(Figure S21.1). With GALigandDock, the lowest ten interface energy poses were between 10.6-12.4 550 \nÅ RMSD, with four sets of unique poses identified. The first set, containing three of the ten poses, is 551 \nin a slightly lower channel depth and in a similar orientation to the structural coordinates, but rotated 552 \napproximately 180 degrees about the pore center. The second set is one pose and is in a similar 553 \nlocation as the first set but is rotated such that the phenyl ring and tertiary butylamine positions are 554 \nflipped. The third group, containing four poses, is similar to the second set but rotated approximately 555 \n180 degrees about the pore center. The last group, containing two poses, is in a similar channel depth 556 \nand location to the structural coordinates, but the phenyl ring and the tertiary butylamine positions 557 \nare flipped (Figure S21.2). 558 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 22 \n559 \nFigure 10. Z944 docking to human CaV3.1 central pore cavity (6KZP). (A) RosettaLigand 560 \n(RMSDMin = 0.73 Å, PNear = 0.62) and (B) GALigandDock (RMSDMin = 0.84 Å, PNear = 1.1•10-4). 561 \nLeft: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 562 \ndistribution of the top 10% of total energy poses. Grey dots: < 1 Å, blue dots: 1-2 Å, yellow dots: 2-4 563 \nÅ, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 564 \nin dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-565 \nPauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 566 \n 567 \nFlecainide (channel blocker) has 7 rotatable bonds and is resolved bound to rat NaV1.5. 568 \nRosettaLigand docking resulted in a RMSDMin of 1.2 Å, while GALigandDock resulted in a 569 \nRMSDMin of 0.94 Å (Figure 11, Table 2). The calculated PNear values suggest no energy funnel 570 \nconvergence with either method (RosettaLigand PNear = 0.12, GALigandDock PNear = 2.9•10-5) 571 \n(Table 3). The interaction energy distribution from RosettaLigand produced an energy minimum 572 \nbetween 2-4 Å and 7-10 Å, while from GALigandDock produced energy minima near 10-12 Å 573 \n(Figure 11). With RosettaLigand, the lowest ten interface energy poses were between 2.1-4.0 Å 574 \nRMSD, with all ten poses in the same channel depth as the structural coordinates, but with a tilt such 575 \nthat the trifluoromethyl groups are level in channel depth, rather than the trifluoromethyl group 576 \nfacing the pore being lower in depth. In one pose, the piperidin-2-yl-methylamine was positioned 577 \nlower into the channel than the rest of the ligand and the structural coordinates, while the rest of the 578 \nligand pose was in a similar position to other poses (Figure S23.1). With GALigandDock, the lowest 579 \nten interface energy poses were between 9.7-11.9 Å RMSD, with all ten poses positioned at the 580 \nfenestration in the same channel depth as the structural coordinates, but oriented outside of the pore 581 \nwith one of the trifluoromethyl groups pointing towards the pore center but positioned at the exterior 582 \nof the channel fenestration (Figure S23.2). 583 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n23 \n584 \nFigure 11. Flecainide docking to rat NaV1.5 central pore cavity (PDB 6UZ0). (A) RosettaLigand 585 \n(RMSDMin = 1.2 Å, PNear = 0.12) and (B) GALigandDock (RMSDMin = 0.94 Å, PNear = 2.9•10-5). Left: 586 \ntransmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 587 \ndistribution of the top 10% of total energy poses. Blue dots: 1-2 Å, yellow dots: 2-4 Å, green dots: > 588 \n4 Å. Carbon atoms of RosettaLigand molecule is indicated in cyan, GALigandDock in dark pink, and 589 \nnative structure coordinates in dark grey. Non-carbon atoms match the Corey-Pauling-Koltun 590 \ncoloring scheme. Black star indicates RMSDMin pose from the entire population. 591 \n 592 \nThe transient receptor potential vanilloid (TRPV) subfamily of cation channels broadly play roles 593 \nin thermo-sensation and thermoregulation in response to heat (>53 °C), notably in noxious heat 594 \nsensing (Nilius, 2007; Gees et al., 2012). TRPV2 channels are widely expressed, with identification 595 \nin dorsal root ganglion neurons, as well as brain, heart, and smooth muscle tissue, among others 596 \n(Gees et al., 2012). The TRPV2 channel has been implicated in thermal pain sensing, muscular 597 \ndystrophy, and cardiomyopathy, among other diseases (Nilius, 2007; Gees et al., 2012). Structurally, 598 \nTRPV2 channels contain six transmembrane spanning domains and commonly assemble as a 599 \nhomomer with four identical subunits (Zheng and Trudeau, 2015). They resemble canonical voltage-600 \ngated ion channels with a pore domain and a weak voltage sensor domain, while their gating is 601 \nregulated by heat and a diverse set of agonists such as 2-aminoethoxydiphenyl borate and cannabidiol 602 \n(Gees et al., 2012; Zheng and Trudeau, 2015; Pumroy et al., 2019). 603 \nCannabidiol (channel activator) has 6 rotatable bonds and is resolved bound to rat TRPV2. 604 \nDocking with RosettaLigand resulted in a RMSDMin of 1.0 Å, while GALigandDock resulted in a 605 \nRMSDMin of 0.90 Å (Figure 12, Table 2). The calculated PNear values suggest no energy funnel 606 \nconvergence for both methods (RosettaLigand PNear  = 4.3•10-2, GALigandDock PNear = 3.8•10-3) 607 \n(Table 3). For RosettaLigand, the interaction energy distribution did not demonstrate a clear 608 \nminimum, with minima seen between 2-8 Å (Figure 12). For GALigandDock, there is a clear 609 \nminimum between 6-7 Å (Figure 12). With RosettaLigand, the lowest ten interface energy poses 610 \nwere between 2.8-7.0 Å RMSD, with two groups of conformations in the same channel depth and 611 \ninterfaced with the S5 and S6 helices of adjacent TRPV2 monomers like the structural coordinates. In 612 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 24 \none group, cannabidiol is positioned parallel to the S6 helical segment, with the cyclohexene at the 613 \nlowest channel depth. In the second group, the poses are positioned in the same orientation as the 614 \nstructural coordinates, with the pentyl chain facing away from the pore center, but the overall 615 \nposition of cannabidiol is slightly elevated relative to the structural coordinates (Figure S22.1). With 616 \nGALigandDock, the lowest ten interface energy poses were either 6.4 Å or 6.7 Å RMSD, with all ten 617 \nposes positioned at the fenestration in the same channel depth and interfacing with the S5 and S6 618 \nhelices of adjacent TRPV2 monomers like the structural coordinates, but with the ligand rotated 619 \naround the cyclohexene plane such that the pentyl chain is oriented parallel to the S6 helix, rather 620 \nthan pointing outward (Figure S22.2). 621 \n622 \nFigure 12. Cannabidiol docking to rat TRPV2 central pore cavity (PDB 6U88). (A) 623 \nRosettaLigand (RMSDMin = 1.0 Å, PNear = 4.3•10-2) and (B) GALigandDock (RMSDMin = 0.90 Å, 624 \nPNear = 3.8•10-3). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs 625 \ninterface energy distribution of the top 10% of total energy poses. Blue dots: 1-2 Å, yellow dots: 2-4 626 \nÅ, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 627 \nin dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-628 \nPauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 629 \n 630 \nDiscussion 631 \nPrevious studies have underscored the importance of ligand docking methods for generating ion 632 \nchannel structure-based hypotheses (Yang et al., 2013; Shim et al., 2019). Furthermore, ligand 633 \ndocking methods when combined with high-resolution structures can aid in rational drug design 634 \n(Wang et al., 2007; Liu et al., 2018; Wulff et al., 2019; Maia et al., 2020). Notably, RosettaLigand 635 \nhas been extensively used to predict the molecular mechanisms of ligand – ion channel interactions 636 \n(Yang et al., 2013; Yang et al., 2015; Nguyen et al., 2017; Nguyen et al., 2019; Shim et al., 2019; 637 \nCraig 2nd et al., 2022; Vu et al., 2020; Maly et al., 2022; Pumroy et al., 2022). While GALigandDock 638 \nhas not yet been tested on ion channel structure-function relationships, it utilizes a new generalized 639 \nenergy function tailored for small molecules while sampling ligand conformations with a genetic 640 \nalgorithm, making it an attractive complement to RosettaLigand. 641 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n25 \nIndeed, this study demonstrates that high-resolution structures paired with RosettaLigand and 642 \nGALigandDock can be useful tools for formulating structural hypotheses and predicting binding 643 \nmodes for drug discovery. Notably, both Rosetta methods can produce ligand poses near the native 644 \nligand structure coordinates. Using a standard 2 Å RMSDMin cutoff (Maia et al., 2020; Park et al., 645 \n2021), the RosettaLigand docking set yielded an average RMSDMin of 1.2 ± 0.6 Å, and the 646 \nGALigandDock docking set yielded an average RMSDMin of 1.1 ± 0.5 Å (Table 2). However, the 647 \nability to discriminate the RMSDMin pose remains challenging, and highlights important practical 648 \nconsiderations when applying Rosetta ligand docking methods to a chosen ion channel target. 649 \nWhen performing Rosetta ligand docking, the features of the ligand conformer library, the size of 650 \nthe receptor site, and prior knowledge of critical functional residues, need to be considered to 651 \ndetermine the appropriate amount of pose sampling. For our study, we generated either 20,000 poses 652 \nor 100,000 poses per docking run. While we did notice a statistically significant difference in 653 \nRMSDMin between 20,000 and 100,000 total poses, the difference was within a few sub-Angstroms. 654 \nFurthermore, another RosettaLigand benchmarking study of the CASF-2016 dataset generated 1,000 655 \ntotal poses and was able to sufficiently sample a RMSDMin below 2 Å (Su et al., 2019; Smith and 656 \nMeiler, 2020). Thus, for ion channel docking, we suggest sampling in the ones to tens of thousands, 657 \nespecially when docking in large receptor sites like the pore-forming domain.  658 \nOverall, our docking data was able to achieve a PNear greater than 0.5 for two test cases for each 659 \ndocking method when using the top 10% of total-scoring poses. Both methods achieved PNear greater 660 \nthan 0.5 for Bay K8644 to the rabbit CaV1.1 central pore cavity, while RosettaLigand achieved PNear 661 \ngreater than 0.5 for Z944 to the human CaV3.1 central pore cavity, and GALigandDock achieved 662 \nPNear greater than 0.5 for GX-936 docking to hNaV1.7-NaVAb voltage sensing domain four (Table 3). 663 \nThe reasons for a low PNear value in most test cases are possibly due to 1) improper scoring by 664 \nRosetta score functions to discriminate native-like poses by energy, 2) multiple favorable ligand 665 \nbinding states in the receptor site that have not been structurally resolved, and/or 3) insufficient pose 666 \nsampling. 667 \nCurrently, we suggest that Rosetta score functions are unable to sufficiently score near-native 668 \nposes accurately in ion channel docking. For each docking run, comparing the RMSDMin pose and the 669 \nlowest interface energy pose indicates that the RMSDMin pose is under-scored (Tables S2.1-2.5). 670 \nConversely, a previous RosettaLigand study using the CASF-2016 dataset (containing 285 crystal 671 \nstructures of protein – ligand complexes with an overall resolution < 2.5 Å and an R-factor < 0.25) 672 \nand 1,000 total poses per docking run was able to frequently achieve PNear values between 0.8 and 1.0 673 \n(Su et al., 2019; Smith and Meiler, 2020), suggesting that the Rosetta score functions can be utilized 674 \nfor other docking studies, but should be verified for accuracy with test cases. Furthermore, the 675 \nCASF-2016 dataset does not contain ion channels, while the sample size per docking run is 676 \nsufficiently lower in the RosettaLigand docking study using CASF-2016 than those in our study, 677 \nsuggesting that sampling is not the predominant issue, but rather favorable scoring of near-native 678 \nposes is the primary issue.  679 \nFor the voltage-sensing domains, GX-936 in complex with the VSD4 of hNaV1.7-NaVAb chimera 680 \n(PDB: 5EK0) was the only small-molecule docked yet that was consistently near or below 1 Å 681 \nRMSDMin for each docking set regardless of the method (Tables S1.1-1.2). This may be due to the 682 \nVSD4 receptor site being the narrowest binding pocket tested, thereby limiting the number of binding 683 \nconfigurations. This suggests that VSDs are generally well-suited for Rosetta ligand docking since 684 \nthe receptor is constrained to allow shape-complementary between ligand and protein while reducing 685 \nthe required sampling compared to pore-forming receptor sites. Notably, the Rosetta ligand docking 686 \nmethods employed did not use any implicit membrane parameters, while GX-936 in a biologically 687 \nrealistic context is partially exposed to a lipid head group (Ahuja et al., 2015). Thus, in this case, 688 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 26 \ndocking without an implicit membrane energy function was still able to achieve the experimental 689 \nstructurally resolved position, however, artifactual, low interface energy poses where GX-936 would 690 \nbe overlapping the phospholipid space are present (Figure S13.1-13.2). Further, GALigandDock was 691 \nconsistently able to achieve a PNear > 0.5 when using a padding of 2.0 Å, 4.0 Å, or 7.0 Å, suggesting 692 \nthat ligands docked to VSDs could be screened by energy funnel convergence (Tables S4.5-4.10). 693 \nMore small molecules structurally resolved and docked to the VSD are needed to validate this 694 \nobservation. 695 \nFor the pore-forming domain, TTX and STX were docked to hNaV1.7, as they are canonical 696 \nsmall-molecule channel blockers used when studying the NaV family. TTX and STX had the fewest 697 \nrotatable bonds in this study, suggesting that little conformer sampling would be needed to achieve a 698 \nnear-native pose. Both TTX and STX docking achieved sub-Angstrom RMSDMin, except STX 699 \ndocking with GALigandDock at a padding of 7 Å (RMSDMin =2.3 Å, padding 7 Å/100,000 poses vs 700 \nRMSDMin = 1.0 Å, padding 4 Å/100,000 poses; Table S1.2). However, all docking methods 701 \nemployed were unable to achieve a PNear > 0.5, suggesting that 1) STX and TTX could have 702 \nalternative binding conformations, 2) the docking methods employed have wrongfully scored 703 \nalternative low-energy conformations, or 3) that inherent loop flexibility in the pore-forming domain 704 \nis a requirement for docking an energy-optimized, induced-fit conformation. It should be noted that 705 \nthe hNaV1.7 selectivity filter is lined with polar residues that could contribute to a vast hydrogen 706 \nbonding network, in addition to water hydrating the selectivity filter and the filter opening. Thus, it is 707 \nunclear if the lowest energy poses, with potential salt bridge interactions, are possible alternate 708 \nbinding modes. It is thus possible that the sum of hydrogen bonding interactions and water-bridging 709 \neffects could bias the energetic potential to a certain state, such as the one structurally resolved. Thus, 710 \nfurther experimental characterization is needed to test these structural hypotheses. 711 \n 712 \nFor the central pore cavity, while there should be little to no lipid interactions (excluding the 713 \nfenestration regions), docking of the central pore cavity produced the most pose variability, likely 714 \ndue to it being the widest receptor site compared to the pore-forming domain and voltage sensing 715 \ndomain. For example, the second orientation of verapamil positioned primarily in the central cavity 716 \nof rabbit CaV1.1 achieved a RMSDMin of 1.4 Å for RosettaLigand and 2.0 Å for GALigandDock 717 \n(Table 2). However, Z944 bound to hCaV3.1, with the wide aromatic end of Z944 in the narrow 718 \nfenestration and the narrow amide end in the wide cavity achieved a RMSDMin of 0.73 Å for 719 \nRosettaLigand and 0.84 Å for GALigandDock (Table 2). Further, RosettaLigand was able to achieve 720 \na PNear value greater than 0.5 in some docking conditions for Z944 (Table 3; Table S4.1-4.4), 721 \nsuggesting that the Rosetta docking methods could prove useful when docking similar ligands that 722 \nbridge between the fenestration and central pore, and target the narrow fenestration with an aromatic 723 \nmoiety. 724 \nSmall molecules docked in the central cavity were bound to the central pore (6JPA, orientation 2, 725 \n6JPB), the channel fenestration region (6JP5, 6JP8, 6U88), or both regions (6JPA orienation1, 6KZP, 726 \n6UZ0). It appears that those bound in the channel fenestration produced a lower RMSDMin, those 727 \nprimarily bound in the pore center produced the largest RMSDMin, while those bound to both regions 728 \nproduced intermediate RMSDMin values, with Z944 docked to human CaV3.1 central pore cavity 729 \n(6KZP) being an exception. The same trend is roughly observed for PNear, where fenestration-only 730 \ncases have PNear values orders of magnitude greater than pore center cases (Table S5). Due to the 731 \nlimited number of cases for each sub-classification, further studies will need to be performed to 732 \nassess a correlation.  733 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n27 \nIt appears that molecules with predominantly planar aromatic rings and space-filling, “bulky” 734 \nstructures targeting the fenestration can be scored and ranked effectively; both Bay K 8644 (6JP8) 735 \nand Z944 (6KZP) possess these features in contrast to other small molecules, that possess aromatic or 736 \naliphatic rings but are generally more linear and “floppy” targeting the fenestration, like flecainide 737 \n(6UZ0) and the first orientation of verapamil (6JPA). (Table S5). This general trend in increased 738 \nRMSDMin for non-aromatic containing small molecules when the receptor site is solely the central 739 \npore could be due to 1) greater ligand flexibility within a larger area, which requires increased 740 \nsampling, 2) a bulkier ligand having fewer conformations and orientations relative to the channel, 3) 741 \nthe absence of explicit water molecules and ions that would crowd the pore or directly bind to the 742 \nligand (Tikhonov and Zhorov, 2008). Furthermore, if the ligand is expected to bind to a central cavity 743 \nmotif present on multiple subunits, then post-hoc tools implementing symmetry to discriminate 744 \nunique binding modes will be necessary to calculate an appropriate adjusted RMSDMin. 745 \n 746 \nConclusions 747 \nIn this study, we aimed to assess if 1) the RosettaLigand and GALigandDock molecular docking 748 \nmethodologies can recapitulate structurally resolved ion channel – small-molecule binding 749 \norientations with known pharmacological significance, 2) if their scoring functions can be used to 750 \naccurately rank small molecule binding orientation for the purpose of blindly screening small 751 \nmolecules, and 3) if there are practical considerations when docking to specific domains of ion 752 \nchannels. With 2.0 Å RMSDMin as a performance cutoff (Maia et al., 2020; Park et al., 2021), our 753 \nresults demonstrate that both RosettaLigand and GALigandDock can frequently sample the 754 \nexperimentally resolved ligand binding mode with less than 2.0 Å RMSD. However, when using an 755 \nestimate of the Boltzmann probability for energy “funnel-likeness” (PNear) as a scoring function 756 \nassessment, we currently perceive Rosetta score functions as unable to sufficiently score near-native 757 \nposes accurately in ion channel docking; from this study, small molecules targeting voltage-sensing 758 \ndomains and bulky small molecules primarily composed of aromatic rings targeting fenestration 759 \nregions appear to be most suited for score-based ranking. Thus, when performing an ion channel 760 \nvirtual drug discovery campaign, special consideration should be given to sufficient pose sampling to 761 \naccount for multiple rotameric and conformational states, to identify the size of the sampling required 762 \nfor sufficient interface scoring of the receptor site, to identify the appropriate state of the ion channel, 763 \nto identify inherent channel flexibility that could influence ligand binding, and to potentially identify 764 \nspecific chemical functional groups known experimentally to influence binding to the target when 765 \nselecting candidate conformations. 766 \n 767 \nAcknowledgements 768 \nWe thank members of the VY-Y, FD, and HW laboratories for helpful discussions. We also thank 769 \nIgor Vorobyov for helpful discussions and Jerome Manera for help with generating 2d structures of 770 \nligands. The work in VY-Y lab was supported by NIH grants R61NS127285, R01HL128537, 771 \nR01HL159304, R01GM132110, R56NS9706, and R01NS128180. BJH was supported by NIH F31 772 \nPredoctoral Fellowship F31NS124337 and NIH T32 Predoctoral Fellowship GM007377. We are 773 \ngrateful to OpenEye Scientific for OpenEye Academic License. 774 \n 775 \n  776 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 28 \nReferences 777 \nAhuja, S., Mukund, S., Deng, L., Khakh, K., Chang, E., Ho, H., et al., and Payandeh, J. (2015) 778 \nStructural basis of NaV1.7 inhibition by an isoform-selective small-molecule antagonist. Science. 779 \n350(6267), aac5464. https://doi.org/10.1126/science.aac5464 780 \n 781 \nAlexander, S. P., Kelly, E., Mathie, A., Peters, J. A., Veale, E. L., Armstrong, J. F., et al., and Verri, 782 \nT. (2021). The concise guide to pharmacology 2021/22: transporters. Br. J. Pharmacol. 178, S412–783 \nS513. https://doi.org/10.1111/bph.15543 784 \n 785 \nAlford, R. F., Leaver-Fay, A., Jeliazkov, J. R., O'Meara, M. J., DiMaio, F. P., Park, H., Shapovalov, 786 \nM. V., Renfrew, P. D., Mulligan, V. K., Kappel, K., Labonte, J. W., Pacella, M. S., Bonneau, R., 787 \nBradley, P., Dunbrack, R. L., Jr, Das, R., Baker, D., Kuhlman, B., Kortemme, T., and Gray, J. J. 788 \n(2017). The Rosetta all-atom energy function for macromolecular modeling and design. J. Chem. 789 \nTheory Comput., 13(6), 3031–3048. https://doi.org/10.1021/acs.jctc.7b00125 790 \n 791 \nBagal, S. K., Brown, A. D., Cox, P. J., Omoto, K., Owen, R. M., Pryde, D. C., Sidders, B., Skerratt, 792 \nS. E., Stevens, E. B., Storer, R. I., and Swain, N. A. (2013). Ion channels as therapeutic targets: a 793 \ndrug discovery perspective. J. Med. Chem., 56(3), 593–624. https://doi.org/10.1021/jm3011433 794 \n 795 \nBhardwaj, G., Mulligan, V. K., Bahl, C. D., Gilmore, J. M., Harvey, P. J., Cheneval, O., et al., and 796 \nBaker, D. (2016). Accurate de novo design of hyperstable constrained peptides. Nature, 538(7625), 797 \n329–335. https://doi.org/10.1038/nature19791 798 \n 799 \nBennett, D. L., Clark, A. J., Huang, J., Waxman, S. G., and Dib-Hajj, S. D. (2019). The Role of 800 \nVoltage-Gated Sodium Channels in Pain Signaling. Physiological reviews, 99(2), 1079–1151. 801 \nhttps://doi.org/10.1152/physrev.00052.2017 802 \n 803 \nBerman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., Shindyalov, I. N., 804 \nBourne, P. E. (2000). The Protein Data Bank. Nucleic Acids Res. 28(1), 235-242. 805 \nhttps://doi.org/10.1093/nar/28.1.235 806 \n 807 \nCase, D. A., Aktulga, H. M., Belfon, K., Ben-Shalom, I. Y., Brozell, S. R., Cerutti, D. S., et al., and 808 \nKollman, P. A. (2021). Amber 2021. University of California, San Francisco. 809 \n 810 \nCatterall, W. A. (1995) Structure and function of voltage-gated ion channels. Annu. Rev. Biochem 64, 811 \n493–531. https://doi.org/10.1146/annurev.bi.64.070195.002425 812 \n 813 \nCatterall, W. A. (2011). Voltage-gated calcium channels. Cold Spring Harb. Perspect. Biol. 3(8), 814 \na003947. https://doi.org/10.1101/cshperspect.a003947 815 \n 816 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n29 \nCatterall, W. A., Lenaeus, M. J., and Gamal El-Din, T. M. (2020) Structure and pharmacology of 817 \nvoltage-gated sodium and calcium channels. Annu. Rev. Pharmacol. Toxicol. 60, 133–154. 818 \nhttps://doi.org/10.1146/annurev-pharmtox-010818-021757 819 \n 820 \nCombs, S., Kaufmann, K., Field, J. R., Blakely, R. D., and Meiler, J. (2011). Y95 and E444 821 \ninteraction required for high-affinity S-citalopram binding in the human serotonin transporter. ACS 822 \nChem Neurosci., 2(2), 75–81. https://doi.org/10.1021/cn100066p 823 \n 824 \nCraig, R. A., 2nd, Garrison, C. E., Nguyen, P. T., Yarov-Yarovoy, V., & Du Bois, J. (2020). 825 \nVeratridine: A janus-faced modulator of voltage-gated sodium ion channels. ACS Chem Neurosci, 826 \n11(3), 418–426. https://doi.org/10.1021/acschemneuro.9b00621 827 \n 828 \nDavis, I. W., Baker, D. (2009). RosettaLigand docking with full ligand and receptor flexibility. J. 829 \nMol. Biol. 385(2), 381-392. https://doi.org/10.1016/j.jmb.2008.11.010 830 \n 831 \nDeLuca, S., Khar, K., Meiler, J. (2015). Fully flexible docking of medium sized ligand libraries with 832 \nRosettaLigand. PLoS One. 10(7), e0132508. https://doi.org/10.1371/journal.pone.0132508 833 \n 834 \nGees, M., Owsianik, G., Nilius, B., and Voets, T. (2012). TRP channels. Compr. Physiol., 2(1), 563–835 \n608. https://doi.org/10.1002/cphy.c110026 836 \n 837 \nFeske, S., Wulff, H., and Skolnik, E. Y. (2015). Ion channels in innate and adaptive immunity. Annu. 838 \nRev. Immunol., 33, 291–353. https://doi.org/10.1146/annurev-immunol-032414-112212 839 \n 840 \nHagen, N. A., Cantin, L., Constant, J., Haller, T., Blaise, G., Ong-Lam, M., du Souich, P., Korz, W., 841 \nand Lapointe, B. (2017). Tetrodotoxin for moderate to severe cancer-related pain: a multicentre, 842 \nrandomized, double-blind, placebo-controlled, parallel-design trial. Pain Res. Manag. 2017:7212713. 843 \nhttps://doi.org/10.1155/2017/7212713 844 \n 845 \nHalgren, T. A. (1996). Merck molecular force field. I. Basis, form, scope, parameterization, and 846 \nperformance of MMFF94. J. Comput. Chem. 17, 490-519. doi: 10.1002/(SICI)1096-847 \n987X(199604)17:5/6<490::AID-JCC1>3.0.CO;2-P 848 \n 849 \nHalgren, T. A. (1996a). Merck molecular force field. II. MMFF94 van der Waals and electrostatic 850 \nparameters for intermolecular interactions. J. Comput. Chem.17520–552. doi: 10.1002/(SICI)1096-851 \n987X(199604)17:5/6<520::AID-JCC2>3.0.CO;2-W 852 \n 853 \nHalgren, T. A. (1996b). Merck molecular force field. III. Molecular geometries and vibrational 854 \nfrequencies for MMFF94. J. Comput. Chem. 17, 553–586. doi: 10.1002/(SICI)1096-855 \n987X(199604)17:5/6<553::AID-JCC3>3.0.CO;2-T 856 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 30 \n 857 \nHalgren, T. A. (1996c). Merck molecular force field. V. Extension of MMFF94 using experimental 858 \ndata, additional computational data, and empirical rules. J. Comput. Chem. 17:616–641. doi: 859 \n10.1002/(SICI)1096-987X(199604)17:5/6<616::AID-JCC5>3.0.CO;2-X 860 \n 861 \nHalgren, T. A. and Nachbar, R. B. (1996). Merck molecular force field. IV. Conformational energies 862 \nand geometries for MMFF94. J. Comput. Chem. 17, 587–615. doi: 10.1002/(SICI)1096-863 \n987X(199604)17:5/6%3C587::AID-JCC4%3E3.0.CO;2-Q 864 \n 865 \nHalgren, T. A. (1999a). MMFF VI. MMFF94s option for energy minimization studies. J. Comput. 866 \nChem. 20, 720–729. doi: 10.1002/(SICI)1096-987X(199905)20:7<720::AID-JCC7>3.0.CO;2-X 867 \n 868 \nHalgren, T. A. (1999b). MMFF VII. Characterization of MMFF94, MMFF94s, and other widely 869 \navailable force fields for conformational energies and for intermolecular-interaction energies and 870 \ngeometries. J. Comput. Chem. 20, 730–748. doi: 10.1002/(SICI)1096-987X(199905)20:7<730::AID-871 \nJCC8>3.0.CO;2-T 872 \n 873 \nHanwell, M. D., Curtis, D. E., Lonie, D. C., Vandermeersch, T., Zurek, E., Hutchison, G. R. (2012). 874 \nAvogadro: an advanced semantic chemical editor, visualization, and analysis platform. J. 875 \nCheminform. 4(1), 17. https://doi.org/10.1186/1758-2946-4-17. 876 \n 877 \nHawkins, P. C. D., Skillman, A. G., Warren, G. L., Ellingson, B. A., Stahl, M. T. (2010). Conformer 878 \ngeneration with OMEGA: Algorithm and validation using high quality structures from the Protein 879 \nDatabank and Cambridge Structural Database. J. Chem. Inf. Model. 50(4), 572-584. 880 \nhttps://doi.org/10.1021/ci100031x 881 \n 882 \nHarris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., et al., 883 \nand Oliphant, T. E. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. 884 \nhttps://doi.org/10.1038/s41586-020-2649-2 885 \n 886 \nHille, B. (1977). Local anesthetics: hydrophilic and hydrophobic pathways for the drug-receptor 887 \nreaction. J. Gen. Physiol. 69, 497–515. https://doi.org/10.1085/jgp.69.4.497 888 \n 889 \nHille, B. (2001). Ion channels of excitable membranes (3rd ed., pp. 503–537). Sunderland, MA, 890 \nUSA. Sinauer Associates, Inc., Sunderland. 891 \n 892 \nHinman, A., and Du Bois, J. (2003). A stereoselective synthesis of (-)-tetrodotoxin. J. Am. Chem. 893 \nSoc., 125(38), 11510–11511. https://doi.org/10.1021/ja0368305 894 \n 895 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n31 \nHockerman, G. H., Peterson, B. Z., Johnson, B. D., and Catterall, W. A. (1997). Molecular 896 \ndeterminants of drug binding and action on L-type calcium channels. Annual review of 897 \npharmacology and toxicology, 37, 361–396. https://doi.org/10.1146/annurev.pharmtox.37.1.361 898 \n 899 \nJiang, D., Shi, H., Tonggu, L., Gamal El-Din, T. M., Lenaeus, M. J., Zhao, Y., Yoshioka, C., Zheng, 900 \nN., Catterall, W. A. (2020) Structure of the cardiac sodium channel. Cell. 180, 122-134.e10. 901 \nhttps://doi.org/10.1016/j.cell.2019.11.041 902 \n 903 \nKaufmann, K. W., and Meiler, J. (2012). Using RosettaLigand for small molecule docking into 904 \ncomparative models. PloS One, 7(12), e50769. https://doi.org/10.1371/journal.pone.0050769 905 \n 906 \nKim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., Li, Q., Shoemaker, B. A., Thiessen, P. A., 907 \nYu, B., Zaslavsky, L., Zhang, J., Bolton, E. E. (2023). PubChem 2023 update. Nucleic Acids Res. 908 \n51(D1):D1373-D1380. https://doi.org/10.1093/nar/gkac956. 909 \n 910 \nKwong, K., and Carr, M. J. (2015). Voltage-gated sodium channels. Curr. Opin. Pharmacol. 22, 911 \n131–139. https://doi.org/10.1016/j.coph.2015.04.007 912 \n 913 \nLemmon, G., Kaufmann, K., and Meiler, J. (2012). Prediction of HIV-1 protease/inhibitor affinity 914 \nusing RosettaLigand. Chem Biol Drug Des., 79(6), 888–896. https://doi.org/10.1111/j.1747-915 \n0285.2012.01356.x 916 \n 917 \nLiu, H., Hofmann, J., Fish, I., Schaake, B., Eitel, K., Bartuschat, A., Kaindl, J., Rampp, H., Banerjee, 918 \nA., Hübner, H., Clark, M. J., Vincent, S. G., Fisher, J. T., Heinrich, M. R., Hirata, K., Liu, X., 919 \nSunahara, R. K., Shoichet, B. K., Kobilka, B. K., and Gmeiner, P. (2018). Structure-guided 920 \ndevelopment of selective M3 muscarinic acetylcholine receptor antagonists. Proc Natl Acad Sci USA, 921 \n115(47), 12046–12050. https://doi.org/10.1073/pnas.1813988115 922 \n 923 \nMaia EHB, Assis LC, de Oliveira TA, da Silva AM, Taranto AG. Structure-Based Virtual Screening: 924 \nFrom Classical to Artificial Intelligence. Front. Chem. 2020 Apr 28;8:343. 925 \nhttps://doi.org/10.3389/fchem.2020.00343. 926 \n 927 \nMaly, J., Emigh, A.M., DeMarco, K.R., Furutani, K., Sack, J.T., Clancy, C.E., Vorobyov, I., Yarov-928 \nYarovoy, V. (2022) Structural modeling of the hERG potassium channel and associated drug 929 \ninteractions. Front Pharmacol. Sep 16;13:966463. https://doi.org/10.3389/fphar.2022.966463. 930 \n 931 \nMcCormack, K., Santos, S., Chapman, M. L., Krafte, D. S., Marron, B. E., West, C. W., Krambis, M. 932 \nJ., Antonio, B. M., Zellmer, S. G., Printzenhoff, D., Padilla, K. M., Lin, Z., Wagoner, P. K., Swain, 933 \nN. A., Stupple, P. A., de Groot, M., Butt, R. P., & Castle, N. A. (2013). Voltage sensor interaction 934 \nsite for selective small molecule inhibitors of voltage-gated sodium channels. Proc. Natl. Acad. Sci. 935 \nU S A. 110(29), E2724–E2732. https://doi.org/10.1073/pnas.1220844110 936 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 32 \n 937 \nMeiler, J., and Baker, D. (2006). RosettaLigand: protein – small molecule docking with full side-938 \nchain flexibility. Proteins. 65(3), 538–548. https://doi.org/10.1002/prot.21086 939 \n 940 \nMoretti, R., Bender, B. J., Allison, B., and Meiler, J. (2016). Rosetta and the design of ligand binding 941 \nsites. Methods Mol. Biol. (Clifton, N.J.), 1414, 47–62. https://doi.org/10.1007/978-1-4939-3569-7_4 942 \n 943 \nMulligan, V. K., Workman, S., Sun, T., Rettie, S., Li, X., Worrall, et al., Baker, D. (2021). 944 \nComputationally designed peptide macrocycle inhibitors of New Delhi metallo-β-lactamase 1. Proc. 945 \nNatl. Acad. of Sci. USA, 118(12), e2012800118. https://doi.org/10.1073/pnas.2012800118 946 \n 947 \nNanou, E., and Catterall, W. A. (2018). Calcium channels, synaptic plasticity, and neuropsychiatric 948 \ndisease. Neuron 98(3), 466–481. https://doi.org/10.1016/j.neuron.2018.03.017 949 \n 950 \nNerbonne, J. M., and Kass, R. S. (2005). Molecular physiology of cardiac repolarization. Physiol. 951 \nRev. 85(4), 1205–1253. https://doi.org/10.1152/physrev.00002.2005 952 \n 953 \nNguyen, H. M., Singh, V., Pressly, B., Jenkins, D. P., Wulff, H., & Yarov-Yarovoy, V. (2017). 954 \nStructural insights into the atomistic mechanisms of action of small molecule inhibitors targeting the 955 \nKCa3.1 channel pore. Mol. Pharmacol., 91(4), 392–402. https://doi.org/10.1124/mol.116.108068 956 \n 957 \nNguyen, P. T., DeMarco, K. R., Vorobyov, I., Clancy, C. E., & Yarov-Yarovoy, V. (2019). Structural 958 \nbasis for antiarrhythmic drug interactions with the human cardiac sodium channel. Proc. Natl. Acad. 959 \nof Sci. USA, 116(8), 2945–2954. https://doi.org/10.1073/pnas.1817446116 960 \n 961 \nNguyen, P. T., and Yarov-Yarovoy, V. (2022). Towards Structure-Guided Development of Pain 962 \nTherapeutics Targeting Voltage-Gated Sodium Channels. Frontiers in pharmacology, 13:842032. 963 \nhttps://doi.org/10.3389/fphar.2022.842032 964 \n 965 \nNilius B. (2007). TRP channels in disease. Biochim. Biophys. Acta., 1772(8), 805–812. 966 \nhttps://doi.org/10.1016/j.bbadis.2007.02.002 967 \n 968 \nNivón, L.G., Moretti, R., Baker, D. (2013). A Pareto-optimal refinement method for protein design 969 \nscaffolds. PLoS One. 8(4), e59004. https://doi.org/10.1371/journal.pone.0059004 970 \n 971 \nPark, H., Zhou, G., Baek, M., Baker, D., and DiMaio, F. (2021). Force field optimization guided by 972 \nsmall molecule crystal lattice data enables consistent sub-Angstrom protein-ligand docking. J. Chem. 973 \nTheory Comput 17(3), 2000-2010. https://doi.org/10.1021/acs.jctc.0c01184 974 \n 975 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n33 \nPajouhesh, H., Beckley, J. T., Delwig, A., Hajare, H. S., Luu, G., Monteleone, D., Zhou, X., Ligutti, 976 \nJ., Amagasu, S., Moyer, B. D., Yeomans, D. C., Du Bois, J., and Mulcahy, J. V. (2020). Discovery of 977 \na selective, state-independent inhibitor of NaV1.7 by modification of guanidinium toxins. Sci Rep. 978 \n10(1):14791. https://doi.org/10.1038/s41598-020-71135-2 979 \n 980 \nPumroy, R. A., Samanta, A., Liu, Y., Hughes, T. E., Zhao, S., Yudin, Y., Rohacs, T., Han, S., 981 \nMoiseenkova-Bell, V. Y. (2019). Molecular mechanism of TRPV2 channel modulation by 982 \ncannabidiol. Elife. 8e48792. https://doi.org/10.7554/eLife.48792 983 \n 984 \nPumroy, R.A., Protopopova, A.D., Fricke, T.C., Lange, I.U., Haug, F.M., Nguyen, P.T., Gallo, P.N., 985 \nSousa, B.B., Bernardes, G.J.L., Yarov-Yarovoy, V., Leffler, A., Moiseenkova-Bell, V.Y. (2022). 986 \nStructural insights into TRPV2 activation by small molecules. Nat Commun. 13(1):2334. 987 \nhttps://doi.org/.1038/s41467-022-30083-3. 988 \n 989 \nSalomon-Ferrer, R., Case, D. A., Walker, R. C. (2012). An overview of the Amber biomolecular 990 \nsimulation package. WIREs Comput. Mol. Sci. 3, 198-210. https://doi.org/10.1002/wcms.1121 991 \n 992 \nShen, H., Liu, D., Wu, K., Lei, J., and Yan, N. (2019). Structures of human NaV1.7 channel in 993 \ncomplex with auxiliary subunits and animal toxins. Science. 363(6433), 1303–1308. 994 \nhttps://doi.org/10.1126/science.aaw2493 995 \n 996 \nShim, H., Brown, B. M., Singh, L., Singh, V., Fettinger, J. C., Yarov-Yarovoy, V., and Wulff, H. 997 \n(2019). The Trials and Tribulations of Structure Assisted Design of KCa Channel Activators. Front 998 \nPharmacol., 10, 972. https://doi.org/10.3389/fphar.2019.00972 999 \n 1000 \nSmith, S. T., Meiler, J. (2020). Assessing multiple score function in Rosetta for drug discovery. PLoS 1001 \nOne. 15(10), e0240450. https://doi.org/10.1371/journal.pone.0240450 1002 \n 1003 \nStevens, M., Peigneur, S., and Tytgat, J. (2011). Neurotoxins and their binding areas on voltage-1004 \ngated sodium channels. Front Pharmacol., 2(71). https://doi.org/10.3389/fphar.2011.00071 1005 \n 1006 \nSu, M., Yang, Q., Du, Y., Feng, G., Liu, Z., Li, Y., & Wang, R. (2019). Comparative Assessment of 1007 \nScoring Functions: The CASF-2016 Update. J. Chem. Inf. Model., 59(2), 895–913. 1008 \nhttps://doi.org/10.1021/acs.jcim.8b00545 1009 \n 1010 \nTikhonov, D. B., and Zhorov, B. S. (2008). Molecular modeling of benzothiazepine binding in the L-1011 \ntype calcium channel. J. Biol. Chem., 283(25), 17594–17604. 1012 \nhttps://doi.org/10.1074/jbc.M800141200 1013 \n 1014 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n 34 \nThomas-Tran, R., Du Bois, J. (2016). Mutant cycle analysis with modified saxitoxins reveals specific 1015 \ninteractions critical to attaining high-affinity inhibition of hNaV1.7. Proc. Natl. Acad. Sci. U S A. 1016 \n113(21), 5856-61. https://doi.org/10.1073/pnas.1603486113 1017 \n 1018 \nVallat, R. (2018). Pingouin: statistics in Python. Journal of Open Source Software, 3(31), 1026. 1019 \nhttps://doi.org/10.21105/joss.01026 1020 \n 1021 \nVirtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al., and 1022 \nvan Mulregt, P. (2020). SciPy 1.0: fundamental algorithms for scientific computing in 1023 \nPython. Nature methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2 1024 \n 1025 \nVu, S., Singh, V., Wulff, H., Yarov-Yarovoy, V., Zheng, J. (2020). New capsaicin analogs as 1026 \nmolecular rulers to define the permissive conformation of the mouse TRPV1 ligand-binding pocket. 1027 \nElife. Nov 9;9:e62039. https://doi.org/10.7554/eLife.62039. 1028 \n 1029 \nWang, S. Y., Tikhonov, D. B., Mitchell, J., Zhorov, B. S., and Wang, G. K. (2007). Irreversible block 1030 \nof cardiac mutant Na+ channels by batrachotoxin. Channels (Austin), 1(3), 179–188. 1031 \nhttps://doi.org/10.4161/chan.4437 1032 \n 1033 \nWulff, H., Christophersen, P., Colussi, P., Chandy, K. G., and Yarov-Yarovoy, V. (2019). Antibodies 1034 \nand venom peptides: new modalities for ion channels. Nat. Rev. Drug Discov. 18(5), 339-357. 1035 \nhttps://doi.org/10.1038/s41573-019-0013-8 1036 \n 1037 \nYang, T., Smith, J. A., Leake, B. F., Sanders, C. R., Meiler, J., and Roden, D. M. (2013). An 1038 \nallosteric mechanism for drug block of the human cardiac potassium channel KCNQ1. Mol 1039 \nPharmacol., 83(2), 481–489. https://doi.org/10.1124/mol.112.081513 1040 \n 1041 \nYang, F., Xiao, X., Cheng, W., Yang, W., Yu, P., Song, Z., Yarov-Yarovoy, V., & Zheng, J. (2015). 1042 \nStructural mechanism underlying capsaicin binding and activation of the TRPV1 ion channel. Nat. 1043 \nChem. Biol., 11(7), 518–524. https://doi.org/10.1038/nchembio.1835 1044 \n 1045 \nZhang, Y., Wang, K., and Yu, Z. (2020). Drug Development in Channelopathies: Allosteric 1046 \nModulation of Ligand-Gated and Voltage-Gated Ion Channels. J. Med. Chem., 63(24), 15258–15278. 1047 \nhttps://doi.org/10.1021/acs.jmedchem.0c01304 1048 \n 1049 \nZheng, J., and Trudeau, M.C. (2015). Chapter 28. TRPV channels. In Handbook of Ion Channels (1st 1050 \ned., pp. 427-431). Boca Raton, FL, USA. CRC Press. https://doi.org/10.1201/b18027 1051 \n 1052 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint \n\n  Running Title \n \n35 \nZhao, Y., Huang, G., Wu, J., Wu, Q., Gao, S., Yan, Z., Lei, J., and Yan, N. (2019a). Molecular basis 1053 \nfor ligand modulation of a mammalian voltage-gated Ca2+ channel. Cell. 177(6), 1495-1506.e12. 1054 \nhttps://doi.org/10.1016/j.cell.2019.04.043 1055 \n 1056 \nZhao, Y., Huang, G., Wu, Q., Wu, K., Li, R., Lei, J., Pan, X., Yan, N. (2019b). Cryo-EM structures 1057 \nof apo and antagonist-bound human Cav3.1. Nature. 576(7787), 492–497. 1058 \nhttps://doi.org/10.1038/s41586-019-1801-3 1059 \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}