compound library will bind to the target’s binding site. When a ligand-receptor complex is available, either from an X-ray structure or an experimentally verified model, a structure-based pharmacophore model describing the possible interaction points between the ligand and the receptor can be generated using different algorithms and later used for screening compound libraries [27]. In ligand-based VLS procedures, the pharmaco- phore is generated via superposition of 3D structures of several known active ligands, followed by extracting the common chemical features responsible for their biological activity. This approach is often used when no reliable structure of the target is available [28]. In this study, we analyzed known active small-molecule antagonists of hPKRs vs. inactive compounds to derive ligand- based pharmacophore models. The resulting highly selective pharmacophore model was used in a VLS procedure to identify potential hPKR binders from the DrugBank database. The interactions of both known and predicted binders with the modeled 3D structure of the receptor were analyzed and compared with available data on other GPCR-ligand complexes. This supports the feasibility of binding in the TM-bundle and provides testable hypotheses regarding interacting residues. The potential cross-reactivity of the predicted binders with the hPKRs was discussed in light of prospective ’off-target’ effects. The challenges and possible venues for identifying subtype-specific binders are addressed in the discussion section. Materials and Methods Homology Modeling and Refinement All-atom homology models of human PKR1 and PKR2 were generated using the I-TASSER server [29], which employs a fragment-based method. Here a hierarchical approach to protein structure modeling is used in which fragments are excised from multiple template structures and reassembled, based on threading alignments. Sequence alignment of modeled receptor subtypes and the structural templates were generated by the TCoffee server [30]; this information is available in the Supporting Information as figure S1. A total of 5 models per receptor subtype were obtained. The model with the highest C-score (a confidence score calculated by I-Tasser) for each receptor subtype, was exported to Discovery Studio 2.5 (DS2.5; Accelrys, Inc.) for further refinement. In DS2.5, the model quality was assessed using the protein report tool, and the models were further refined by energy minimization using the CHARMM force field [31]. The models were then subjected to side-chain refinement using the SCWRL4 program [32], and to an additional round of energy minimization using the Smart Minimizer algorithm, as implemented in DS2.5. The resulting models were visually inspected to ensure that the side chains of the most conserved residues in each helix are aligned to the templates. An example of these structural alignments appears in figure S2. For validation purposes, we also generated homology models of the turkey b1 adrenergic receptor (b1adr) and the human b2 adrenergic receptor (b2adr). The b1adr homology model is based on 4 different b2adr crystal structures (PDB codes – 3SN6, 2RH1, 3NY8, and 3d4S); the b2adr model is based on the crystal structures of b1adr (2VT4, 2YCW), the Dopamine D3 receptor (3PBL), and the histamine H1 receptor (3RZE). The models were subjected to the same refinement procedure as previously described, namely, deletion of loops, energy minimization, and side chain refinement, followed by an additional step of energy minimization. Sometimes the side chain rotamers were manually adjusted, following the aforementioned refinement procedure. Small-Molecule Binders of Prokineticin Receptors Throughout this article, receptor residues are referred to by their one-letter code, followed by their full sequence number in hPKR1. TM residues also have a superscript numbering system according to Ballesteros-Weinstein numbering [33]; the most conserved residue in a given TM is assigned the index X.50, where X is the TM number, and the remaining residues are numbered relative to this position. Identification of a 7TM-bundle binding site The location of a potential small-molecule-TM binding cavity was identified based on (1) identification of receptor cavities using the "eraser" and "flood-filling" algorithms [34], as implemented in DS2.5 and (2) use of two energy-based methods that locate energetically favorable binding sites – Q-SiteFinder [35], an algorithm that uses the interaction energy between the protein and a simple Van der Waals probe to locate energetically favorable binding sites, and SiteHound [36], which uses a carbon probe to similarly identify regions of the protein characterized by favorable interactions. A common site that encompasses the results from the latter two methods was determined as the TM-bundle binding site for small molecules. SAR Analysis A dataset of 107 small-molecule hPKR antagonists was assembled from the literature [18,19]. All ligands were built using DS2.5. pKa values were calculated for each ionazable moiety on each ligand, to determine whether the ligand would be charged and which atom would be protonated at a biological pH of 7.5. All ligands were then subjected to the "Prepare Ligands" protocol, to generate tautomers and enantiomers, and to set standard formal charges. For the SAR study, the dataset was divided into two parts: (1) active molecules, with IC50 values below 0.05 mM, and (2) inactive molecules, with IC50 values above 1 mM. IC50 values were measured in the calcium mobilization assay [18,19]. When possible, the molecules were divided into pairs of active and inactive molecules that differ in only one chemical group, and all possible pharmacophore features were computed using the "Feature mapping" protocol (DS 2.5). These pairs were then compared to determine those pharmacophore features’ impor- tance for biological activity. Ligand-Based Pharmacophore Models The HipHop algorithm [37], implemented in DS2.5, was used for constructing ligand-based pharmacophore models. This algorithm derives common features of pharmacophore models using information from a set of active compounds. The two most active hPKR antagonists (the lowest IC50 values in the Janssen patent [19,20]) were selected as ‘reference compounds’ from the data set described above, and an additional antagonist molecule with a different scaffold was added from a dataset recently published [38], and were used to generate the models (figure S3). Ten models in total were generated, presenting different combinations of chemical features. These models were first evaluated by their ability to successfully recapture all known active hPKR antagonists. An enrichment study was performed to evaluate the pharmacophore models. The dataset contains 56 active PKR antagonists seeded in a random library of 5909 decoys retrieved from the ZINC database [39]. The decoys were selected so that they will have general and chemical properties similar to the known hPKR antagonists (by filtering the ZINC database according to the average molecular properties of known hPKR antagonists 6 4 Standard Deviation range). In this way, enrichment is not simply achieved by separating trivial features PLoS ONE | www.plosone.org 3 November 2011 | Volume 6 | Issue 11 | e27990