Pharmacophore Modeling and Virtual Screening for the Discovery of New type 4 cAMP Phosphodiesterase (PDE4) Inhibitors

Type 4 cAMP phosphodiesterase (PDE4) inhibitors show a broad spectrum of anti-inflammatory effects in almost all kinds of inflamed cells, by an increase in cAMP levels which is a pivotal second messenger responsible for various biological processes. These inhibitors are now considered as the potential drugs for treatment of chronic inflammatory diseases. However, some recently marketed inhibitors e.g., roflumilast, have shown adverse effects such as nausea and emesis, thus restricting its use. In order to identify novel PDE4 inhibitors with improved therapeutic indexes, a highly correlating (r = 0.963930) pharmacophore model (Hypo1) was established on the basis of known PDE4 inhibitors. Validated Hypo1 was used in database screening to identify chemical with required pharmacophoric features. These compounds are further screened by using the rule of five, ADMET and molecular docking. Finally, twelve hits which showed good results with respect to following properties such as estimated activity, calculated drug-like properties and scores were proposed as potential leads to inhibit the PDE4 activity. Therefore, this study will not only assist in the development of new potent hits for PDE4 inhibitors, but also give a better understanding of their interaction with PDE4. On a wider scope, this will be helpful for the rational design of novel potent enzyme inhibitors.


Introduction
Type 4 cAMP-specific phosphodiesterase (PDE4) are a family of low k m 3',5'-cyclic adenosine monophosphate (cAMP)specific phosphodiesterases containing more than 20 isozymes encoded by four genes (PDE4A, PDE4B, PDE4C, and PDE4D) in mammals [1]. Even though four subfamilies share the conserved catalytic domain, each PDE4 gene plays a very important role in controlling the cell functions. PDE4s are taken as critical regulators of intracellular cAMP levels, cAMP signaling, and signal compartmentalization by their wide tissue distribution as well as differential expression and regulation among various cell types [1]. Thus many PDE4 inhibitors have showed remarked anti-inflammatory potential, by increasing cAMP levels. Recently the use of some newly marketed PDE4 inhibitors such as roflumilast, have been restricted due to their nausea and emesis. Therefore, the major pharmaceutical research focus in the field of chronic inflammatory diseases treatments, is to develop novel PDE4 inhibitors with high therapeutic index [1,2]. In our study, we successfully used pharmacophore modeling, database screening, and molecular docking approaches in identifying lead candidates to be used in potent PDE4 inhibitor design and thereby devising a new class of safer and effective anti-inflammatory agents.

Results and Discussion
Pharmacophore modeling A set of ten pharmacophore models was generated by a training set containing 28 compounds. Structures of the training set compounds are shown in Figure 1. The total cost values of ten pharmacophore models ranged from 106.849 to 120.562 ( Table 1). The cost difference between the null cost and total cost must be greater and it should be smaller between fixed cost and total cost values for a good pharmacophore model. In the present work, the first pharmacophore model (Hypo1) is basically composed of four features: two hydrogen bond acceptors, one hydrophobic region and one aromatic ring feature ( Figure 2). Hypo1 was developed with a fixed cost value of 99.761 and a null cost value of 204.947. Among the total cost values of ten pharmacophore models, Hypo1 scored the closest value to the fixed cost value than other models. The cost difference for the first pharmacophore model was 98.098. A cost difference value above 60 implies that the pharmacophore model correlates the estimated and experimental activity values more than 90% [9,10]. Therefore, Hypo1 could be considered as a good model. Based on the correlation coefficient, ten pharmacophore models were further evaluated. The correlation values of the generated pharmacophore models were greater than 0.91, and the values for the first three pharmacophore models were even higher, i.e., above 0.950. These results implied the capability of the pharmacophore model to predict the activity of the training set compounds. Hypo1 showed the highest correlation coefficient value of 0.963930, indicating its strong predictive ability. Moreover, RMSD values for ten pharmacophore models were less than 1, further supporting the predictive ability of these models. Among the ten pharmacophore models, Hypo1 was developed with better statistical values, such as higher correlation, larger cost difference and lower RMSD. Based on the experimental activity (IC 50 ) values, training set and test set compounds were categorized in following four groups: Highly active (IC 50 < 10nM, ++++), active (10≤ IC 50 < 200nM, +++), moderately active (200≤ IC 50 < 1000 nM, ++), and inactive (IC 50 ≥1000nM, +) [10]. Table 2 shows that activity values of all 28 compounds in the training set were predicted within their experimental activity scale, indicating the predictability of Hypo1. The pharmacophore mapping of most and least active compounds is shown in Figure 3. The most highly active compound (0.051 nM) mapped all the features of Hypo1, and the least active compound (4000 nM) missed hydrophobic and ring aromatic features. The reliability of Hypo1 has been further revealed.

Validation of the pharmacophore model
The validation process was performed by a test set of 28 compounds with diverse activity classes and different functional groups. Diverse conformers of these test set compounds were built in the same manner as for training set compounds using DS. Based on the geometric fit of these compounds over Hypo1, the estimated activity values were predicted for every test set compound. A correlation coefficient value of 0.948 is shown by the simple regression between the experimental and estimated activity values of training and test set compounds ( Figure 4). Two compounds out of 30 test set compounds were predicted in a different activity scale with a success rate of 93%. (Table 3). Especially, twenty-nine compounds from 30 test set compounds had error values less than 2.8, which was hardly different from the experimental and estimated activity values. The result showed a fairly good correlation between the experimental and estimated IC 50 values, indicating a good predictive capacity of Hypo1. Fischer's randomization method was additionally performed on the training set compounds to validate the statistical robustness of Hypo1. In this validation process, the experimental activities of the training set were scrambled randomly and the resulting training set was used in HypoGen module with the parameters chosen for the original pharmacophore generation [10]. To achieve a 95% confidence level, a set of 19 random spreadsheets is generated ( Figure 5). The result clearly indicated that none of the randomly generated pharmacophore models obtained from this validation method was produced with better statistical values than Hypo1. The result of Fischer's randomization test confirmed the statistical confidence of Hypo1. In order to verify whether the correlation between the experimental and estimated activities was mainly dependent on one particular compound in the training set, leave-one-out method was further used to perform the final validation. This was finished by recomputing the pharmacophore mode where one compound was excluded at a time. Under the same conditions which were used in the generation of the original pharmacophore model, 28 HypoGen calculations were carried out on 28 new training sets. The results indicated that compared to Hypo1, all the 28 new models generated by this method didn't have any meaningful difference. This result confirmed the confident level of Hypo1 that its correlation coefficient didn't depend on one particular compound in the training set.

Database screening and drug-likeness prediction
Based on the above validation results, the selected Hypo1 was used as a 3D query to search chemical databases including Specs (135556 molecules), Maybridge (59652 molecules) and NCI (238819 molecules), containing totally 434027 compounds. The inhibitory activity values of these compounds were estimated. A total of 220 compounds were firstly screened by restricting the minimum estimated activity to 1 nM. Then on the basis of Lipinski's rule of five and ADMET properties, these compounds were further screened to a number of 40. Finally, we subjected these 40 drug-like compounds along with the training set compounds to molecular docking study. Figure 6 lists the steps of the database screening procedure.

Docking study
To further refine the retrieved hits, forty drug-like hit compounds along with the training set compounds were docked into the active site of PDE4. The active site was defined based on the bound inhibitor in a crystal structure of PDE4 (PDB entry: 1XON). The docking poses were ranked by the binding free energy calculation. The binding free energy and molecular interactions with the active site residues were considered as important components in selecting the best hit compounds. The most active compound of training set (compound 1) with the binding free energy of -8.897 kcal/mol, has formed hydrogen bond interactions with Gln369 and His160 and ionic interactions with Zn 2+ and Mg 2+ ( Figure 7A). Moreover, this compound also exhibited a very important π-π interaction with the benzene ring of Phe372 and hydrophobic interactions with Met357, Thr333 and Gln369. On the basis of the molecular interaction of compound 1 and PDE4, we  selected twelve drug-like compounds with a binding free energy lower than -8.897 kcal/mol of compound 1 as final hits for further evaluation process. Intriguingly, these hits were obtained from Maybridge and Specs databases. Table 4 shows the list of twelve hits along with estimated activity values. The estimated activity value of every hit was less than 0.2 nM and the binding free energy was lower than -9.0 kcal/mol. Particularly, the first hit, PD00519, obtained from Maybridge database, had an estimated activity value of 0.091 and the lowest free energy of -11.671 kcal/mol. The molecular interaction between PDE4 and PD00519 is shown in Figure 7B. The carboxyl group of PD00519 formed ionic interactions with metal ions (Mg 2+ and Zn 2+ ) while the dimethyl group formed hydrophobic interactions with Met357, Thr333 and Gln369. PD00519 had also formed a vital π−π interaction with Phe372 and hydrogen bond interactions with Gln369 and His160 in the active site of PDE4. The molecular docking study indicated that compared with compound 1 in training set, the benzotrifluoride of PD00519 that was fitted adequately into the hydrophobic pocket of PDE4 has formed hydrophobic interactions with Met273, Leu229 and Ser208 in the active site. The understanding of this interaction between PDE4 and PD00519 will be beneficial to develop new design for novel PDE4 inhibitors. Figure 8A shows overlay of compound 1 and PD00519 in the active site of PDE4 while Figure 8B indicates their docking positions in the crystal structure of PDE4. Figure  9 represents a fairly good pharmacophore mapping of PD00519 and compound 1 on Hypo1 as predicted by molecular docking. The oxygen atoms of PD00519 and compound 1 that formed important interactions with metal ions in the pocket overlaid two HBA features of Hypo1. Moreover, the benzene rings of PD00519 and compound 1 having π−π interactions with Phe372 mapped RA features of Hypo1, while the hydrophobic groups mapped the HY features of Hypo1. The pharmacophore mapping of twelve hits on Hypo1 is depicted in Figure 10. Every hit compound has mapped four pharmacophoric features of Hypo1. Figure 11 shows the superimposition of twelve hit compounds on the Hypo1. The results indicated that hit compounds can produce perfect mapping with Hypo1. The benzene rings and oxygen atoms of hit compounds that overlaid the RA and HBA features of Hypo1, respectively, enabled considerable hydrophobic and polar interactions with the important amino acids in the active site. In addition, hydrophobic groups of these compounds that mapped the HY feature of Hypo1 interacted hydrophobically with the amino acids. Thus, in the design of potent inhibitors of PDE4, all twelve hit compounds which showed good results with respect to following properties, such as estimated activity, calculated drug-like properties and scores can be proposed as potential leads. Novelty search for compounds using SciFinder Scholar and PubChem search had also ascertained that these hits were not reported earlier for PDE4 inhibition. Therefore, we suggest that the identified compounds are novel and potent virtual leads for PDE4 inhibitor design.

Pharmacophore model generation
HypoGen module of Discovery Studio program (DS), version 2.5, from Accelrys (San Diego, USA) was used to perform all pharmacophore modeling calculations. While 13 compounds were selected from one article [3], other fourty-three structurally diverse compounds with inhibitory activity (IC 50 ) data were chosen as the training and test set compounds from similar articles [4][5][6][7][8] reported by same group of researchers. All the 56 compounds were tested for inhibition activity against PDE4 prepared from U937 cells (a cell line derived from human monocytes) by using the same experimental conditions. Based on the diversity of chemical structures and experimental activity values, we selected 28 compounds with wide activity range (0.051 to 4000 nM) as the training set. The two-dimensional (2D) chemical structures of all the compounds were built and subsequently converted to 3D structures in Discovery Studio program 2.5 (DS). For the compounds in the the data set, Pharmacophore Modeling of PDE4 Inhibitors PLOS ONE | www.plosone.org CHARMM forcefield was used to perform energy minimization process. Poling algorithm was used to generate a maximum of 230 diverse conformations with the energy threshold of 15 kcal mol -1 above the calculated energy minimum for every compound in the dataset [9,10]. Diverse Conformer Generation protocol running with Best/Flexible conformer generation option was applied to generate multiple conformers. By performing a more rigorous energy minimization in both torsional and cartesian space, this method ensures the best coverage of conformational space [11]. All the 28 compounds of training set were submitted to the HypoGen module of DS. The minimum and maximum count for all the features in the hypothesis run were set of 1 and 6, respectively. Uncertainty value was set to 2 and the minimum inter-feature distance was set to 2.5 Å from the default value of 2.97 Å. Hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic (HY) and ring aromatic (RA) features were used to generate ten pharmacophore models using 3D QSAR pharmacophore generation of DS [9,10]. All other parameters used in HypoGen module were kept at their default settings [9,10]. In this study, the top 10 hypotheses returned by the hypotheses generation process were selected for further calculations.

Pharmacophore model evaluation
Based on cost functions and other statistical parameters which were calculated by HypoRefine module during hypothesis generation, the quality of the generated pharmacophore models was evaluated. The best pharmacophore model should have a high correlation coefficient, low RMSD values and total cost that should be away from the null cost and close to the fixed cost [9,10]. All of these cost values are reported, and a difference of 40-60 bits between the total and null costs suggests a 75-90% chance of representing a true correlation in the data [9,10]. To investigate the ability to estimate the activity of new compounds, the selected pharmacophore model was further validated by three methods including test set method, Fischer's randomization test and leave-one-out method. 28 diverse compounds were used as the test set to validate the pharmacophore model. For Fischer's randomization test, 95% confidence level was chosen in this validation study and the 19 random spreadsheets were constructed. Finally, we performed the leave-one-out methodology for the cross validation of the model by using the same parameters as used for generating original pharmacophore model, thus 28 pharmacophore models were generated. To ensure the influence of each compound from the training set in the generation of selected pharmacophore model, one compound at a time from 28 compounds was left [10][11][12][13].

Virtual screening
In order to identify novel hit compounds, the best pharmacophore model after validation was used as 3D structural search query to screen three chemical databases including Specs (135556 molecules), Maybridge (59652 molecules) and NCI (238819 molecules), respectively. Search 3D Database protocol with Best/Flexible search option was applied in database screening. The hits identified through database screening were further filtered using estimated activity, Lipinski's rule of five [14], and ADMET properties [15][16][17][18]. A Lipinski-positive compound has (i) a molecular weight < 500; (ii) < 5 hydrogen bond donor groups; (iii) < 10 hydrogen bond acceptor groups and (iv) an octanol/water partition coefficient (Log P) value < 5 [9,10].

Molecular docking
The docking study was performed using the Molecular Operating Environment (MOE) software (Chemical Computing Group Inc.). The crystal structure of PDE4 obtained at a resolution of 1.72 Å was downloaded from the protein data bank (PDB entry: 1XON). This structure was then protonated in the Molecular Operating Environment (MOE) via MMFF94x force field. The hits identified through database screening were subjected to molecular docking studies. The active site was defined with a 6 Å radius around the bound inhibitor. The triangle matcher algorithm of the MOE software packages was used to dock the identified hits into the protein active site. According to this algorithm, different poses were generated by aligning ligand triplets of atoms on triplets of alpha spheres. For all scoring functions, lower scores indicated more favorable poses. The scoring function of these compounds has to obey the following parameters: (1) Specifying ASE Scoring to use for ranking the poses output by the placement stage; (2) Specifying Forcefield Refinement to use to relax the poses, respectively; (3) Specifying Affinity dG Scoring to use for ranking the poses output by the refinement stage. The free energy of binding was calculated from the contributions of hydrogen bond, ionic, hydrophobic and van der Waals interactions between the protein and ligand, intramolecular hydrogen bonds and strains of the ligand. We observed in the S field that the docking poses were ranked by the binding free energy calculation.

Conclusions
In the present work, a highly correlating (r = 0.963930) pharmacophore model (Hypo1) containing two hydrogen bond acceptors, one hydrophobic region and one aromatic ring feature, was selected through various parameters such as total cost, correlation coefficient and cost difference. Further validation was done by using test set prediction, Fischer randomization method and leave-one-out method. Result of these validation tests showed that Hypo1 could accurately predict the active compounds, it has better statistical values compared to other randomly generated pharmacophore models and its correlation coefficient is not solely depended on a single compound. This validated Hypo1 was led to database screening for identifying compounds which can be used as potent PDE4 inhibitor design. Further studying these compounds by drug-like filtrations and molecular docking, also suggested the robustness of Hypo1. In the end, twelve structurally diverse compounds having high estimated activity and strong molecular interactions with key active site amino acids of PDE4 were identified. Therefore, the results of this study will assist, not only in the development of new potent hits for PDE4, but also in providing a better understanding of the interaction between PDE4 and inhibitors. This will in turn be beneficial to the rational design of novel potent enzyme inhibitors.

Supporting Information
Text S1. Pharmacophore modelling and 3D database search.

(DOC)
Text S2. The steps of the procedure in our study.
(DOC) Figure S1. The Chemical structures of twelve hit compounds from databases. (TIF)