A new structure-based QSAR method affords both descriptive and predictive models for phosphodiesterase-4 inhibitors.

We describe the application of a new QSAR (quantitative structure-activity relationship) formalism to the analysis and modeling of PDE-4 inhibitors. This new method takes advantage of the X-ray structural information of the PDE-4 enzyme to characterize the small molecule inhibitors. It calculates molecular descriptors based on the matching of their pharmacophore feature pairs with those (the reference) of the target binding pocket. Since the reference is derived from the X-ray crystal structures of the target under study, these descriptors are target-specific and easy to interpret. We have analyzed 35 indole derivative-based PDE-4 inhibitors where Partial Least Square (PLS) analysis has been employed to obtain the predictive models. Compared to traditional QSAR methods such as CoMFA and CoMSIA, our models are more robust and predictive measured by statistics for both the training and test sets of molecules. Our method can also identify critical pharmacophore features that are responsible for the inhibitory potency of the small molecules. Thus, this structure-based QSAR method affords both descriptive and predictive models for phosphodiesterase-4 inhibitors. The success of this study has also laid a solid foundation for systematic QSAR modeling of the PDE family of enzymes, which will ultimately contribute to chemical genomics research and drug discovery targeting the PDE enzymes.

INTRODUCTION

The cyclic nucleotide phosphodiesterases (or PDE’s) comprise a group of enzymes that degrade the phosphodiester bond in the second messenger molecules cAMP and cGMP. They regulate the localization, duration, and amplitude of cyclic nucleotide signaling within subcellular domains. Thus, PDE’s are important regulators of signal transduction mediated by these second messenger molecules [1-3].

According to their sequence homology, substrate specificity and pharmacological properties, the PDE enzymes are classified into 11 families (i.e., PDE1 - PDE11) [2]. Even though different subtypes of this super-family are functionally related, they do have different substrate specificities. For example, PDE-4, PDE-7 and PDE-8 are cAMP selective, whereas PDE-5, PDE-6 and PDE-9 are cGMP selective. Other members of this family (e.g., PDE-1, PDE-2, PDE-3, PDE-10 and PDE-11) can hydrolyse both cAMP and cGMP.

Phosphodiesterases have proven to be a therapeutically important class of drug targets. Sildenafil is an inhibitor of cGMP-specific phosphodiesterase (PDE-5) and used to treat erectile dysfunction (ED) [4]. Inhibitors of PDE-4 have been explored as potential drugs for chronic obstructive pulmonary disorder (COPD) [5]. Inhibitors of other PDE enzymes are being explored for coronary heart disease [6], dementia, depression, and schizophrenia [7]. Thus, there are continuing interests in discovering novel inhibitors of PDE enzymes that are potent and selective for specific subtypes.

We have been interested in modeling the inhibitors of PDE-4 due to their potential in treating neurodegenerative diseases and cognitive disorders. PDE-4 selectively hydrolyzes cAMP and has been mostly targeted for the treatment of inflammation, and in particular chronic obstructive pulmonary disorder (COPD) [5]. However, selective phosphodiesterase-4 (PDE-4) inhibitors have been demonstrated recently as novel agents to treat memory deficit [8] and neurodegeneration [9]. Preclinical studies indicate that these PDE-4 inhibitors can counteract deficits in long-term memory caused by pharmacological agents, aging or over-expression of mutant forms of human amyloid precursor proteins [9]. Although a number of major advances have been made with respect to PDE-4 inhibition, most of the PDE-4 inhibitors have unacceptable side effects, particularly nausea and emesis [10, 11]. Thus, it is critically important to understand the structural determinants of potent and selective PDE-4 inhibitors, so that we can rationally design new molecules that minimize the undesirable side effects.

Most PDE-4 inhibitors can be classified into three structural families: rolipram-related compounds, xanthine derivatives and nitraquazone analogues [12]. This structural diversity suggests that the binding site of the PDE-4 enzyme presents a number of alternative pharmacophores that are capable of interacting with different classes of small molecule inhibitors. Many investigators have tried to understand the structural basis for PDE-4 inhibition (both potency and selectivity) by constructing 3D QSAR (quantitative structure-activity relationship) models based on known PDE-4 inhibitors [12-16]. Even though these models provided preliminary pharmacophoric insights into PDE-4 inhibitors, the absence of the 3D structure of the target enzyme (i.e., PDE-4) in these models poses a major limitation on the reliability and interpretability of these models. Other traditional 2D QSAR analyses often result in only statistical models with little to offer in terms of structural insights that the medicinal chemist can use to optimize the inhibitors [12, 17]. Thus, it is highly desirable to develop structure-based, both descriptive and predictive models.

Here, we describe the application of a new QSAR method to the analysis and modeling of PDE-4 inhibitors. This method takes advantage of the X-ray structural information of the target enzyme (i.e., PDE4) to characterize its small molecule inhibitors. Other groups have published different methods to conduct conceptually similar QSAR analysis on other targets, notably the work by Hopfinger’s group and Cruciani’s group [18-22]. Our method derives the pharmacophoric centers from the 3D structure of PDE-4’s binding pocket, and uses them as pharmacophore reference. It then generates structure-based pharmacophore key (SB-PPK) descriptors for the inhibitors based on their pharmacophoric matching patterns with the above reference. We have analyzed 35 indole derivative-based PDE-4 inhibitors [15]. Partial Least Square (PLS) analysis has been employed to obtain the predictive models measured by both the training set and test set statistics [23]. We will compare our results with previous studies of the same data set by more conventional methods and highlight some novel aspects of our approach.

METHODS AND MATERIALS

Dataset

35 indole derivative-based PDE-4 inhibitors [15] 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) [24].

Table 1.

Molecular Structures and Inhibition Activities (-logM) of the Data Set

This dataset covers a variety of molecules with various levels of inhibitory activities against PDE-4. Fig. () 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.

The Overall Work Flow of the QSAR Method

Fig. () 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.

Deriving Pharmacophore Feature Pairs from Ligand-Receptor (R-L) Complex

The LigandScout program [25] 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.

Perception of Pharmacophore Feature Pairs for a Ligand

The LigandScout program[25] 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.

Generation of Structure-Based Pharmacophore Key (SB-PPK) Descriptors

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.

Statistical Analysis with Partial Least Square (PLS) Method

The relationship between the molecular descriptors and the biological activities is modeled by the PLS (Partial Least Square) method [26]. 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.

Validation of the QSAR Models

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 [23], 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 [27] 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 (r2) and the test set (R2). While the former reflects the model quality as measured on the training set, the latter (R2) 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.

RESULTS AND DISCUSSION

Pharmacophore Features of the PDE-4 Binding Site

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. () and Fig. (), 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.

Pharmacophore Reference Generated from the PDE-4 Binding Site

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.

Table 2.

Structure-BASED Pharmacophore Feature Pairs of PDE4 Binding Site

AA1a	AA2	AA3	AA4	AA5	AA6	AD1	AD2	AD3
6.34b	10.69	2.72	4.94	5.05	9.91	6.50	2.29	4.24
AD4	AL1	AL10	AL11	AL12	AL13	AL14	AL15	AL16
6.15	4.79	3.98	3.99	7.84	2.31	7.03	8.47	2.93
AL2	AL3	AL4	AL5	AL6	AL7	AL8	AL9	DL1
7.98	8.56	2.86	5.91	3.12	4.48	3.61	10.69	7.47
DL2	DL3	DL4	LL1	LL2	LL3	LL4	LL5	LL6
2.98	3.09	3.68	7.51	9.98	4.88	5.45	5.36	5.98

The code represents pair wise combinations of pharmacophore features in the PDE4 structure. For example, there are six AA pairs (AA1 to AA6) and four AD pairs (AD1 to AD4)where “A” and “D” represent hydrogen bond acceptor and donor features, respectively.

The numbers are the inter-feature distances between the corresponding feature pairs. For example, the distance between “A” and “A” in AA1 is 6.34 (A).

Structure-based Pharmacophore Key (SB-PPK) Descriptors of PDE-4 Inhibitors

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.

QSAR Model Building with Partial Least Square (PLS) Method

We follow the workflows outlined in Fig. () and Fig. () 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.

Effect of Number of Principle Components on the Quality of QSAR Models

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 r2 obtained for the training set monotonically increases while the predictive R2 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 r2 and the predictive R2. Also, to compare with other literature work which uses number of principle components of 6 [15], 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.

Effect of Dataset Partitioning on QSAR Model Quality

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 r2 ranges from 0.69 to 0.89, and the predictive R2 values range from 0.05 to 0.75. The splitting which renders the highest r2 and R2 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.

The Best SB-PPK Models have Better Predictive Ability than those Obtained with other Methods

Fig. () and Fig. () 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 [15], CoMSIA [15], 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 R2 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.

Table 3.

Comparison of Model Performance in Terms of r2 and Predicted R2 when Different Set of Descriptors are Used

Molecular Descriptors	Regression Method	Number of Components	Correlation Coefficient (r2)	Predicted R2
SBPPK	PLS d	6	0.747	0.624
MOE-2D a	PLS d	6	0.664	0.579
CoMFA b	PLS	6	0.986	0.560
CoMSIA c	PLS	6	0.967	0.590

2D descriptors generated using the MOE software (Chemical Computing Group, Toronto, CA).

CoMFA analysis of the PDE4 inhibitors is cited from [15].

CoMSIA analysis of the PDE4 inhibitors is cited from [15].

Partial Least Square (PLS) program implemented in the MOE package.

Critical Pharmacophore Features have been Revealed by the SB-PPK QSAR Analysis

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. () 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. (). 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.

CONCLUSIONS

We have successfully applied a new structure-based QSAR formalism to study the SAR (Structure-activity Relationship) of 35 PDE-4 inhibitors. This new QSAR formalism describes molecular structures based on the matching of their pharmacophore feature pairs with those (the reference) of the target binding pocket. This reference is derived from the X-ray crystal structures of the target under study. Thus, these descriptors are target-specific and easy to interpret. The PLS (partial least square) method is then employed to reveal the correlation between these target-specific descriptors and the biological activity, and help identify the critical pharmacophore features for PDE-4 inhibition.

To ensure the predictive ability of the built QSAR models, we have adopted a rigorous model validation protocol that emphasizes rational design of the training and test sets, and uses the statistics calculated for both the training set (r2) and test set (R2) as the model selection criteria. Consequently, our models display highly balanced predictive ability for both the training and the test sets (Table 3). In contrast, CoMFA and CoMSIA models have shown signs of over fitting, with high r2 for the training set but much lower R2 for the test set. Thus, this model validation protocol based on both training and test set statistics is critical for the success of this study.

The new structure-based descriptors overcome the drawbacks of traditional descriptors that ignore the binding pocket information. These target-specific descriptors could offer structural insights into the critical features responsible for the potency of the inhibitors. The feature importance analysis proposed in this study has helped fulfill this task. Thus, our models afford both predictive and descriptive models for PDE-4 inhibitors.

Since the structures of many other PDE enzymes have also been solved (e.g., PDE5: 2h44, 2h42, 2h40; PDE-9: 2hd1; PDE3B: 1soj, 1so2; PDE2A: 1z1l and PDE1B: 1taz), the success of this study has also laid a solid foundation for studying the inhibition mechanisms of other PDE enzymes using this new structure-based QSAR formalism. Such a systematic study of the whole family of PDE enzymes will eventually afford insights into the structural requirements for both the potency and selectivity of PDE inhibition, and ultimately contributes to chemical genomics research and therapeutic development targeting the PDE family of enzymes.