Structural insight into selective phosphodiesterase 4B inhibitors: pharmacophore-based virtual screening, docking, and molecular dynamics simulations

Phosphodiesterase (PDE4) is the predominantly expressing family of PDE enzyme in immune and inflammatory cells. Inhibition of PDE4 has been reported to suppress a diverse spectrum of inflammatory responses both in vitro and in vivo. Many PDE4 inhibitors in development are shown to be efficacious in animal models of various inflammatory disorders, such as asthma, COPD, psoriasis, inflammatory bowel diseases, and rheumatoid arthritis, as well as in clinical trials for asthma and COPD. Therefore, PDE4 might prove as potential therapeutic targets in inflammatory and autoimmune diseases. The PDE4 enzyme family consists of four members (named PDE4A, PDE4B, PDE4C, and PDE4D) that are highly expressed in neutrophils, monocytes, central nervous system, and smooth muscles of the lung. Non-selective PDE4 inhibitors like Rolipram are associated with severe side effects like nausea and emesis. Second-generation PDE4 inhibitors like Roflumilast are reported with improved side effects, but narrow therapeutic window. The behavioral correlation of emesis in mice by deleting PDE4D encoding gene confirmed that inhibition of PDE4D, but not PDE4B is responsible for emetic effect of non-selective PDE inhibitors (Jin, Ding, & Lin, 2012). Therefore, selective PDE4B inhibitors can provide anti-inflammatory efficacy without any side effect. However, conserved active site residues of PDE4B and PDE4D make it difficult to design selective PDE4B inhibitors. In the line of identifying selective PDE4B inhibitors, a combined receptorand ligand-based approach can improve the binding affinity and selectivity of PDE4B inhibition. Thirty-four PDE4B structures are available in RCSB-PDB. Hitherto, no structural classification of available co-crystallized PDE4B inhibitors is reported. An in-depth analysis of available PDE4B-inhibitor complexes may help to find empirical parameters for selective enzyme inhibition. Chen et al. (2010) developed multiple pharmacophore models based on the 18 crystal structures of phosphodiesterase 4 (PDE4) with co-crystal ligand for discovering potential PDE4 inhibitors. Srivani, Usharani, Jemmis, and Sastry (2008) performed the crystal structure analysis of PDE4B and PDE4D which led to find significant variations in the M-loop region, which is the integral part of the active site of PDE4B and PDE4D. Ke (2004) performed the similar crystal structure analysis for PDE4 and PDE5 with respect to the metal binding site. In present study, first we retrieved the information from available PDE4B crystal structures and grouped the crystal ligands in different classes based on chemical scaffolds. We have employed structure and ligand-based dual approach for identifying potential and selective molecules for PDE4B. Structure-based approach presents advantage over ligand-based approach by incorporating spatial distribution of hotspot. At the same time ligandbased approach helps to identify key pharmacophoric features in ligands essential for bioactivity. We combined the chemical features of reported crystal ligands of PDE4B in RCSB PDB data bank and developed ligandbased pharmacophores. Using these queries, SPECS database of commercially available small molecules was screened to find new PDE4B inhibitors. To predict the selectivity in PDE4B over PDE4D, the screened molecules were docked in PDE4B (PDB ID: 3G45) and PDE4D (PDB ID: 3G4G) crystal structures using Glide. The stability of top-scoring ligand in PDE4B active site was confirmed by a 10 ns MD simulations (Figure 1). Our results provided critical insight for designing selective PDE4B inhibitors. This work is organized as follows. The crystal structure analysis, pharmacophore modeling, virtual screening (VS), docking and ADME studies, MD simulation materials, and reproducible methods are provided in next section. Further part includes is the analysis and discussion of crystal structure analysis and pharmacophore

structures of phosphodiesterase 4 (PDE4) with co-crystal ligand for discovering potential PDE4 inhibitors. Srivani, Usharani, Jemmis, and Sastry (2008) performed the crystal structure analysis of PDE4B and PDE4D which led to find significant variations in the M-loop region, which is the integral part of the active site of PDE4B and PDE4D. Ke (2004) performed the similar crystal structure analysis for PDE4 and PDE5 with respect to the metal binding site.
In present study, first we retrieved the information from available PDE4B crystal structures and grouped the crystal ligands in different classes based on chemical scaffolds. We have employed structure and ligand-based dual approach for identifying potential and selective molecules for PDE4B. Structure-based approach presents advantage over ligand-based approach by incorporating spatial distribution of hotspot. At the same time ligandbased approach helps to identify key pharmacophoric features in ligands essential for bioactivity. We combined the chemical features of reported crystal ligands of PDE4B in RCSB PDB data bank and developed ligandbased pharmacophores. Using these queries, SPECS database of commercially available small molecules was screened to find new PDE4B inhibitors. To predict the selectivity in PDE4B over PDE4D, the screened molecules were docked in PDE4B (PDB ID: 3G45) and PDE4D (PDB ID: 3G4G) crystal structures using Glide. The stability of top-scoring ligand in PDE4B active site was confirmed by a 10 ns MD simulations (Figure 1). Our results provided critical insight for designing selective PDE4B inhibitors.
This work is organized as follows. The crystal structure analysis, pharmacophore modeling, virtual screening (VS), docking and ADME studies, MD simulation materials, and reproducible methods are provided in next section. Further part includes is the analysis and discussion of crystal structure analysis and pharmacophore modeling and also results of VS, docking, ADME, and MD simulations. The last section is a concluding remark on this developed pharmacophore models and resulted selective PDE4B selective inhibitors.

Software and workstation
All computations were performed on Fujitsu linux workstation (xeon quad-core E3-1220 processor). Phase 3.9, Glide 6.3, LigPrep 3.0, Prime 3.6, Impact 6.3, and Qik-Prop 4.0 of Maestro 9.8 (Schrödinger, LLC, 2014) were used to generate pharmacophore hypothesis, docking studies and ADME studies. Molecular dynamics (MD) study was done with AMBER 14 on 320 processor SUN Microsystems cluster. The analyses of MD simulations were carried out by CPPTRAJ module of AMBER 14 and Visual Molecular Dynamics 1.9.1 tool (VMD). The crystal structures of human PDE4B and PDE4D were retrieved from RCSB Protein Data Bank (PDB) and were utilized for molecular modeling studies.

Comparative crystal structure analyses
The coordinates of different PDE4B structures were retrieved from PDB. Proteins were prepared using protein preparation wizard of Schrödinger.Inc. The residues of chain A of the homomeric structures were considered for analyses. The structures with reasonable resolution (≤2 Å) were shortlisted. 1Y2H with resolution of 2.4 Å was included as all other group member has low resolution (1Y2J: 2.55 Å). 3HMV was considered for the analyses as it might provide significant depth to the comparative study as it contains the C-terminal and nitrobenzamide substituent in crystal ligand. In addition, structures with extra resolved residues (4X0F and 3G45) were also considered. Shortlisted monomeric structures were prepared using protein preparation wizard (Impact 6.3) of Schrödinger.Inc. All shortlisted and prepared structures were aligned to apo-protein conformation (PDB ID: 1F0J) using Pymol software (DeLano, 2002). H-bonds and hydrophobic interaction between protein and co-crystal ligands were evaluated in PDB ligand explorer 4.1.0 (Moreland, Gramada, Buzko, Zhang, & Bourne, 2005).

Structure-guided de novo design
Shortlisted PDE4B structures were grouped based on the chemical substructures of the bound inhibitors. Within each group, a representative molecule (Rn) was designed that possessed different possible chemical features of all inhibitors in that group. To predict binding affinities, these newly designed representative molecules were docked in previously prepared protein structures. Likewise, using chemical features of representative molecules of each group, next single representative molecule (F) was designed by trial and error and its binding affinity in PDE4B (PDB ID: 3G45) was assessed using molecular docking.
Next, this final representative molecule F was submitted to Phase 3.9 to develop the ligand-based pharmacophore hypotheses. Thus, generated diverse, ligand-based pharmacophoric hypotheses were used to perform VS.
To validate the generated model, a common data-set of 30 active molecules from reported literature and 254 inactive molecules were used. Inactive molecules were chosen such that they consists of similar physiochemical descriptors (molecular weight, number of rotational bonds, hydrogen bond donor count, hydrogen bond acceptor count and the octanol-water partition coefficient) to active molecules, but deprived of any of the chemical descriptors of the active ligands. Various statistical parameters such as accuracy, precision, sensitivity, specificity and enrichment factor (E value) were calculated for each hypothesis.
Furthermore, for the analysis of results, E value score was calculated using the following formula E ¼ TP Â D Ht Â A where D, A, Ht, and TP represent the total number of compounds of the database, total number of actives, total number of hits and total number of true actives in the hits, respectively.

Pharmacophore-based VS
The novel and potential leads can be predicted from large chemical databases using VS techniques (Kumar, Garg, & Bharatam, 2015;Sengupta, Roy, & Bandyopadhyay, 2015). The 3D pharmacophore model represents the chemical functionalities responsible for the bioactivities, and hence the validated 3D queries (pharmacophore hypotheses developed in previous step) were used. SPECS database containing 212,783 molecules was used for VS. In first stage of VS, molecules were prescreened based on simple physico-chemical properties, i.e. the molecules with number of rotatable bonds of ≤7 and not obeying Lipinski's rule of five for rejected from current analysis. These filtered molecules were and then submitted to Phase 3.9 to pregenerate conformers to fasten screening process. These pre-selected, prepared molecules were next screened against previously generated pharmacophore hypotheses and hit molecules were docked into the PDE4B (PDB ID: 3G45) and PDE4D (PDB ID: 3G4G) structures.

Molecular docking
Crystal structures were prepared using Protein Preparation Wizard (Impact 6.3, Schrodinger) as follows: bond orders were assigned, hydrogen atoms were added, and formal charges were treated. Co-crystallized water molecules were not conserved in different PDE structures and therefore all water molecules were deleted. Hydrogen bonding network was optimized using the exhaustive sampling option. The protein was minimized to an RMSD limit from starting structure of 0.3 Å using Impref module of Impact with OPLS_2005 force field. Missing side chains were structured using Prime 3.6.
Ligands were prepared using LigPrep 3.0 with Epik 2.8 to expand protonation and tautomeric states at 7.0 ± 2.0 pH units. Specified chiralities were retained and at most 32 stereoisomers per ligand were generated keeping low energy ring conformation 1 per ligand. Molecules and proteins were minimized using the OPLS 2005 force field.
Docking grids were generated with the default settings in Glide using the co-crystal ligand to define the center of the grid box (20 Å × 20 Å × 20 Å). Default parameters for van der Waals radius scaling factor 1.0, partial charge cut-off 0.25, van der Waals radius scaling per atom scaling factor 1.0, charge scale factor 1.0 were used and no constraints were included during grid generation.
Prepared small molecules were docked into the protein structure using Glide 6.3 XP docking. "Write XP descriptor information" was chosen during the docking run which deduced energy terms such as hydrogen bond interactions, electrostatic interaction, hydrophobic enclosure, and π-π stacking interactions and the rest of the parameters like "dock flexibly," "add epik state penalties to the docking score" were kept default. The 3D complex structures of all hits docked in PDE4B and PDE4D were analyzed for Glide score and H-bonding interactions. The molecules with considerable difference in docking score and H-bonding interaction in PDE4B and PDE4D were selected and evaluated for binding free energy calculation of docked complex using Prime-MM-GBSA-3.6, employing VSGB continuum dielectric model as solvent model.

Absorption, distribution, metabolism, and excretion
The analysis of ADME (Absorption, Distribution, Metabolism and Excretion) properties are crucial determinants for the successful development of new drugs, and are now imperative for rational drug design. The failure of which might lead to rejection of a drug in the later stages of drug process. All the shortlisted hits were further processed for ADME properties analysis using Qik-Prop 4.0 tool of Schrodinger. QikProp is built using experimental details of 710 compounds including 500 drugs and heterocyclic compounds. It calculates properties like molecular weight, molecular volume, no. of H-bond donors, no. of H-bond acceptors, polar surface area, predicted octanol/water partition coefficient (QPlogPo/w), and violations related to Lipinski's "Rule of 5" and Jorgensen's "Rule of 3" to filter out compounds with clear cut undesirable properties.

Molecular dynamics simulation
Among newly identified hit molecules, the stability of protein-ligand interactions of docked AO-022/43390796 (AO) in PDE4B (PDB ID: 3G45) was evaluated by molecular dynamics simulations. Topology and parameter files for the protein residues were generated using "ff99SB" force field in AMBER 14 (Lindorff-Larsen et al., 2010). The missing parameters for the small molecule, Zn and Mg were generated using the general amber force field (GAFF) and the Parmchk module of Antechamber. All systems were simulated in a standard protocol of energy minimization which is 2500 steps of steepest descent followed by 1000 conjugate gradient steps. Thus, minimized system was gradually heated from 0 to 298 K in 200 ps. Following equilibration, MD simulation was carried out for 10 ns with periodic boundary conditions in an isothermal-isobaric (NPT) ensemble at a temperature of 298 K with Berendsen temperature coupling and a constant pressure of 1 atm with isotropic molecule-based scaling. SHAKE algorithm was applied to fix all covalent bonds containing hydrogens and particle-mesh-Ewald method was used for long-range electrostatic interactions. The analyses of MD simulation were carried out by CPPTRAJ module of AMBER 14 and VMD visualizer (Humphrey, Dalke, & Schulten, 1996). All MD simulations were performed on 320 processor SUN Microsystems cluster. The active site residues are mostly conserved between PDE4B and PDE4D. Beside catalytic Gln615 (Gln535 in PDE4D), other H-bonding residues in PDE4B active site like Tyr405 (Tyr325 in PDE4D), His406 (His326 in PDE4D), His410 (His330 in PDE4D), Asp447 (Asp367 in PDE4D), Asp564 (Asp484 in PDE4D), Asn567 (Asn487 in PDE4D) are also conserved between two closely related homologs. Likewise, the hydrophobic residues are also conserved (Table 1).

Results and discussion
Besides structurally resolved catalytic domain, full PDE4 proteins (all PDE4A-4D) also consist of upstream regulatory regions (UCR1-2). Recently, Burgin et al. (2010) showed that UCR2 closing over to PDE4B active site region prevent access to substrate ATP Tyr274 of PDE4B-UCR2 is occupied by Phe196 of PDE4D-UCR2 which could make H-bond with the bound inhibitor and therefore provide selective PDE4B inhibition, which can be explored to design selective PDE4B inhibitors. Therefore, in our study we used structures with these additional resolved residues (PDE4B: 3G45:: PDE4D: 3G4G) to perform docking and molecular dynamics study.

Structure-guided de novo design
The diverse chemical scaffolds of different PDE4B co-crystal inhibitors were used to design ligand-based pharmacophore. Among 34 structures, crystal structures with resolution of ≤2 Å were considered. In addition, four low-resolution crystal structures were included: 4X0F, 3G45 with additional N-terminus; 1Y2H with resolution of 2.4 Å within all group members (1Y2J-2.55Å) and 3HMV (additional C-terminus and nitrobenzamide substituent of crystal ligand) that might provide significant depth to the comparative crystal study. All 14 refined PDE4B structures ( Figure 2) were aligned to apo-protein conformation (PDB ID: 1F0J) to critically analyse the distinct changes in the conformations of the bound inhibitors.
To identify and enumerate all possible pharmacophoric features of active molecules, different co-crystal ligands were divided into six groups on the basis of their chemical substructures. Within each group, a representative molecule (R 1-6 ) was designed which showed all possible interactions of different grouped ligands (Figures 2  and 3). The interactions of designed representative molecules in PDE4B structures were evaluated by molecular docking and compared with in-place docking score of crystal ligands. These group representatives are as follows: Group 1 (R 1 ): This group consisted of two co-crystal ligands: 5RM (PDB ID: 1XM6) and ROL (PDB ID: 4X0F). Similar to these ligands, representative molecule, R 1 possessed conserved phenyl ring which makes π-π interaction with Phe618. The conserved H-bond with Gln615 was maintained by keeping two alkoxy groups. Further to enhance interactions in hydrophobic region (Ile582, Phe586, Leu682, and Tyr405) a terminal cyclopentyl was incorporated. Additional 4-pyrrolidine-2-onering at p-anisole position assisted H-bonding with Tyr274.
Group 3 (R 3 ): 6DE (PDB ID: 1Y2H), the only member in the group was considered as a group representative molecule. 6DE consisted of a pyrazole ring that captured π-π interaction with Phe618. The carbonyl group of ethyl carboxylate substituent at C4 made Hbond with Gln615. 2-chloro-1-phenyl-substitution at N2 form hydrophobic interactions particularly with Ile582.
Group 4 (R 4 ): The common chemical feature of this group ligands was a six-membered, 2-3 nitrogen containing aromatic ring such as pyrimidine in 19T (PDB ID: 4MYQ) and triazine in 1S1 (PDB ID: 4KP6). Such ring systems were favorable to make π-π interactions with Phe618 and therefore, pyrimidine ring was incorporated in R 4 . 2-pyrolyl substitution at C2 promoted H-bond with Gln615. Substituting ethyl at C6 and amino phenylacetamide at C5 (instead of phenyl acetic acid as it was in 1S1) which forms hydrophobic interactions with Ile582 and Phe586 were included in R 4 also. C4 position was substituted with amino group that captured H-bond with Asn567.
Group 5 (R 5 ): This group members HBT (PDBID: 3HMV) and 3OJ (PDBID: 3O0J) consisted a hetero-bicyclic ring system: one aromatic and other aliphatic ring with heteroatom (S, B or O) that made π-π interaction with Phe618. Therefore, R 5 was designed to have a similar bicyclic ring system. However, we incorporated Natom in the aromatic ring to make H-bonding with Gln615. We incorporated CONH 2 group at C5 to make H-bond with His406 or Asn567. Further, benzamide substitution at C6 position captured hydrophobic ambience with Ile582 and Phe586. Position C2 was substituted with methoxybenzene to capture an additional H-bond Gln615.
Adding up all features of six representative ligands (R 1-6 ) from six groups, we next designed a final representative molecule F (Figure 4). From R 1 , cyclopentyl methyl ether substituent was considered. Pyrazolo pyridine bicyclic moiety substituted at C3 was grasped from R 2 . 1H-Pyrazole was considered from R 3 . Aminobenzamide substituent modified as 3,4-dihydro-2H-pyrrol-3-yl[4-(methylamino)phenyl] methanol substituent at C5 was adopted from R 4 . Methoxypyridine was gained from R 5 . Pyridine ring, amido substituent at C3 and methylamino phenyl at C5 was considered from R 6 . Molecular interactions of this final representative molecule in PDE4B (PDB ID: 3G45) are shown in Figure 5. This final representative molecule F when docked in PDE4B structure (PDB ID: 3G45), made H-bonds with catalytic Gln615. Other H-bonding residues include Tyr405, Thr517, and Asn567. Further, F showed H-bond with PDE4B isoform specific Tyr274. In addition, F made H-bonds with three co-crystal water molecules. Metal interaction with bound Mg 2+ was also observed. Glide score of F in PDE4B was −13.9.

Ligand-based pharmacophore modeling and VS
Structurally different molecules having diverse pharmacophoric features can exert similar biological activity. The pharmacophore searching is one of the rational approaches for drug design based on the maximum information available. Based on this assumption, ligand-based pharmacophore hypotheses were generated using protocol adopted from Phase module, (Maestro 9.8), followed by VS of commercial databases based on these hypotheses.
Prior to VS, generated hypotheses were assessed using data-set of active (30) and inactive (254) ligands. Data-set from literature was employed to validate the generated pharmacophore hypotheses. To each models, different statistical parameters like accuracy, precision, sensitivity, and specificity of the best pharmacophore models were calculated. Furthermore, an E value of .74, 1.31, and 1.05 was calculated for model 1, 2, and 3, respectively (Supplementary Table 3). From the overall validation results, we assured that hypotheses can differentiate between the actives and inactive molecules.
These three hypotheses were individually submitted to find hit molecules in Phase against Specs database. In an early filter, "non-hit" molecules were screened which saved the computational cost during molecular docking. Remaining 203,325 (12,458 omitted) SPECS molecules were screened against the three hypotheses. The maximum hit number per hypothesis search was set to 1000. The 3000 hits obtained from the 3 VS, were docked in PDE4B (PDB ID: 3G45) and PDE4D (PDB ID: 3G4G).

Molecular docking and ADME analyses
Docking study was performed using Glide program on previous hit molecules to identify the key interactions with PDE4B structure. N-terminal residues of PDE4B provide selectivity to the bound inhibitors. Therefore, all molecules were docked in PDE4B structure with resolved N-terminus residues (PDB ID: 3G45). Glide computes the binding energy (BE) in terms of the Glide docking score with respect to the docked ligands. Among all hits, molecules that have Glide score of less than -8 were only considered for further analyses.
Most of the current PDE4B inhibitors are associated with PDE4D related adverse effects. To minimize PDE4D related toxicity, molecules obtained with Glide score <−8 were next docked in PDE4D structure (with resolved N-terminus residue, PDB ID: 3G4G) and molecules with similar (or high) docking score in PDE4D structures were rejected from further analyses. Molecules with minimum Glide score difference of >2 (PDE4D < PDE4B) were only chosen as selective PDE4B inhibitors.
Whereas the intrinsic ligands cAMP, and reported inhibitors such as rolipram, roflumilast, and 3-isobutyl-1methylxanthine (IBMX) showed similar Glide score, H-   bond, and hydrophobic ambience in both PDE4B and PDE4D, the screened ligands showed better selectivity profile for PDEB in terms of Glide score, H-bond and π-π interaction ( Table 2).
Twenty-seven ligands from hypothesis-1, 25 ligands from hypothesis-2, and 28 ligands from hypothesis-3 were finalized on the basis of considerable difference in docking score of both PDE4B and PDE4D. Nine ligands from hypothesis-1, 14 ligands from hypothesis-2, and 10 ligands from hypothesis-3 were selected for considerable difference in prime binding energy of PDE4B and PDE4D-docked complexes.
The ADME properties of all the obtained ligands were assessed using Qikprop 4.0, and these ADME properties of best hits are listed in Supplementary  Table 4. The 20 of the hits showed optimum ADME properties while 13 were rejected. The ADME property of these hits makes them promising candidates for as PDE4B inhibitors. A schematic of VS approach is shown in Figure 7. The selective inhibitors, their Glide scores along with the hydrogen bonding of key residues are presented in Table 2 and Supplementary Table 5.

Molecular dynamics simulation
Protein-ligand (PDE4B: AO) complex was simulated for 10 ns. Figure 8(A) shows the temporal RMSD plot of protein residues and ligand for the last 4 ns MD simulation run. The curve shows that the RMSD values of both protein and ligand have very small variation. During the earlier phase of the simulation the variation is more which reduces and remains approximately same 6000 ps onwards. The deviation is minimal during last 4 ns and therefore, the trajectories between this time periods (7-10 ns) are used to perform MD analyses and to extract the minimum energy conformation. Figure 8(B) shows minimum energy structure of PDE4B-AO complex for  last 4 ns run. Comparison of starting docked conformation to this minimum energy conformation revealed that H-bonds with Gln615 and Tyr274 remain stable during MD simulation (Figure 8(B), Supplementary Figures 1  and 2). During simulation, two intermediate duration (~2 ns) H-bonds with side chain of Tyr405 were also formed. Apart from H-bonds, hydrophobic interaction with Phe618 was consistent during MD simulation (Figure 8(B) and Supplementary Figures 3-5). The stable protein-ligand interactions during 10 ns MD simulation support our previous docking results and confirmed the role of new molecules as potent PDE4B-inhibitors.

Conclusion
With the well-established role of PDE4B inhibitors in COPD, selective PDE4B inhibitor drug research has become crucial. This work aims to contribute in this direction. Based on the knowledge of available PDE4Binhibitor complex structures from RCSB protein data bank, pharmacophore models were built representing chemical feature for PDE4B inhibition. Such pharmacophore models were used to screen lead-like molecular library of SPECS database. Virtually screened molecules were further subjected to docking and ADME analysis. Twenty molecules with better binding affinities to PDE4B than PDE4D were selected as probable lead-like molecules for PDE4B inhibitor design. Top scoring ligand AO-022/43390796 in complex with PDE4B was subjected to 10 ns MD simulation. Results of the simulation studies showed stable hydrogen bonds and hydrophobic interaction with the catalytic residues highlighting stable binding of the ligands to the protein.
Results reported in this article are likely to help the medicinal chemists to develop new drugs for COPD.

PDE Phosphodiesterase cAMP
Cyclic adenosine monophosphate COPD Chronic obstructive pulmonary disease PDB Protein data bank VS Virtual screening MD Molecular dynamics ADME Absorption distribution metabolism excretion RMSD Root mean square deviation

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2016.1183520.