(such as mass, overall charge, etc.). These properties included AlogP (a log of the calculated octanol-water partition coefficient, which measures the extent of a substance hydrophilicity or hydrophobic- ity), molecular weight, formal charge, the number of hydrogen bond donors and acceptors, and the number of rotatable bonds. All molecules were prepared as previously described, and a conforma- tional set of 50 "best-quality" low-energy conformations was generated for each molecule. All conformers within 20 kcal/mol from the global energy minimum were included in the set. The dataset was screened using the "ligand pharmacophore mapping" protocol (DS2.5), with the minimum interference distance set to 1 A˚ and the maximum omitted features set to 0. All other protocol parameters were maintained at the default settings. To analyze enrichment results and select the best pharmaco- phore model for subsequent virtual screening, ROC curves were constructed for each model, where the fraction of identified known binders (true positives, representing sensitivity) was plotted against the fraction of identified library molecules (false positives; 1- specificity). Based on this analysis, the best pharmacophore model was selected for virtual screening purposes. Generation of the DrugBank data set and virtual screening The DrugBank database [40] (release 2.0), which contains ,4900 drug entries, including 1382 FDA-approved small- molecule drugs, 123 FDA-approved biotech (protein/peptide) drugs, 71 nutraceuticals, and over 3240 experimental drugs, was used for Virtual Screening. The database was filtered, based on the average molecular properties of known hPKR antagonists 6 4SD (standard deviation). These properties included AlogP, molecular weight, the number of hydrogen bond donors and acceptors, the formal charge, and the number of rotatable bonds. The liberal 64SD interval was chosen because the calculated range of molecular properties of the known antagonists was very narrow. Molecules were retained only if their formal charge was neutral or positive, since the known compounds were positively charged. This resulted in a test set containing 432 molecules. All molecules were prepared as previously described, and a set of 50 "best-quality" low-energy conformations was generated for each molecule; all conformations were within 20 kcal/mol from the global energy minimum. The data set was screened against the pharmacophore model (chosen from the ROC analysis) using the "ligand pharmacophore mapping" protocol in DS2.5. All protocol settings were maintained at default settings except for minimum interference distance, which was set to 1 A˚ and the maximum omitted features was set to 0. To prioritize the virtual hits, fit values were extracted, to reflect the quality of molecule mapping onto the pharmacophore. Only molecules with fit values above the enrichment ROC curve cutoff that identifies 100% of the known PKR antagonists (FitVa- lue>2.85746) were retained as virtual hits for further analysis. The similarity between the virtual hits and known small- molecule PKR antagonists was evaluated by calculating the Tanimoto coefficient distance measure using the ’Find similar molecules by fingerprints’ module in DS2.5, which calculates the number of AND bits normalized by the number of OR bits, according to SA/(SA+SB+SC), where SA is the number of AND bits (bits present in both the target and the reference), SB is the number of bits in the target but not the reference, and SC is the number of bits in the reference but not the target. Small-Molecule Docking Molecular docking of the small-molecule hPKR antagonists dataset (active and inactive molecules), as well as of virtual hits, to Small-Molecule Binders of Prokineticin Receptors the hPKR1 homology model, was performed using LigandFit [34] as implemented in DS2.5. LigandFit is a shape complementary- based algorithm that performs flexible ligand-rigid protein docking. In our experiments, the binding site was defined as a 284.8 A˚3 TM cavity area, surrounded by binding site residues identified using the energy-based methods described above. Default algorithm settings were used for docking. The final ligand poses were selected based on their empirical LigScore docking score [41]. Here we used the (default) Dreiding force field to calculate the VdW interactions. All docking experiments were conducted on a model without extracellular and intracellular loops. Loop configurations are highly variable among the GPCR crystal structures [42]. Therefore, deleting the loops in order to reduce the uncertainty stemming from inaccurately predicted loops is a common practice in the field [43,44,45]. To further validate our protocol, we also performed molecular redocking of the small-molecule partial inverse agonist carazolol and the antagonist cyanopindolol to their original X-ray structures from which loops were deleted, and to loopless homology models of b1adr and b2adr using LigandFit, as previously described. As in the case of docking to the hPKR1 model, this procedure was performed on loopless X-ray structures and models. The binding site was identified from receptor cavities using the "eraser" and "flood-filling" algorithms, as implemented in DS2.5. The highest scoring LigScore poses were selected as the representative solutions. The ligand-receptor poses were compared to the corresponding X-ray complexes by (1) calculating the root mean square deviation (RMSD) of heavy ligand atoms from their respective counterparts in the crystallized ligand after superposi- tion of the docked ligand-receptor complex onto the X-ray structure; (2) calculating the number of correct atomic contacts in the docked ligand-receptor complex compared with the X-ray complex, where an atomic contact is defined as a pair of heavy ligand and protein atoms located at a distance of less than 4 A˚; and by (3) comparing the overall number of correctly predicted interacting residues in the docked complex to the X-ray complex (where interacting residues are also defined as residues located less than 4 A˚ from the ligand). Small-molecule docking analysis The resulting ligand poses of the known hPKR antagonists were analyzed to identify all ligand-receptor hydrogen bonds, charged interactions, and hydrophobic interactions. The specific interactions formed between the ligand and binding site residues were quantified to determine the best scoring pose of each ligand (active and inactive). For each ligand pose, a vector indicating whether this pose forms a specific hydrogen bond and/or hydrophobic p interaction with each of the binding site residues was generated. The data were hierarchically clustered using the clustergram function of the bioinformatics toolbox in Matlab version 7.10.0.499 (R2010a). The pairwise distance between these vectors was computed using the Hamming distance method, which calculates the percentage of coordinates that differ. For a m-by-n data matrix X, which is treated as m (1-by-n) row vectors x1, x2, …, xm, the distance between the vector xs and xt is defined as follows: dst~ # xsj=xtj =n where # is the number of vectors that differ. The poses of the virtual hits ligands were further filtered using structure-based constraints derived from analyzing the interactions between known PKR antagonists and the receptor, obtained in the PLoS ONE | www.plosone.org 4 November 2011 | Volume 6 | Issue 11 | e27990