Three crystal structures of pharmaceutically relevant protein target, HER2 in complex with their co-crystal ligands i.e. inhibitors (Scheme 2), having a resolution of less than 2.5Ǻ (PDB-IDs; 3W33 (1.70Å), 3W32 (1.80Å) and 3PP0 (2.25Å)), were considered in the study [13]. Out of these protein-ligand complexes, 3W32 protein-ligand complex was selected on the basis of crystal resolution, re-docking score and RMSD value, for the generation energy optimized structure based pharmacophore model.

Scheme 2. Structure of HER2 co-crystal ligands with their respective PDB-IDs.

Schrödinger offers the Phase Commercially Available Compound Database (CACDB), containing unique structural records of 612551 molecules. Also, freely available NCI-Open-2012 (265242 molecules as anticancer agents) database was downloaded and used for the VS against an E-pharmacophore model. A widely accepted benchmark dataset of decoys having molecular weight 400Da was downloaded from the Schrödinger’s website [14]. This decoy set consists of 1000 molecules with similar properties but dissimilar topology to the active compounds.

Coordinates for each crystal structures were taken from the RCSB Protein Data Bank (PDB) [15] and prepared using the Protein Preparation Wizard, located in Maestro Software Package [16]. Bond orders and formal charges were added for hetero-groups, and hydrogens were added to all atoms in the system. To optimize the hydrogen bond network, His tautomers and ionization states were predicted, 180° rotations of the terminal angle of Asn, Gln, and His residues were assigned, and hydroxyl and thiol hydrogens were sampled. Water molecules in all structures were removed as none of them was found to establish a stable interaction with the protein and the inhibitor. For structures with missing side-chain atoms, the refinement module in Prime [17] was used to predict their conformations. For each structure, a brief relaxation was performed using an all-atom constrained minimization carried out with the Impact Refinement module [18] using the OPLS-2005 force field to alleviate steric clashes that may exist in the original PDB structures. The minimization was terminated when the energy converged or the RMSD reached a maximum cutoff of 0.30 Å.

All the screened hits were prepared by Ligprep module using Epik in Schrödinger package [19] to expand protonation and tautomeric states at 7.0±2.0 pH units, while bioactive conformations of co-crystal ligands were used in the study. Conformational sampling was performed on all database molecules using the ConfGen search algorithm. Confgen with OPLS 2005 force field was applied for generation of conformers. The duplicate poses were eliminated if the RMSD was less than 1.0 Å. A distance dependent dielectric constant of 4 and maximum relative energy difference of 10 kcal/mol was applied as suggested by Salam et.al [20].

Glide energy grids were generated for each of the prepared complexes, and the generated grid file from the prepared receptor was used for docking calculations considering the ligands as flexible but treating the receptor as a rigid structure. The binding site was defined by a rectangular box surrounding the X-ray ligand. Ligands were docked into their respective binding sites in the protein crystal structure. The “Glide XP” protocol [21] was chosen during the docking run which deduces energy terms such as hydrogen-bond interactions, electrostatic interaction, hydrophobic enclosure, and pi–pi stacking interaction and the rest of the parameters were kept at default for the scoring [22].

Energy-optimized structure based pharmacophores, i.e. E-pharmacophores were generated through docking, post-processing, E-pharmacophore option situated in scripts menu bar of Maestro Software Package. Glide_XP output (Xpdes (mae.pv)) file, obtained after docking of selected 3W32 HER2 co-crystal ligand to its native protein crystal structure, was used to generate the E-pharmacophore. Pharmacophore sites were automatically generated with Phase [23], using the default set of six chemical features: hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R). Each pharmacophore feature site was first assigned an energetic value equal to the sum of the Glide XP contributions of the atoms forming the site. This allows the sites to be quantified and ranked on the basis of these energetic terms involved in their bioactivity towards the target protein.

Three methods of VS protocol were applied in the current study; E-pharmacophore search using Phase [24], docking using Glide [22, 25], and 3D shape similarity search using Phase [26].

After running Phase find matches, a total of 1608 hits were retrieved from ASINEX and NCI databases, which were prepared by LigPrep, and filtered applying Lipinski’s rule and ADME (QikProp) filter criteria. Phase “shape screening” protocol was applied to find out shape similarity of these hits. The active compound 3W32 co-crystal ligand having six pharmacophoric features, along with excluded volume, was selected to create the 3D shape query as the same compound was used to generate the E-pharmacophore in the Phase find matches [26, 27].

Many metrics are currently used to evaluate the performance of ranking methods in VS, for instance, the average rank of actives, the enrichment factor (EF), the area under the receiver operating characteristic curve (ROC), Boltzmann-enhanced discrimination of receiver operating characteristic (BEDROC), the area under the accumulation curve (AUAC), the robust initial enhancement (RIE), and were used to determine the robustness of the hypothesis. Evaluating the performance of VS methods is useful to select the method which performs best in a given situation eventually, retrieving active compounds out of a mixed set of active compounds and decoys (compounds presumably inactive against the examined target) [28].

The area under the ROC curve is used by influential groups to measure VS performance in part because it does possess desirable statistical behavior. A value of 1/2 shows that the ranking method does not do better than random picking. It can be interpreted as the probability that an active will be ranked before an inactive. It has a value between 0 (worst performance attainable) and 1 (best performance) [25, 28, 29].

Enrichment factor (EF) is defined as the ratio of the probabilities of searching an active compound in top X% of the data set [30]. A maximum enrichment is therefore 100 if all of the actives (A=1) are found within the top 1% of the decoys (D=0.01).

Where, H_a: Number of actives in the hit list (true positives); H_t: Number of the hits retrieved; A: Number of active molecules in the database; D: Total number of database compounds; % A: The ratio of the actives retrieved in the hit list (precision); % Y: The yield of actives (recall); EF: Enrichment factor (i.e. enrichment of the concentration of the actives by the model relative to random screening without any pharmacophoric approach); GH: The Guner-Henry score.

Fig. 1. Superposition of the co-crystal ligand 3W32 (Grey) with its docked pose (Sky blue).

Fig. 2a. Energy-optimized pharmacophore (E-pharmacophore) hypotheses ADHRRR (Pink sphere/circle: hydrogen bond acceptor, green sphere/circle: hydrophobic group, orange ring: aromatic ring, light-blue: hydrogen bond donors).

Fig. 2b. Site measurement for the selected E-pharmacophore hypothesis ADHRRR (Pink sphere/circle: hydrogen bond acceptor, green sphere/circle: hydrophobic group, orange ring: aromatic ring, light-blue: hydrogen bond donors).

Figure 3a. ROC Plot between 1-specificity and sensitivity.

Figure 3b. ROC Plot between percent screen and percent actives found.

Fig. 4a. 2D Interactions of top scoring screened hit 355682-ASINEX at the binding site of HER2 protein, 3W32.

Fig. 4b. 3D Interactions (H-bond; Purple colour dotted lines) of top scoring screened hit 355682-ASINEX at the binding site of HER2 protein, 3W32.

Fig. 5a. 2D Interactions of top scoring combinatorial hit 355682-ASINEX_1 at the binding site of HER2 protein 3W32.

Fig. 5b. 3D Interactions (H-bond; Purple colour dotted lines) of top scoring combinatorial hit 355682-ASINEX_1 at the binding site of HER2 protein, 3W32.

Despite the early recognition problem, the EF has some problems ignoring complete ranking of the whole dataset molecules [30]. A method that is superior to random selection of compounds has EF > 1. To address the problem of EF, as discussed in various reports [31], Sheridan et al. [26] developed an exponential weighted scoring scheme RIE which gives heavier weight in “early recognized” hits. A second enrichment metrics, the BEDROC [28], was also used as a way to ensure that the results and conclusions were significant. BEDROC is a generalization of the ROC that addresses the “early scoring problem” by Boltzmann weighting the hits based on how early they are retrieved. BEDROC can be defined as the probability that an active is ranked before a randomly selected compound was exponentially distributed with parameter α and its values ranges in between 0-1. The value of α= 20.0 is suggested as a reasonable choice for virtual screening evaluations and corresponds to 80% of the total score being accounted for in the top 8% of the database. BEDROC and RIE have a linear relationship.

For the validation of the screening methods, 3 co-crystal ligands, i.e.. inhibitors, present in the crystal structures of HER2 proteins: PDB-IDs; 3W33, 3W32 and 3PP0 (Scheme 2), 23 combinatorial hit molecules (Table 5) and 1000 decoys making a total of 1026 compounds were used in the study. Phase find matches and Glide XP docking was applied selecting these 1026 molecules as input. The output files of these screening methods were chosen to calculate the enrichments using the above said 26 molecules as actives.

The hit molecules obtained from the above screening methods were submitted to Canvas to generate binary fingerprints (dendritic: linear and branched), using Tanimoto Coefficient (0≤T≤1), to measure the similarity between two fingerprints. The resulted binary fingerprints were used for the clustering of the obtained hits based on their structural similarity. A clustering (hierarchical: non overlapping clustering method where cluster size increases from single molecule per cluster) algorithm divides a group of objects into clusters where molecules within each cluster are more closely related to one another than objects assigned to different clusters. The screened hits (1061) were differentiated into 15 clusters and sorted on the basis of fingerprints and docking XP gscores. The best scoring seven cluster representative molecules were selected, which represent a comprehensive structural diversity amongst the resulted hit molecules.

Nearly 40% of drug candidates fail in clinical trials due to poor ADME (absorption, distribution, metabolism, and excretion) properties. These late-stage failures contribute significantly to the rapidly escalating cost of new drug development. The ability to detect problematic candidates early can dramatically reduce the amount of wasted time and resources, and streamline the overall development process. QikProp [32] program was used for in silico prediction of pharmacokinetic properties of the seven screened hit molecules.

After binary fingerprinting and hierarchical clustering, the selected best scoring cluster representative was submitted to Combiglide for maximizing the number of hits. A combinatorial library was generated using “Interactive Enumeration and Docking” option present in the Combiglide module of the Schrödinger Software. First, a combinatorial library was created with the help of best scoring template molecule obtained from VS through R group variation, and then it was docked into the crystal structure of the same protein without affecting the previous binding pose. The molecules obtained after combinatorial designing, having better XP gscores than the template molecule, were selected for further studies.

Docking post-processing, interaction fingerprint option situated in scripts menu bar of Maestro user interface was

used to analyze the structural interactions of the obtained hits (retrieved hits and Combiglide generated combinatorial hits) with the HER2 protein crystal structure.

The main purpose of virtual screening is to find a novel scaffold and the possible lead compounds suitable for further development. A combination of structure and ligand based approaches, in association with ligand similarity (shape based) practices, were applied in a stepwise filtering approach to obtain the potential hits.

It is a hybrid approach of ligand and a structure-based technique which uses docking energy score for finding the bioactive component of ligands against the receptor. Energy optimized six site (one HBdonar, one HBacceptor, one HPhobic, and three ring aromatics) hypothesis (Table 2, Fig. 2) was generated and submitted as a query to retrieve the potential hits from ASINEX and NCI databases. Entire 1608 retrieved hits were subjected to filtering by applying Lipinski’s rule and ADME properties. Database hits were ranked in order of their fitness score (a measure of how well the aligned ligand conformer matches the hypothesis based on RMSD site matching, vector alignments, and volume terms) and selected, having fitness value ≥1. The fitness scoring function is an equally weighted composite of these three terms and ranges from 0 to 3, as implemented in default database screening in Phase [20].

The screening validation of the VS process is needed to check out wether the used method is efficient in retrieving the actives from the databases and ranking them early or not. Different validation parameters were calculated in this study and the results are presented in Table 3. The enrichment results showed that the used screening protocol was satisfactory in retrieving the active compounds from the molecular databases. Three out of the five ranked actives, i.e. in more than 50% yield, were retrieved at 1% of the results with remarkable ROC (0.98), AUAC (0.98), RIE (11.66), GH (0.43) values and a high EF of 50 (at1% of sample size). The ROC plot (Fig. 3a) and % screen plot ROC plot (Fig. 3b) demonstrates that the used method was sensitive and specific in recognizing the active molecules. A substantial value of BEDROC at alpha=20.0, and alpha= 160.9, shows the early recognition of the actives amongst the database compounds.

The six structural features (A1, D7, H12, R13, R14, and R15) of 3W32 HER2 protein co-crystal ligand, including excluded volume, were taken as a 3D-shape query. The shape similarity of seven good scoring compounds is presented in the Table 4, which indicate towards the structural diversity amongst the retrieved hits.

The main aim of docking study is to find the binding affinity between protein–ligand complexes. The hit molecules, retrieved from the virtual screening technique, were further refined by molecular docking using Glide XP protocol to check whether these chemical features mapped with structure-based interaction mode or not. A total of 1061 binding poses were obtained from 1608 screened molecules, which were sorted on the basis of XP gscores (Table 4).

The fingerprints encode the absence or presence of a set of structural fragments (or chemical features) which can be useful in clustering of molecules, diversity analysis and selection of molecules, QSAR model building and similarity search of a database. Out of 15 clusters, 7 structurally diverse representative molecules (Scheme 3) with highest XP gcores (>-10) among the clusters were selected. The structural diversity of the retrieved hits emphasizes the scaffold hoping capability of the E-pharmacophore ADHRRR, when used in the mining of databases of drug like compounds. Also, it has been observed that the entire seven screened hit molecules obtained, satisfy the prerequisites of ADME properties (Table 5).

Scheme 3. 2D Structure of seven selected top scoring screened hit molecules.

The top scoring hit molecule obtained from ASINEX database (355682-ASINEX) was submitted to Combiglide for the interactive enumeration and docking. Variation in R group resulted in the generation of a combinatorial library of more than 10000 combinatorial hits. Out of these, 22 molecules (355682-ASINEX_01 to 355682-ASINEX_22), obtained from varying R₁ and R₂ substituents, were found to show better binding interactions and XP gscore than the template molecule (355682-ASINEX). The results are summarized in Table 6.

Table 6. Top scoring combinatorial hits (Combiglide generated hits) obtained through R group variation.

Comps. Code	Hits obtained from R group variation (Combiglide hits)				Interaction Fingerprints
	R1	R2	Structure	XP gscore	SIFT Any Contact	SIFT Sidechain Interaction	SIFT HBaccept
355682- ASINEX (Template Molecule)	-H	-H		-12.9541	2	2	0
355682-  ASINEX_1	-NH2	-C6H5		-14.230	3	3	0
355682-  ASINEX_2	-NH2	-COCH3		-14.108	3	3	0
355682-  ASINEX_3	-OH	-CSCH3		-13.948	3	3	0
355682-  ASINEX_4	-OH	-COCH3		-13.926	2	2	0
355682-  ASINEX_5	-Cl	-OCH3		-13.913	3	3	0
355682-  ASINEX_6	-OCH3	-CONH2		-13.790	3	3	0
355682-  ASINEX_7	-OH	-OCH3		-13.720	2	2	0
355682-  ASINEX_8	-CONH2	-C2H5		-13.706	3	3	1
355682-  ASINEX_9	-NH2	-OCH3		-13.705	3	3	1
355682-  ASINEX_10	-CH3	-H		-13.698	3	3	0
355682-  ASINEX_11	-Cl	-OH		-13.615	2	2	0
355682-  ASINEX_12	-OCH3	-CH3		-13.588	2	2	0
355682-  ASINEX_13	-CH3	-CSCH3		-13.558	3	3	0
355682-  ASINEX_14	-NH2	-OH		-13.527	2	2	0
355682-  ASINEX_15	-CONH2	-OCH3		-13.509	2	2	0
355682-  ASINEX_16	-CONH2	-CSCH3		-13.432	2	2	0
355682-  ASINEX_17	-Cl	-C6H5		-13.329	2	2	0
355682-  ASINEX_18	-C2H5	-C6H5		-13.234	2	2	0
355682-  ASINEX_19	-C2H5	-CSCH3		-13.217	2	2	0
355682-  ASINEX_20	-NH2	-H		-13.115	3	3	0
355682-  ASINEX_21	-NH2	-CONH2		-13.052	3	3	0
355682-  ASINEX_22	-NH2	-C2H5		-13.019	3	3	0

The selected hit molecules were further analyzed by visual inspection and SIFT analysis (Tables 4 and 6) to find out orientation and critical interactions between the receptor and ligand. All the seven screened hit molecules have shown XP gscores of >-10. The best scoring hit molecule, 355682-ASINEX interacts (Fig. 4) with essential amino acid residue Met 793 (backbone), similar to that of the 3W32 co-crystal ligand in the complex with the HER2 protein crystal structure. This molecule was also found to be involved in two additional H-bond interactions with Gln791 (backbone) and Thr854 (side chain) amino acid residues. The best scoring combinatorial hit molecule, 355682-ASINEX_1 generated through Combiglide program, showed the highest XP gscore of -14.230, which was comparable to the binding score of 3W32 HER2 protein co-crystal ligand (-14.874). Furthermore, this molecule showed two new H-bond interactions with Arg841 and Phe856 (Fig. 5) in addition to the above interactions.