35 indole derivative-based PDE-4 inhibitors [15Chakraborti AK, Gopalakrishnan B, Sobhia ME, Malde A. 3D-QSAR studies of indole derivatives as phosphodiesterase IV inhibitors Eur J Med Chem 2003; 38(11-12): 975-82.] were used in this study. The molecular structures and their biological activities against PDE-4 are given in Table 1. The compounds contain an indole moiety which replaces the “rolipram-like” 3-methoxy-4-cyclopentoxy motif. The activities were determined from measurement of serum TNF-α levels in LPS challenged mice (mouse endotoxemia model) [24Hulme C, Moriarty K, Miller B, et al. The synthesis and biological evaluation of a novel series of indole PDE4 inhibitors I Bioorg Med Chem Lett 1998; 8(14): 1867-72.].

This dataset covers a variety of molecules with various levels of inhibitory activities against PDE-4. Fig. (1) shows the activity distribution of the 35 indole derivatives. To develop and rigorously validate a QSAR model, 28 molecules were selected as the training set, and the remaining 7 molecules were held out as the test set. The specific selection of the training and test sets was determined by either randomly splitting or clustering-based splitting of the original dataset.

The 3D conformations of all the inhibitors were generated using the OMEGA 2.0 program available from OpenEye Scientific (OMEGA 2.0, OE Scientific, NM, USA). In this study, only the best conformer of each molecule was used to generate the pharmacophore points.

Fig. (2) depicts the overall work flow of the structure-based QSAR method employed in this study. The novel aspect of this method is the way by which it generates the molecular descriptors by the SB-PPK (Structure-based Pharmacophore Key) program. Unlike conventional 2D or 3D molecular descriptor generation methods where ONLY the structural information of the inhibitor molecules is used to characterize the molecules, the SB-PPK program uses the structural information of the target binding site to measure the inhibitor molecules. Thus, target-specific molecular descriptors are generated, which reflects not only the inhibitor molecules but also the way these molecules interact with the binding site of the target. The general process of SB-PPK involves: 1) extraction of pharmacophore features from the target binding site; 2) perception of pharmacophore features of small molecules; and 3) generation of the structure-based pharmacophore key (SB-PPK) descriptors. Specifically, it derives the pharmacophoric centers from the 3D structure of a protein-ligand (R-L) complex, and uses them to generate the pharmacophore reference. It then generates structure-based pharmacophore key (SB-PPK) descriptors for the inhibitor molecules based on their pharmacophoric matching patterns with the above reference. Once molecular descriptors are generated, PLS (Partial Least Square) method is employed to establish the relationship between these molecular descriptors and the activity of the inhibitors. Details of the individual steps are described as follows.

Fig. (2) The overall flowchart of the SB-PPK QSAR Method.

The LigandScout program [25Wolber G, Langer T. LigandScout: 3-D pharmacophores derived from protein-bound ligands and their use as virtual screening filters J Chem Inf Model 2005; 45(1): 160-9.] was first employed to derive structure-based pharmacophore centers (Fig. 3a). The advanced 3D pharmacophore model option of the LigandScout program was used, which allows multiple pharmacophore feature types on the same heavy atom. The pharmacophore feature types used in this study include hydrogen bond donor (D), hydrogen bond acceptor (A), positively charged (P), negatively charged (N), and lipophilic (L) centers. Pharmacophore feature pairs are then generated by making all possible pair-wise combinations of the above pharmacophore centers. A specific pharmacophore feature pair is determined by both the feature types involved and the distance between the two centers. For example, AL in Fig. (3a) indicates a hydrogen bond acceptor (A) and lipophilic (L) feature pair; and the distance between them is 4.5 A. Another example, AP, means a hydrogen bond acceptor (A) and a positive (P) feature pair and their inter-feature distance is 3.8 A. All possible feature pairs are generated, and this recorded information reflects the pharmacophoric characteristics of the target binding site, and will be used as the reference for generating the SB-PPK descriptors of the inhibitor molecules. To avoid confusion with traditional pharmacophore fingerprints, we note that we do not need “distance binning” to capture different possible distance ranges because the pharmacophore features pairs are derived from the 3D structure of the target protein. Thus, all the feature pairs have defined distances.

Fig. (3) a. Deriving the pharmacophore feature pairs (the reference) from the receptor-ligand complex. b. Deriving the pharmacophore feature pairs for inhibitor molecules. c. Generating structure-based pharmacophore key (SB-PPK) descriptors.

The LigandScout program[25Wolber G, Langer T. LigandScout: 3-D pharmacophores derived from protein-bound ligands and their use as virtual screening filters J Chem Inf Model 2005; 45(1): 160-9.] is used to perceive the pharmacophoric groups on small organic molecules (Fig. 3b). All pair-wise combinations of the perceived pharmacophore centers are checked for the feature types and interfeature distances. In the example given in Fig. (3b), there is an AL pair and an AP pair, and the inter-feature distances for each pair are 4.6 A and 3.8 A, respectively.

Once the pharmacophore feature pairs for both the receptor and ligand molecules are generated, we can compare them to determine pattern matches. For example, Fig. (3c) shows that AL’s and AP’s from the receptor’s pharmacophore feature pairs (the reference) and small molecule’s pharmacophore pairs match with each other. Thus, the corresponding positions for AL and AP are increased by 1 count. Both fingerprint-like binary descriptors and feature pair counts descriptors can be generated for a given ligand.

The relationship between the molecular descriptors and the biological activities is modeled by the PLS (Partial Least Square) method [26Hokuldsson A. A PLS regression methods J Chemometr 1988; 2: 211-8.]. Specifically, the PLS implementation in the MOE package (Chemical Computing Group, Montreal, Canada) was used in this study. Other types of machine learning methods can be used as well to analyze the dataset.

Model validation is a critical step in any QSAR study. As shown in Fig. (4), we have adopted a standard workflow, similar to what Tropsha group advocate [23Golbraikh A, Shen M, Xiao Z, Xiao YD, Lee KH, Tropsha A. Rational selection of training and test sets for the development of validated QSAR models J Comput-Aided Mol Des 2003; 17(2-4): 241-53.], to build robustly predictive QSAR models. Rational splitting of the dataset into training and test sets as well as obtaining the statistics for both the training set and test sets is critical in this model validation protocol. Thus, a given dataset is split into multiple pairs of training and test sets. Both random shuffling and clustering-based [27Labute P. Quasar-Cluster: A different view of molecular clustering Chem Comp Group J 2002 [Jan 5, 2008]; http://www.chemcomp.com/journal/cluster.htm] splitting of the data set can be used for dataset partitioning. In the random shuffling approach, the first 7 molecules were held out as the test set after the order of molecules in the data set was randomly shuffled. In the clustering-based partitioning approach, the dataset was split in such a way that both the training set (28 molecules) and the test set (7 molecules) had the most coverage of the clusters. For each split, PLS was used to analyze the training set in order to build a PLS regression model. Regression coefficients were then calculated for both the training set (r²) and the test set (R²). While the former reflects the model quality as measured on the training set, the latter (R²) is a more robust indicator of the predictive power of the model. In the end, only those model(s) with both coefficients greater than set thresholds are retained as the final model(s).

Fig. (4) Workflow for rigorous validation of QSAR models.

As described before (Cf. Method), LigandScout was used to analyze the interactions between PDE-4 inhibitors and the binding site of PDE-4, based on the X-ray structure of the inhibitor-bound PDE-4 (1xon.pdb). This analysis has revealed nine pharmacophore centers. The 3D rendering and the 2D depiction of the interactions are shown in Fig. (5a) and Fig. (5b), respectively. The following occupies the four hydrophobic centers: the cyclopentane ring, two chlorine atoms on the pyridine ring, and the central benzene ring. The cyclopentane ring interacts with hydrophobic residues in the PDE-4 binding pocket: MET337, MET357, PHE340, PHE372, and ILE336. Chlorines on the pyridine ring interact with hydrophobic residues MET273, PHE340, ILE336, and LEU319, and the central benzene ring interacts with ILE336, PHE372 and TYR159. The amide nitrogen occupies a hydrogen bond donor center, which interacts with a water molecule that serves as a bridge to the enzyme’s binding pocket. Four hydrogen bond acceptor centers are identified and they are occupied by the following atoms/groups: the methoxy oxygen on the benzene ring, amide oxygen, pyridine nitrogen, and the ether oxygen (connecting the benzene and the cyclopentane rings). The methoxy oxygen and the ether oxygen atoms both interact with GLN369. The amide oxygen and the pyridine nitrogen interact with the small molecule EDO744 and a water molecule, respectively.

Fig. (5) a. 3D rendering of the pharmacophore features derived from Ligand-bound PDE4 structure using the LigandScout program. b. 2D depiction of the pharmacophore features perceived from the Ligand-bound PDE4 structure using the LigandScout program.

Based on the above-generated pharmacophore features, we created the pharmacophore reference in terms of pharmacophore feature pairs (Cf. Method). Each structure based descriptor is determined by the feature types involved in a pharmacophore feature pair and the distance between the two feature centers. Table 2 shows the pharmacophore feature pairs (the reference) derived from the PDE4-inhibitor complex (1xon.pdb). It has six AA pairs with inter-feature distances ranging from 6.34 to 10.69 A. Four AD pairs with inter-feature distances ranging from 2.29 to 6.50 A. Likewise, there are sixteen AL pairs, four DL pairs and six LL pairs with different inter-feature distances. Together, this set of feature pairs uniquely characterizes the binding pocket of the PDE-4 enzyme, and serves as a reference with which candidate inhibitor molecules can be compared to generate their molecular descriptors.

The descriptors of the PDE-4 inhibitors were determined based on the matches between the feature pairs in the reference and those in each individual molecule. Fig. (6a) shows the values of calculated SB-PPK descriptors of PDE-4 inhibitor 3 (mol.3 in Table 1). The average descriptor values over the dataset of 35 PDE-4 inhibitors are presented in Fig. (6b). In inhibitor 3, the value of AL10 is 4, which indicates 4 possible AL10 (inter-feature distance 3.98) matches between the feature pairs of inhibitor 3 and those of the binding site reference. Obviously, different molecules have different number of matches between different feature pairs. In Fig. (6b), the error bars show the variations of each descriptor value around the averages, amongst the different inhibitors in the dataset under study. Those descriptors with larger variations tend to have more important roles in QSAR models.

Fig. (6) a. The SB-PPK descriptor values of a PDE4 inhibitor (mol.3 of Table 1). b. The SB-PPK descriptor values averaged over the 35 PDE4 Inhibitors.

We follow the workflows outlined in Fig. (2) and Fig. (4) to conduct the QSAR model building and model validation. A variety of factors that may have significant effect on QSAR model quality have been thoroughly examined and discussed below. The quality of the final model is compared with previous work.

To obtain predictive QSAR models, we systematically performed PLS analysis with different number of principle components used in the regression model. Fig. (7) shows how the number of principle components affects the quality of the resulting models for the best partitioning of training and test sets. The conventional correlation coefficient r² obtained for the training set monotonically increases while the predictive R² obtained for the testing set reaches a plateau before it starts to decrease. When the number of principle components (PC) is less than 3, the resulting QSAR models are not predictive for either the training set or the test set. When the number of principle components is over 8, the predictive ability for the test set starts to go down. This indicates potential over-training at higher number of principle components. According to our model validation protocol (Fig. 4), the models with the number of principle components between 5 and 8 are predictive for both the training and the test sets. The best models should be a trade-off between conventional r² and the predictive R². Also, to compare with other literature work which uses number of principle components of 6 [15Chakraborti AK, Gopalakrishnan B, Sobhia ME, Malde A. 3D-QSAR studies of indole derivatives as phosphodiesterase IV inhibitors Eur J Med Chem 2003; 38(11-12): 975-82.], we decided to use 6 components to build the final best models.

Fig. (7) Quality of the QSAR models depends on the number of principle components (PC) employed.

To build high quality, predictive QSAR models, we systematically performed dataset splitting into training and test sets using both random shuffling and clustering-based splitting methods. For each partitioning approach, we obtained ten different training/test sets, and built 90 different QSAR models with the number of principal components ranging from 2 to 10. Of all 180 QSAR models, the conventional correlation coefficient r² ranges from 0.69 to 0.89, and the predictive R² values range from 0.05 to 0.75. The splitting which renders the highest r² and R² for corresponding models with number of principle components between 5 and 8 is selected as the partitioning of training and test sets to build the final QSAR models.

Fig. (8a) and Fig. (8b) show the predictive quality of the best QSAR models for PDE-4 inhibitors built with the SB-PPK descriptors and PLS method. Table 3 present the comparison between our best models with those of other methods. We have compared the results obtained with CoMFA [15Chakraborti AK, Gopalakrishnan B, Sobhia ME, Malde A. 3D-QSAR studies of indole derivatives as phosphodiesterase IV inhibitors Eur J Med Chem 2003; 38(11-12): 975-82.], CoMSIA [15Chakraborti AK, Gopalakrishnan B, Sobhia ME, Malde A. 3D-QSAR studies of indole derivatives as phosphodiesterase IV inhibitors Eur J Med Chem 2003; 38(11-12): 975-82.], 2D QSAR (using MOE package) and SB-PPK (our own method). All methods use PLS as the modeling tool with number of principle components set at 6. Both CoMFA and CoMSIA afforded highly predictive models based on the training set; however, the predictive R² for the test set are both significantly lower than that for the training set. Our method, on the other hand, afforded much more balanced predictive quality, for both the training set and the test set. Thus, based on the reported correlation coefficients (for both training and test sets), we believe our method (SB-PPK plus the specific workflow process) afforded the best model both in terms of predictive ability for the training set and, more importantly, in terms of the predictive power on the test set.

Fig. (8) a. QSAR model-predicted activity versus the actual activity for the molecules in the training set. b. QSAR model-predicted activity versus the actual activity for the molecules in the test set.

One of the important tasks of QSAR analysis is to identify critical features that are responsible for explaining the SAR trend in the dataset. In order to do so, we define a score for each SB-PPK descriptor based on the magnitude of the coefficient for that descriptor in the PLS regression models. Specifically, we define the score according to equations 1 and 2.

Score (i) is the score of a given descriptor contributed by a particular model, i, and rank is the rank order of a particular descriptor according to the magnitude of its coefficient in that PLS model. The summation is conducted over all the predictive models obtained. In this study, we base the score on four predictive PLS models, corresponding to the number of principal components 5, 6, 7, and 8. Fig. (9a) shows that descriptor LL1, AL12, AA6, AA1 and AL2 are relatively important in the QSAR models. These descriptors correspond to the pharmacophore feature pairs depicted in Fig. (9b). Thus, the SB-PPK descriptor based QSAR analysis can reveal critical features that correspond to the pharmacophore features of the target binding pocket. This descriptive ability is very useful and often desired by the medicinal chemist working on compound optimizations.

Fig. (9) a. Descriptor significance scores for the most important SB-PPK descriptors. b. Mapping the significant pharmacophore feature pairs according to the descriptor significance scores.