1. Introduction

High throughput virtual screening (HTVS) is a widely used technique for determining the binding affinity of a large number of ligands to a target receptor ^1,2,3. The output of HTVS provides a qualitative comparison of a ligand library, which is to say the ligands can be ranked from best to worst based on receptor binding affinity. The speed and accuracy of HTVS is dependent on how well the binding pocket of the receptor has been characterized. Commonly, the binding pocket for HTVS docking experiments is determined from the crystal structure of a biomacromolecular receptor with bound ligand.

Jacob et al. (2011) provided a step-by-step description of the process of initiating a docking experiment beginning with a pdb file for a ligand-receptor complex ⁴. Briefly, the authors obtained the protein database pdb file for Aplysia californica-acetylcholine binding protein (Ac-AChBP) (accession number 2BR8) containing a-conotoxin (a-CTx) PnIA[A10L:D14K] bound to the ligand binding domain. The Cartesian coordinate information specifying each atom of the ligand was removed from the pdb file of the ligand-receptor complex leaving the remaining pdb file consisting of only apoenzyme information. AutoDock Tools was then used to generate a grid parameter file encompassing this binding pocket, allowing for HTVS docking experiments to determine how new a-CTx ligands bind to the Ac-AChBP. HTVS experiments have been used to discover DNA gyrase inhibitors, improve knowledge of dihydrodipicolinate reductase binding paradigms, discover STAT3 inhibitors, optimizing JAK2 inhibitors, kinase inhibitors, G-quadruplex ligands to demonstrate inert group 9 metals inhibiting protein-protein interfaces, and others ^{5,6,7,8,9,10,11,12}.

In this study, a-CTx peptide ligands were docked to a homology model of the a₃β₂-nicotinic acetylcholine receptor (nAChR) using DockoMatic version 2.0 software suite ¹³. The HTVS was conducted using a-CTxMII as a ligand because experimental data supports a very high affinity of the peptide for the a₃β₂-nAChR isoform (IC₅₀ = 0.5 nM) ¹⁴. The sequence for a-CTxMII is GCCSNPVCHLEHSNLC with disulfide bonds connecting Cys2-Cys8 and Cys3-Cys16. The peptide has a well-established a-helix originated at Pro6 and ending at His12. Coordinate files have been deposited in the Research Collaboratory for Structural Bioinformatics (RCSB) ¹⁵ representing a-CTxMII from two separate groups; in this investigation we used the original work submitted by Shon et al. ¹⁶ rather than Hill et al. ¹⁷. The structure generated by Shon et al. was based on traditional aqueous solution conditions for NMR experiments in contrast to the further refinement in structure gained by Hill et al. as a result of solvent condition manipulation to reduce polarity by addition of either acetonitrile or trifluoroethanol. In their work, Hill et al. acknowledge the global fold of the peptide remains the same as that reported by Shon et al. For our purposes, either system would have sufficed as a starting point for computational studies. The structure of a-CTxMII contains electropositive amino acids at either end of the helix (Asn5 and His12) that have been proposed to be significant for ligand to receptor interaction ^18,19. Two additional analogs of a-CTxMII were also evaluated, a-CTxMII[E11A] and a-CTxMII[N5R:E11A:H12K].

Conotoxins are small, cysteine-rich peptides found in the venom of snails belonging to the genus Conus. a-CTxs selectively bind and inhibit nAChRs leading to a great deal of research into their use as analgesics and treatments for Parkinson’s disease ^20,21,22. α-CTxMII from the venom of Conus magus, is a well-studied molecular probe for α₃β₂-nAChRs. The single point mutant a-CTxMII[E11A] was determined to shift binding selectivity from α₃β₂-nAChR isoforms to preferential binding of α₆β₂β₃-containing nAChRs ²³. The triple mutant, a-CTxMII[N5R:E11A:H12K], was designed to have positively charged amino acids at both ends of the α-helix, which was hypothesized to enhance receptor binding affinity while still retaining the E11A selectivity for α₆β₂β₃-containing nAChRs ¹⁹. Figure 1 shows the structures of the three peptides, where the top structures indicate the larger size of the triple mutant (right), which measures 16 Å from end-to-end of the α-helix, in comparison to either a-CTxMII (left) or a-CTxMII[E11A] (middle), which are both on the order of 11 Å from end-to-end of the α-helix. The bottom set of structures shows the electrostatic topography of each peptide for visualization of the removal of negative charge (red) in position 11 for the a-CTxMII[E11A] mutant, and the enhancement of positive charge (blue) character for the triple mutant.

Figure 1.Top, from left to right: distance between amino acid positions 5 and 12 of a-CTxMII (~11 Å), a-CTxMII[E11A] (~11 Å), and a-CTxMII[N5R:E11A:H12K] (~16 Å). Bottom, from left to right: electrostatic topography of a-CTxMII, a-CTxMII[E11A], and a-CTxMII[N5R:E11A:H12K] showing the decrease in negative charge (red) and increase in positive charge (blue) with each successive amino acid mutation

DownLoad: Full-Size ImgPowerPoint

The aim of this study was to identify the reproducibility and limitations of HTVS for each ligand individually and for the three ligands compared to one another for binding to a homology model of the a₃β₂-nAChR isoform.

2. Materials and Methods

DockoMatic is a comprehensive open source computational application with a graphical user interface (GUI) designed to facilitate HTVS; the docking engine used by DockoMatic is AutoDock v 4.2 ²⁴. DockoMatic was installed and made accessible on a high performance computing (HPC) cluster consisting of 32 AMD processors per node with 16 nodes available.

DockoMatic was then used to generate mutant analogs of a-CTxMII beginning with the pdb file (accession number 1M2C) originating from nuclear magnetic resonance spectroscopy structure elucidation ¹⁶. The two mutations to a-CTxMII, a-CTxMII[E11A] and a-CTxMII[N5R:E11A:H12K], were generated using the Treepack program, subordinate to the DockoMatic GUI ⁴. These three peptides were docked to the α₃β₂-nAChR for a total of ten trials each. Every trial consisted of 100 AutoDock cycles, resulting in 1000 docking trials per peptide ²⁵.

3. Results and Discussion

The results from HTVS provide a qualitative ranking of ligands using a binding energy in kcal/mol as the unit of measurement. Analysis of the lowest free energy for binding (DG) resultant from the ten trials for each of the three peptides showed significant variation (Table 1).

Table 1.Lowest free energy of binding value for a-CTxMII, a-CTxMII[E11A], and a-CTxMII[N5R:E11A:H12K] docking to a homology model for the a₃β₂-nAChR over 10 trials

Conotoxin	trial1	trial2	trial3	trial4	trial5	trial6	trial7	trial8	trial9	trial10
MII	−14.96	−16.23	−15.32	−13.78	−16.27	−15.39	−14.53	−14.14	−15.32	−14.69
E11A	−15.98	−17.58	−16.22	−16.25	−16.26	−15.92	−17.58	−16.61	−7.14	−16.55
Triple	−17.3	−3.9	−8.21	−9.02	−16.95	−17.4	−17.34	−17.74	−7.43	−17.92

 | Show Table

DownLoad: CSV

The average free energy of binding for the most favorable bound ligand, resultant from 100 docking cycles over all ten trials, for a-CTxMII was −15.063 ± 0.816 kcal/mol, a-CTxMII[E11A] was −15.069 ± 3.034 kcal/mol, and a-CTxMII[N5R:E11A:H12K] was −13.321 ± 5.483 kcal/mol (Table 2). From these data, one could conclude that more favorable interaction with the receptor is achieved for a-CTxMII and a-CTxMII[E11A] as compared to a-CTxMII[N5R:E11A:H12K], but the large standard deviation for the triple mutant indicated docking reproducibility was a problem and more stringent analysis of the results was required. We speculate that the triple mutant is subject to greater variability in docking result for two reasons: (1) the peptide is 5 Å larger than either a-CTxMII or a-CTxMII[E11A] resulting in greater difficulty to bind to the same receptor binding site, and (2) the triple mutant has greater charge character at either end of the peptide, lending itself to a larger number of binding opportunities that may constitute local energy wells rather than global energy minimum binding interactions.

Table 2.Correlation between average energy for docked ligand-receptor complex and binding pose similarity across ten docking trials consisting of 100 cycles per trial

Peptide	Superimposed Docked Peptide	Average ΔG (kcal/mol)	Standard Deviation	Root Mean Square Deviation (Å)
MII	1, 6, 8, 9, 10	−15.063	±0.816	0.510
E11A	1, 3-8, 10	−15.609	±3.034	0.526
Triple	1, 4, 5, 7, 8, 10	−13.321	±5.483	3.459

 | Show Table

DownLoad: CSV

From the data provided in Table 2, it can be seen that the root mean square deviation (RMSD) between similar binding pose peptides for a-CTxMII and a-CTxMII[E11A] were roughly comparable and within specifications that are considered generally acceptable (e.g. ≤ 0.5 Å). These data are in marked contrast to the RMSD for the triple mutant, which is much higher (3.459 Å), indicating greater variability in ligand binding pose between docking trials. It should be noted that trial 2 for the triple mutant was omitted from calculation of average free energy of binding because it was well outside of the standard deviation for calculated free energy, and the bound structure did not converge on a common binding site as was observed for the other ligand structures. As stated previously, the larger size and charge character of the triple mutant is suspected to affect convergence across AutoDock trials. Although AutoDock was the docking engine used in this investigation, we suspect other docking engines will be subject to the same limitations. The reason for this is because HTVS protocols limit the number of rotatable bonds (e.g. AutoDock default value = 32) allowed per docking trial to achieve efficient computational processing. To remove the rotatable bond limitation, requires more computationally intensive investigation which then restricts the number of screening experiments that can be run.

Protein data base files for the thirty most energetically favorable ligand binding poses (top results over ten trials for the three peptides) were visualized to identify whether inconsistency in binding energy correlated to significant differences in bound ligand positioning. Figure 2 shows the ligand-receptor complexes where the image to the left represents all thirty peptide binding poses for the three ligands and the image to the right provides a view of reproducible docking pose that is conserved for nineteen of the docked structures derived from HTVS.

Figure 2.HTVS results for all thirty docking trials (left); nineteen of the thirty trials identified a single conserved binding pose (right)

DownLoad: Full-Size ImgPowerPoint

Close visual inspection of the nineteen ligands with conserved binding pose to the receptor was performed to identify that five of the peptides are a-CTxMII, eight are a-CTxMII[E11A], and the remaining six ligands are a-CTxMII[N5R:E11A:H12K] (Figure 3). In order to limit computational time, AutoDock generates these binding poses in an efficient process by limiting the number of rotatable bonds to 32 for the ligand, and not allowing dynamic motion for the receptor ²⁶. This computational docking method is limited in its ability to optimize side chain interactions between the ligand and the receptor for large systems like those explored in this exercise. As a result, the calculated binding energy for peptide ligands is largely variable because of side chain immobility inherent to docking engines, such as AutoDock.

The conclusion from this investigation is that HTVS should be performed over multiple trials, followed by visualization of the docking results to identify commonalities in binding pose. The qualitative ranking of ligands in this exercise did not extend between ligands. Thus, it was only through visual inspection of ligand binding pose and validation of binding pose across a minimum of ten trials, that the reproducibility of results could be assessed. To achieve ten trials for each ligand-receptor docking, the authors used DockoMatic, drawing on the docking engine AutoDock. Alternatively, AutoDock Vina, which also takes advantage of parallel processing across computing clusters, could have been used for efficient docking ²⁷. Regardless of the docking engine, it is our finding that manual inspection of results and visual comparison of binding pose across trials is required for confident binding pose identification. When data are available in the RCSB, visual inspection of docking results should be validated by comparison to ligand-receptor structure obtained experimentally by x-ray crystallography. This study suggests that results obtained through HTVS should segue into energy calculation using molecular dynamics simulations where ligand and receptor dynamics are permitted. HTVS can be a valuable tool to achieve a starting point for rigorous and time-intensive molecular dynamics simulation experiments or other more thorough methodologies ²⁸.

Figure 3.Five of the nineteen peptides with consistent binding pose represent a-CTxMII (left), eight of them are a-CTxMII[E11A] (middle), and the remaining six ligands are a-CTxMII[N5R:E11A:H12K] (right)

DownLoad: Full-Size ImgPowerPoint

Acknowledgements

Molecular graphics and analyses were performed with the UCSF Chimera package. Chimera was developed by the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco (supported by NIGMS P41-GM103311). Visualization and structure analysis was also performed with the PyMOL Molecular Graphics System, Version 1.8 Schrödinger, LLC. This material is based in part upon work supported by the National Science Foundation under Grant No. 1229709

Conflict of interest

The authors declare no conflicts of interest in the reporting of this study.