Ligand generation

Ligands were extracted as Structure Data Files (.sdf) from PubChem (Kim et al, 2023). Using the Avogadro software (Hanwell et al, 2012), each ligand structure underwent bond-correction, protonation at pH 7.4, and energy minimization using the Merck Molecular Force Field (Halgren, 1996a; b; c; Halgren and Nachbar, 1999; Halgren, 1999a; b). The resulting models were saved as Tripos Mol2 (.mol2) files. The protonation and bond order of saxitoxin and tetrodotoxin was matched to experimentally reported work (Hinman and DuBois, 2003; Thomas-Tran and Du Bois, 2016). Both experimentally resolved structures of verapamil docked to rabbit Ca_V1.1 were tested (Zhao et al, 2019).

Next, using Amber Tools’ Antechamber protocol, the partial atomic charge, atom, and bond type assignments for each ligand were AM1-BCC corrected (Case et al, 2021; Salomon-Ferrer, Case, and Walker, 2012; Appendix S1). AM1-BCC correction is commonly used in Rosetta-based ligand docking protocols (Smith and Meiler, 2020; Park et al, 2021) and has demonstrated a similar performance correlation with other RosettaLigand input preparation protocols (Smith and Meiler, 2020).

The AM1-BCC corrected ligands were used for subsequent steps specific to each method. For RosettaLigand, an in-house script (Appendix S2) using the OpenEye Omega toolkit (Hawkins et al, 2010) was used to generate the conformer library, followed by Rosetta to generate the associated ligand parameters file. For GALigandDock, the input conformer was generated using the RosettaGenFF crystal structure prediction protocol (Park et al, 2021; Appendix S3), taking the lowest energy packing arrangement as input.

Ion channel preparation

Ion channel structures were downloaded from the Protein Data Bank (Berman et al, 2000). Prior to RosettaLigand docking, structures were relaxed with backbone constraints using the RosettaRelax protocol (Nivón, Moretti, and Baker, 2013). This protocol allows the repacking of protein sidechains and minimization of the structure into the Rosetta score function for comparison between poses. The lowest energy pose from 100 relaxed poses was used for docking.

RosettaLigand docking

RosettaLigand docking was performed using previously described RosettaScripts protocols (Davis and Baker, 2009; DeLuca, Khar, and Meiler, 2015; Appendix S4-S7). Briefly, the initial placement of all ligands into their respective ion channels was performed by superimposing the initial ligand poses onto the experimental ligand coordinates.

RosettaLigand uses grid-based sampling to score the ion channel – ligand interface, with a low-resolution and high-resolution sampling phase. For an unbiased sampling of the local space, an initial transformation during the low-resolution sampling phase is performed on the initial ligand position using Rosetta’s Transform mover, The low-resolution Transform mover is a grid-based Monte Carlo simulation, where the ligand is translated up to 0.2 Å and rotated up to 20 degrees per iteration, for a total of 500 iterations. A box size of 7-8 Å restrains the ligand center, while all ligand atoms are constrained to a grid with a user defined width to prevent ligand scoring outside the target site. This scoring grid width was calculated uniquely for each ligand to ensure all ligand conformers would not automatically fail the Transform step. As used previously, the scoring grid width was calculated as the maximum conformer atom-atom distance plus twice the box size value used in the Transform mover (Moretti, Bender, Allison, and Meiler, 2016; Appendix S4). The pose with the lowest score is then used as the starting pose for high-resolution docking. For high-resolution docking using Rosetta’s HighResDocker mover, six docking cycles of rotamer sampling are performed with repacking of sidechains every third iteration. Lastly, the ion channel – ligand complex is minimized using Rosetta’s FinalMinimizer mover, and the interface scores are reported using the InterfaceScoreCalculator.

For all RosettaLigand docking runs, the ligand_soft_rep and hard_rep scoring functions were reweighted based on previous work assessing Rosetta score functions with the Comparative Assessment of Score Function 2016 (CASF-2016) dataset (Su et al., 2019; Smith and Meiler, 2020). Specifically, reweights of Coulombic electrostatic potential (fa_elec), Lennard-Jones repulsive energy between atoms in the same residue (fa_intra_rep), sidechain-backbone hydrogen bond energy (hbond_bb_sc), sidechain-sidechain hydrogen bond energy (hbond_sc), and Ramachandran preferences (rama) were applied for the soft-repulsive and hard-repulsive docking phases (Appendix S6).

For each docking run, either 20,000 poses or 100,000 poses were generated to assess if there are statistically significant differences in the lowest recorded pose RMSD (RMSD_Min) to experimental ligand coordinates. In RosettaLigand, a ligand interface is defined either by a representative ligand atom (a “neighbor atom”, defined as the geometric center of mass by default), or all ligand atoms, relative to all ion channel C_β atoms within a specified radius from the ligand (commonly default to 6 or 7 Å). For RosettaLigand, two mutually exclusive ligand area interface modes are available for scoring the pose interface: the ligand neighbor atom cutoffs mode (add_nbr_radius=”True”), and the all ligand atom cutoffs mode (all_atom_mode=”True”). In previous work, both modes were used (Moretti, Bender, Allison, and Meiler, 2016; Smith and Meiler, 2020), thus, both modes were evaluated for any differences in performance. Four individual RosettaLigand docking sets were performed for each PDB structure, by combining different pose totals and ligand area interface modes, for performance comparison.

GALigandDock docking

As described in the original study, docking was performed using the RosettaScripts’ GALigandDock mover in DockFlex mode (Park et al, 2021; Appendix S8-S10). Replicate runs of GALigandDock were performed in parallel for each structure evaluated. Each run consisted of 20 generations with a pool of 100 poses, where each generation updates the pool by total energy. By default, GALigandDock outputs the top 20 structures from the final generation, however, for this study, the entire pool of poses were used. For each docking run, 20,000 poses (1,000 runs) or 100,000 poses (5,000 runs) were output with a padding value of 2 Å, 4 Å, or 7 Å to test for statistical differences in the RMSD_Min. In total, six individual GALigandDock docking sets were performed for each PDB structure, by combining different pose totals and padding sizes, for performance comparison.

P_Near calculation

P_Near is a quantitative metric to evaluate an energy funnel’s quality from a population of poses, with the lowest-scoring pose as the converged minima, or reference state. P_Near was calculated as previously defined (Bhardwaj, Mulligan, Bahl, et al., 2016), and applied as previously described for ligand docking (Smith and Meiler, 2020), with the reference state being the experimental ligand coordinates. P_Near ranges from 0 (protein will not converge to reference state) to 1 (protein will always converge to the reference state).
The numerator is the Boltzmann probability of an individual pose being near the reference state, governed by the “likeness” parameter (λ), and the thermal energy, the product of the Boltzmann constant and absolute temperature (k_BT). The denominator is the partition function of a canonical ensemble. For small molecule docking, the root mean square deviation (RMSD) of the pose’s ligand coordinates relative to experimental ligand coordinates when the protein pose is superimposed to the entire protein structure, or the receptor site is used. The energy scoring metric (E_i) is the interface energy: the energy solely composed of protein-ligand interactions.

A previous Rosetta small-molecule docking study (Park et al., 2021) categorized reference-like poses at 1 Å or 2 Å RMSD without calculating P_Near, while another study evaluating RosettaLigand performance calculated P_Near using λ =1.5 Å, and k_BT =0.62 (Meiler and Smith, 2020). Therefore, we calculated P_Near with reference-like poses defined by λ=2.0 Å and k_BT =0.62 Rosetta energy units (REU), however we calculated P_Near using all previously reported values for the parameters λ and k_BT (Tables S4.1-4.10). A specific P_Near cutoff indicative for drug discovery pipelines has not been established, hence we will refer to P_Near ≥ 0.5 as a ‘first pass’ cutoff for this study when evaluating energy funnel convergence to the experimental ligand coordinates.

Statistical tests

Tests for normality, heteroscedasticity, and Pearson’s correlation between covariates were performed in Python using NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), and Pingouin (Vallat, 2018) with a significance level (α) of 0.05. Population data consisted of the RMSD_Min for each docking set. Not all RMSD_Min data for each docking set fit a normal distribution. Since RMSD_Min data should skew towards zero Å, a lognormal base 10 transformation was applied to all RMSD_Min population data when comparing across docking sets. Shapiro Wilk tests for normality and Levene heteroscedasticity tests were performed to ensure transformed data were normal and of equal variance, respectively.

Covariates for both methods were the number of rotatable bonds of the ligand, the number of ligand heavy atoms, and the resolution of the structure. Statistical tests for RosettaLigand also included the number of conformers generated, transformed with a lognormal base 10, as a covariate. The ligand molecular weight was initially included as a covariate but was discarded after identifying strong Pearson’s correlation with the number of ligand heavy atoms (Table S3, Figure S1.).

To assess RMSD_Min across all docking sets for a given method, a repeated measures factorial analysis of variance (ANOVA) with covariates were performed using IBM SPSS (version 29). For RosettaLigand, the two factors were sample size and ligand area interface mode, while for GALigandDock the two factors were sample size and padding value. Mauchly’s test for sphericity was performed for levels of padding (p = 0.63).