known binders docking section of this work. The constraints included (1) an electrostatic interaction between the ligand and Glu1192.61, (2) at least one hydrogen bond between the ligand and Arg1443.32, and/or Arg3076.58, and (3) at least two hydrophobic interactions (p-p or p-cation) between the ligand and Arg1443.32 and/or Arg3076.58. Evolutionary selection analysis Evolutionary selection analysis of the PKR subtypes’ coding DNA sequences was carried out using the Selecton server (version 2.4) [46,47]. The Selecton server is an on-line resource which automatically calculates the ratio (v) between non-synonymous (Ka) and synonymous (Ks) substitutions, to identify the selection forces acting at each site of the protein. Sites with v.1 are indicative of positive Darwinian selection, and sites with v,1 suggest purifying selection. As input, we used the homologous coding DNA sequences of 13 mammalian species for each subtype, namely, human, rat, mouse, bovine, rabbit, panda, chimpanzee, orangutan, dog, gorilla, guinea pig, macaque and marmoset. We used the default algorithm options and the obtained results were tested for statistical significance using the likelihood ratio test, as implemented in the server. Results SAR analysis highlights molecular features essential for small-molecule antagonistic activity A review of the literature revealed a group of non-peptidic compounds that act as small-molecule hPKR antagonists, with no apparent selectivity toward one of the subtypes [17,18,19,20,38]. The reported compounds have either a guanidine triazinedione or a morpholine carboxamide scaffold. We decided to perform structure-activity relationship (SAR) analysis of the triazine-based compounds, owing to the more detailed pharmacological data available for these compounds [17,18,19,20]. SAR analysis of the reported molecules with and without antagonistic activity toward hPKR provides hints about the geometrical arrangement of chemical features essential for the biological activity. By comparing pairs of active and inactive compounds that differ in only one functional group, one can determine the activity-inducing chemical groups at each position. To this end, we constructed a dataset of 107 molecules identified by high-throughput screening. This included 51 molecules that we defined as inactive (Ca2+ mobilization IC50 higher than 1 mM), and 56 molecules defined as active (IC50 below 0.05 mM). All compounds share the guanidine triazinedione scaffold (see figure 2), which includes (a) a heterocyclic ring baring three nitrogen atoms and two oxygen atoms, and (b) a guanidine group, which is attached to the main ring by a linker (position Q in figure 2). Where possible, the dataset was divided into pairs of active and inactive molecules that differ in only one functional group. This resulted in 13 representative pairs of molecules that were used to determine which specific chemical features in these molecules are important for antagonistic activity, in addition to the main triazine ring and guanidine group. As shown in figure 2, the four variable positions in the scaffold - A1, D, L2, and Q, were compared among the 13 pairs, and the activity-facilitating chemical groups at each position were determined. These include the following features: (1) Positions A1 and D require an aromatic ring with a hydrogen bond acceptor in position 4 of the ring. (2) Position L2 may only accept the structure -NH(CH2)-. Small-Molecule Binders of Prokineticin Receptors (3) Position Q may include up to four hydrogen bond donors, a positive ionizable feature, and an aromatic ring bearing a hydrogen bond acceptor. In conclusion, the SAR analysis revealed 2D chemical features in the molecules, which may be important for receptor binding and activation. Next, these features will be used to generate ligand- based pharmacophore models for virtual screening (next section) and in docking experiments to determine the plausible ligand- receptor contacts (see below). Ligand-based virtual screening for novel PKR binders To identify novel potential hPKR binders, we utilized a ligand- based procedure in which molecules are evaluated by their similarity to a characteristic 3D fingerprint of known ligands, the pharmacophore model. This model is a 3D ensemble of the essential chemical features necessary to exert optimal interactions with a specific biological target and to trigger its biological response. The purpose of the pharmacophore modeling procedure is to extract these chemical features from a set of known ligands with the highest biological activity. The two most potent (IC50,0.02 mM for intracellular Ca2+ mobilization) hPKR antagonists were selected from the dataset described in the previous section, to form the training set (compounds 1 and 2, figure S3). In addition, we also incorporated data from a third compound published recently (compound 3 in figure S3), to ensure good coverage of the available chemical space [38]. The HipHop algorithm [37] was used to generate common features of pharmacophore models. This algorithm generated 10 different models, which were first tested for their ability to identify all known active hPKR triazine-based antagonists (data not shown). During the pharmacophore generation and analysis procedure, we also projected the knowledge generated during our 2D SAR analysis onto the 3D pharmacophore models, and chose those that best fit the activity-facilitating chemical features identified in the 2D SAR analysis previously described. The two best models, which recaptured the highest number of known active hPKR binders and included all required 2D features deduced from the SAR analysis, were chosen for further analysis. The 3D spatial relationship and geometric parameters of the models are presented in figure 3A. Both models share a positive ionizable feature and a hydrogen bond acceptor, corresponding to the N3 atom and O1 atoms on the main ring, respectively (figure 2). However, the models vary in the degree of hydrophobicity tolerated: model 2 is more restrictive, presenting one aromatic ring feature and one hydrophobic feature, whereas model 1 is more promiscuous, presenting two general hydrophobic features. The aromatic/hydrophobic features correspond to positions A1 and D of the scaffold (figure 2). Figure 3A also shows the mapping of one of the training set molecules onto the pharmacophore model. All four features of both models are mapped well, giving a fitness value (FitValue) of 3.602 and 3.378 for hypotheses 1 and 2, respectively. The fitness value measures how well the ligand fits the pharmacophore. For a four-feature pharmacophore the maximal FitValue is 4. Next, we performed an enrichment study to ultimately evaluate the pharmacophore model’s performance. Our aim was to verify that the pharmacophores are not only able to identify the known antagonists, but do so specifically with minimal false positives. To this end, a dataset of 56 known active hPKR small-molecule antagonists was seeded in a library of 5909 random molecules retrieved from the ZINC database [39]. The random molecules had chemical properties (such as molecular weight and formal charge), similar to the known PKR antagonists, to ensure that the PLoS ONE | www.plosone.org 5 November 2011 | Volume 6 | Issue 11 | e27990