Potent bace-1 inhibitor design using pharmacophore modeling, in silico screening and molecular docking studies

Dataset collection

In a computerized pharmacophore generation process the accurate choice of the training set is a key issue. The built pharmacophore hypothesis can be as good as the input data information. The following criteria should be considered during the selection of data set in order to achieve a significant pharmacophore hypothesis. (1) All compounds used in the training set have to bind to the same receptor in roughly the same fashion. Compounds having more binding interaction with the receptor are more active than those with fewer; (2) the data set must be widely populated covering an activity range of at least 4 orders of magnitude; (3) the most active compounds should inevitably be included in the training set and (4) all biologically relevant data should be obtained by homogenous procedures [11]. Every individual feature in the resulting hypotheses will invade a certain weight that is proportional to its relative contribution to biological activity.

By taking these criteria into account, we have collected a total number of 320 BACE-1 inhibitors from various literature resources [12–25] and a database has been created. The 2D structure of the compounds were built using ChemSketch program version 12, and subsequently converted in to 3D structures using Discovery studio 2.5 (DS) [26]. In the next step, 60 compounds were selected as final dataset as the BACE-1 inhibitory activities of these 60 compounds were studied under same biological assay condition. Based on the principle of structural diversity and wide coverage of activity range, 20 compounds were carefully selected as training set compounds and the rest were used as test set in model validation. Here, the IC₅₀ value of the training set compounds was taken into account, the inhibitory activity values of the training set compounds span over a range of four orders of magnitude, from 4 nM to 37 000 nM. The chemical structure and experimental activity of the training set compounds are shown in Figure 1.

Diverse conformation generation

Prior to the generation of pharmacophore hypotheses, the training set compounds, which were converted to 3D structure, were used to generate diverse conformations. Diverse Conformation Generation protocol implemented in DS was used to generate conformations using the Best conformation model generation method with CHARMM force field and Poling algorithm to ensure the energy-minimized conformation for each compound. The parameters like maximum number of 250 conformers, the ‘best conformational analysis’ method, and an energy threshold of 20 kcal/mol above the global energy minimum were chosen during conformation generation.

Pharmacophore modeling

The training set comprises of 20 compounds was used in pharmacophore hypothesis generation. The HypoGen algorithm available in 3D QSAR Pharmacophore Generation protocol of DS tries to generate hypotheses with features common amongst active molecules and do not reflect the inactive molecules of the training set. The training set compounds were predicted for their inherent chemical features using Feature Mapping protocol implemented in DS. During pharmacophore hypothesis generation a minimum of 4 and a maximum of 5 pharmacophoric features like hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), positive ionizable (PI), ring aromatic (RA) and hydrophobic (HY) were included. These features were selected based on the feature mapping results. All parameters were set to their default values except uncertainty value, which has been changed to 2 instead of 3. An uncertainty value of 2 was more convenient for our dataset because the activity values of the training set spanned exactly the required 4 orders of magnitude; this choice has been further confirmed by preliminary calculations and by other literature evidence [27]. The uncertainty value represents the ratio of the uncertainty range of the actual activity against measured biology activity for each compound.

The HypoGen algorithm

With the full range of training set compounds from active to inactive the pharmacophore hypotheses were generated by HypoGen algorithm implemented in DS. This algorithm constructs and ranks the pharmacophore hypotheses that correlate best between 3D spatial arrangement of features in a given training set compounds and their respective experimental activities. This process is accomplished in three steps: the constructive phase, the subtractive phase and the optimization phase [28].

The constructive phase identifies the hypotheses that are common amongst the active compounds. HypoGen enumerates all possible pharmacophore configurations using all combinations of pharmacophore features for each of the conformations of the most active compound. In order to consider the left over most active compounds the hypotheses must fit a minimum subset of its features. Hence, a large database of pharmacophore configurations will be generated at the end of the constructive phase.

The subtractive phase will remove the pharmacophore configurations that are present in the least active compounds. All compounds whose activity is by default 3.5 orders of magnitude less than that of the most active compound are considered to represent the least active compounds. The value 3.5 is adjustable depending on the activity range of the training set. The optimization phase improves the hypothesis score. These scores of the generated hypotheses depend on the errors in activity estimation from regression and complexity. The optimization involves a variation of features and/or locations to optimize activity prediction via a simulated annealing approach. The total cost parameter will be calculated for every new hypothesis. The HypoGen will quit and reports the 10 top-scoring hypotheses when there is no improvement in the hypothesis score.

Data analysis

The quality of a pharmacophore hypothesis is best determined by two theoretical cost calculations, which are represented in bit units [29]. One is the ‘fixed cost’ also known as cost of an ideal hypothesis, which represents the simplest model that fits all the data perfectly. The second cost is the ‘null cost’, which represents the highest cost of a pharmacophore with no features that estimates every activity to be the average of the activity data of the training set compounds.

The total cost of any pharmacophore hypothesis should always be close to the fixed cost and away from the null cost to be the significant model. The cost difference between fixed and null cost values should be larger for a meaningful pharmacophore hypothesis. A value of 40-60 bits in a pharmacophore hypothesis indicates that it has 75-90% probability of representing a true correlation in the data.

The hypotheses are also evaluated based on other cost components. The cost value for every individual hypothesis is the summation of three cost components: the error cost (E), the weight cost (W) and the configuration cost (C). The error cost is the value represents the root-mean-squared difference (RMSD) between experimental and estimated activity value of the training set compounds. The weight cost is a value that increases in a Gaussian form as this function weights in a model deviate from the ideal value of two. The configuration cost or entropy cost measures the entropy of the hypothesis space. If the input training set compounds are too multiplex, e.g. because of too flexible training set compounds, this will result in an effusive number of hypotheses as an outcome of the subtractive phase. This configuration cost should always be less than a maximum value of 17 [30]. The correlation coefficient of the pharmacophore hypothesis should be close to 1.

Pharmacophore validation

The generated pharmacophore hypothesis was validated using test set, Fischer randomization, decoy set and leave-one-out methods.

Database screening

The validated pharmacophore hypothesis, Hypo 1, was used as a 3D query for screening four different chemical databases. The purpose of this screening is to retrieve novel and potential leads suitable for further development. The chemical databases used were Maybridge, Chembridge, NCI and Asinex. Conformers were generated for each molecule in the database using best conformer generation method that allows a maximum energy of 15 kcal/mol above that of the most stable conformation. The database screening was carried out using Ligand Pharmacophore Mapping protocol implemented in DS with Best/Flexible Search option. The retrieved compounds were filtered by restricting the estimated activity value less than 100 nM and the obtained compounds were further refined using molecular docking study.

Molecular docking

Pharmacophore modeling normally firmly associated with docking procedure, which in a first step flexibly aligns the ligand molecule into a rigid macromolecule environment and then estimates the tightness of the interaction by different scoring functions [35]. The Docking takes all the information from a rigid protein environment and scores several possible interaction modes for different alignments. There are many docking programs available for molecular docking studies. In this study, we used GOLD (Genetic Optimisation for Ligand Docking), a docking program [36] that uses genetic algorithm for docking and performs automated docking with full acyclic ligand flexibility, partial cyclic ligand flexibility and partial protein flexibility in the neighborhood of the protein active site. The crystal structure of BACE-1 complexed with an inhibitor SC7 (PDB ID: 2QP8) was used in molecular docking studies. The inhibitor SC7 was extracted from the active site and the retrieved database hits were docked based on the ligand SC7 coordinates, in to the active site of BACE-1. The water molecules were removed prior to docking because they were not found to play any important roles in BACE1-ligand interaction. The early termination option parameter in GOLD was changed from 3 to 5 and the maximum save conformations was set to 10. All the other parameters were set at their default values.

Pharmacophore generation

Statistical data analysis

The generated hypotheses were subjected to cost analysis. The two main values used for cost analysis are the difference between fixed and null cost and another is the difference between the null cost and the total cost (Δcost). The fixed cost of the run was 74.77 bits, which was well separated from the null cost of 203.22 bits and close to the total cost of 81.24 bits. The large difference (128.45 bits) observed between the fixed cost and null cost value indicates that Hypo 1 has more than 90% statistical significance to be a significant model. All the 10 hypotheses were subjected to further evaluation for their capability to predict the activity of the training set compounds. Configuration cost or entropy value must be less than 17 for which a value of 15.59 was obtained in this study. All hypotheses have scored RMSD values lower than 1.5 and ranging from 0.804 to 1.111, this characterization further emphasizing the good predictive quality of these hypotheses. Based on the rule to select a hypothesis with a lowest total cost, high correlation coefficient, large cost difference and significantly low RMSD value, Hypo 1 gave the best statistical values among other hypotheses. Hence, Hypo 1 with one HBD, one PI, one RA and two HY was chosen as the best hypothesis for further analysis. The inter-feature distance constraints were observed for this five-featured pharmacophore hypothesis (Hypo 1) (Figure 2).

Activity prediction and mapping of training set compound on Hypo1

To verify the predictive ability of Hypo1 on training set compounds, the activity of each training set compound is estimated by regression analysis. The experimental activities of training set compounds were classified into four groups: most active (IC₅₀ ≤ 100 nM, ++++), active (100 nM < IC₅₀ ≤ 1000 nM, +++), moderately active (1000 nM < IC₅₀ ≤ 10 000 nM, ++) and inactive (IC₅₀ > 10 000 nM, +). The estimated activity values of training set compounds based on Hypo 1 and the corresponding error values are calculated (Table 2). The error value is the ratio between the estimated and experimental activities. The positive error value indicates that the estimated IC₅₀ value is higher than the experimental activity, whereas the negative error value indicates that the estimated IC₅₀ value is lower than the experimental activity. An error value of less than 10 signifies the prediction of activity lesser than one order of magnitude. Among 20 training set compounds, only one compound had an error value of greater than 3. From Table 2 it is clear that the estimated activity values of most of the training set compounds was predicted with the same activity scale as the experimental activity. Among 20 training set compounds, one most active compound (++++) was estimated as active (+++), one active compound (+++) was estimated as moderately active (++), one moderately active compound (++) was estimated as inactive (+) and two inactive compounds (+) were estimated as moderately active (++). The divergence between the estimated and experimental activity observed in four compounds was only about 1 order of magnitude, which might be an artifact of the program that uses different number of degrees of freedom for these compounds to mismatch the pharmacophore model. Interestingly, for feature fitting, the most active compounds in the training set mapped well on all the chemical features that are one HBD, one PI, one RA and two HY features of Hypo 1 with good fitting score. The active, moderately and inactive compounds have missed at least one of five features. In addition, the most active compounds mapped well on the RA and PI features whereas some active, moderately active and all inactive compounds could not map on the RA and PI features signifying the importance of these two features. The pharmacophore overlay of most active Compound 1 and Hypo 1 has shown a fit value of 9.52. The RA feature corresponds to phenyl ring present in between two amide and a sulfonamide groups, one HBD feature corresponds to nitrogen of amide group located at the branch, PI group corresponds to the only amino nitrogen, one HY feature corresponds to a phenyl group whereas the another HY feature corresponds to alkyl group (Figure 3A). The pharmacophore overlay of least active compound 20 has revealed that it missed two features when mapped on Hypo 1 with a fit value of 5.16. This compound has mapped only the HBD and two HY features in the same manner as most active compound with no mapping over RA and PI features (Figure 3B). Fit value indicates how well the features in the pharmacophore overlaps the chemical features in the molecule and thereby aid in understanding the chemical meaning of the hypothesis [37]. These results emphasized Hypo 1 as a reliable model to accurately estimate the experimental activity of the training set compounds.

Compound	IC50 nM		Errora	Fit valueb	Activity scalec
Experimental	Estimated	Experimental	Estimated
1	4	1.7	-2.4	9.52	++++	++++
2	7	16	+2.3	8.54	++++	++++
3	13	13	-1	8.64	++++	++++
4	52	88	+1.7	7.80	++++	++++
5	90	150	+1.7	7.57	++++	+++
6	115	170	+1.5	7.51	+++	+++
7	240	190	-1.3	7.48	+++	+++
8	348	320	-1.1	7.24	+++	+++
9	450	470	+1.1	7.07	+++	+++
10	674	460	-1.5	7.08	+++	+++
11	730	610	-1.2	6.97	+++	+++
12	980	1400	+1.4	6.60	+++	++
13	1500	880	-1.7	6.80	++	+
14	2800	8100	+2.9	5.84	++	++
15	4600	5200	+1.1	6.03	++	++
16	8000	10000	+1.3	5.73	++	++
17	12000	8000	-1.6	5.84	+	++
18	19000	19000	-1	5.48	+	+
19	25000	6100	-4.1	5.96	+	++
20	37000	38000	+1	5.16	+	+

^aPositive value indicates that the estimated IC₅₀ is higher than the experimental IC₅₀; negative value indicates that the estimated IC₅₀ is lower than the experimental IC₅₀.

^bFit value indicates how well the features in the pharmacophore map the chemical features in the compound.

^cActivity scale: most active, ++++, IC₅₀ ≤ 100 nM; active, +++, 100 nM < IC₅₀ ≤ 1000 nM; moderately active, ++, 1000 nM < IC₅₀ ≤ 10,000 nM; inactive, +, IC₅₀ > 10,000 nM.

Validation of Hypo 1

Hypo 1 was further validated by test set, Fischer randomization test, leave-one-out and decoy set methods.

Test set method

A total of 40 compounds structurally different from the training set compounds were selected as test set. The test set compounds were prepared in the same way training set compounds were prepared. The top-scored 10 hypotheses was regressed against 40 test set compounds and calculated the correlation coefficient values ranging from 0.917 to 0.859 (Table 1) between experimental and estimated activities. Among 10 hypotheses, Hypo 1 has given a correlation coefficient of 0.917 (Figure 4) indicating a good correlation between the estimated and experimental activities. The predictive ability of the Hypo 1 against test set compounds was considered better than other hypotheses and the estimated activity values along with the experimental and error values based on Hypo 1 are tabulated (See additional file 1: Experimental and estimated IC50 values of the test set compounds based on the pharmacophore hypothesis ‘Hypo 1’.). Most of the test set compounds was estimated correctly to their experimental activity. The test compounds were classified into four groups in a similar way as that of training set: most active (IC₅₀ ≤ 100 nM, ++++), active (100 nM < IC₅₀ ≤ 1000 nM, +++), moderately active (1000 nM < IC₅₀ ≤ 10 000 nM, ++) and inactive (IC₅₀ > 10 000 nM, +). A total of 11 out of 12 active (+++) compounds were estimated correctly as active, but 1 compound was estimated as most active (++++). Interestingly all the six most active (++++) compounds were estimated correctly as most active (++++). Total of twelve active (+++) compounds were estimated correctly as active. Out of thirteen moderately active (++) compounds only one compound was over estimated as active (+++), and 3 compounds were under estimated as inactive (+), whereas nine compounds were estimated correctly as moderately active. Among nine inactive (+) compounds one was over estimated as moderately active (++) compound whereas eight were estimated as inactive compounds. A total of 2 inactive (+) compounds out of 9 were over estimated as moderately active (++) whereas 7 were correctly estimated as inactive compounds. These results suggested that Hypo 1 has a good agreement with the experimental data and able to predict the activities of a wide variety of BACE-1 inhibitors.

Database screening

The validated pharmacophore hypothesis, Hypo1, was used as a 3D structural query for retrieving compounds from chemical databases including MayBridge (59 652 compounds), Chembridge (50 000 compounds), NCI (238 819 compounds) and Asinex (213 462 compounds). As a result of first screening 11 578, 590, 5096 and 63 265 compounds were retrieved from Maybridge, Chembridge, NCI and Asinex respectively. Since the active site of BACE-1 is larger in size, the experimentally known most active inhibitors are also larger in size and violate the first rule of Lipinski’s rule of five. Hence, the retrieved hit compounds were filtered based only on the estimated activity values calculated by Hypo 1. The activity range for the most active compounds is <100 nM. Finally 773 compounds were selected by restricting the minimum estimated activity to <100 nM.

Molecular docking

To further refine the retrieved hits and also to remove the false positives, these 773 compounds along with the 20 training set compounds were docked into the active site of BACE-1 using GOLD 4.1 program. There are number of crystal structures for BACE-ligand complexes are available in PDB. The crystal structure of BACE-SC7 (PDB ID: 2QP8) complex was taken based on its high resolution. The GOLD fitness score was calculated for all the 793 compounds, it distinguishes molecules based on their interacting ability. The GOLD fitness score for the most active compound in the training set was 53.035. The compounds for further analysis were selected based on the ligand conformations which can satisfy the necessary interactions at the active site and scoring GOLD fitness score greater than 60. Finally 20 compounds from Maybridge and 15 compounds from Asinex have shown the required interaction with BACE-1 as well as good GOLD fitness scores. The compounds with the same chemical scaffolds were filtered carefully based on the molecular interactions observed at the active site. Finally, two compounds with different scaffolds one from Maybridge (RJC01726) and one from Asinex (Asnx-2) were selected as representative compounds. The binding mode of the final hits and the most active Compound 1 in the training set are shown in Figure 5. Figure 5A represents the binding mode of Compound 1 with a GOLD fitness score of 53.035. It has formed hydrogen bond interactions with D93, G95, T133, Q134, G291 and T293 and hydrophobic interactions with Y132, F169, and T292. The GOLD fitness score of RJC01726 was 68.289 and the mode of binding in the active site (Figure 5B) is similar to Compound 1. It has formed hydrogen bond interactions with T133, Q134 and T293 and hydrophobic interactions with D93, G95, F169 and T292. Asnx-2 has shown hydrogen bond interactions with T133, G291 and T293 as well as hydrophobic interactions with D93, Y132, F169 and T292 with a GOLD fitness score of 62.026 (Figure 5C). Figure 5D represents the overlay of most active Compound 1, RJC01726 and Asnx-2 at their binding modes. The pharmacophore overlay of the final hits compounds are shown in Figure 6 and their 2D representations are shown in Figure 7.

The hydrophobic interactions of the final hits compounds were observed using Ligplot program [39]. The novelty of the two hits compounds were confirmed using SciFinder search [40] and PubChem search [41].