Discover natural compounds as potential phosphodiesterase-4B inhibitors via computational approaches

cAMP, intracellular cyclic adenosine monophosphate, is a ubiquitous second messenger that plays a key role in many physiological processes. PDE4B which can reduce the cAMP level by hydrolyzing cAMP to 5′-AMP has become a therapeutic target for the treatment of human diseases such as respiratory disorders, inflammation diseases, neurological and psychiatric disorders. However, the use of currently available PDE4B inhibitors is restricted due to serious side effects caused by targeting PDE4D. Hence, we are attempting to find out subfamily-selective PDE4B inhibitors from natural products, using computer-aided approaches such as virtual screening, docking, and molecular dynamics simulation. Finally, four potential PDE4B-selective inhibitors (ZINC67912770, ZINC67912780, ZINC72320169, and ZINC28882432) were found. Compared to the reference drug (roflumilast), they scored better during the virtual screening process. Binding free energy for them was −317.51, −239.44, −215.52, and −165.77 kJ/mol, better than −129.05 kJ/mol of roflumilast. The pharmacophore model of the four candidate inhibitors comprised six features, including one hydrogen bond donor, four hydrogen bond acceptors, and one aromatic ring feature. It is expected that our study will pave the way for the design of potent PDE4B-selective inhibitors of new drugs to treat a wide variety of diseases such as asthma, COPD, psoriasis, depression, etc.


Introduction
Intracellular cyclic adenosine monophosphate (cAMP) is an important second messenger involved in signaling pathways. Due to its ubiquitous distribution in human bodies, it plays crucial roles in regulating many physiological processes such as inflammation, cell growth, and sensory signaling (Li et al., 2013). The concentration of cAMP in cells is tightly regulated at the point of synthesis by adenylyl cyclase and at the point of degradation by enzymes from cyclic nucleotide phosphodiesterases (PDEs) (Xu et al., 2000). Type 4 phosphodiesterase (PDE4) enzymes are taken as critical cAMP hydrolyzing enzymes in cells and, as such, are responsible for the inactivation of cAMP via hydrolyzing it to 5′-AMP (Chen et al., 2012). Because of their wide tissue distribution, the PDE4 family acts as a critical regulator of intracellular cAMP levels, cAMP signaling, and signal compartmentalization (Jin, Ding, & Lin, 2012).
PDE4 is a key isozyme expressed by immune cells and plays a critical role in degrading cAMP in, for example, lymphocytes, monocytes, macrophages, and microglial cells (Bender & Beavo, 2006). The cAMPspecific PDE4 comprises four members, namely PDE4A, PDE4B, PDE4C, and PDE4D. They differ in their pattern of expression in the body and differ in their pattern of targeting subcellular compartments (Fox, Burgin, & Gurney, 2014). PDE4B that regulates inflammation via cAMP degradation has been a well-validated target for therapeutic immunomodulation and the focus of extensive drug development programs for treatment of autoimmune diseases (Hatzelmann et al., 2010;McCann et al., 2010). Owing to the increase of PDE4B drives down cAMP levels, causing a rise in production of tumor necrosis factor alpha in response to the inflammatory stimulus mediated by the Toll receptor pathway, elevating the cAMP level attenuates inflammatory responses to treat such disorders (Jin & Conti, 2002;Jin, Lan, Zoudilova, & Conti, 2005). Moreover, the elevated level of cAMP contributes to therapeutics of respiratory diseases, neurological and psychiatric disorders, hinting PDE4B is an attractive target for investigation (Chen et al., 2011).
Indeed, PDE4 as a therapeutic target has gained much attention and there are more PDE4 inhibitors than that have been developed for any other PDE families. These inhibitors have been extensively studied as antiinflammatory therapeutics for treatment of asthma and COPD (chronic obstructive pulmonary disease). Suffer from other diseases such as rheumatoid arthritis, septic shock, and atopic dermatitis can also be alleviated by PDE4B inhibitors (Wang et al., 2007). Besides, the potential use of PDE4 inhibitors for treatment of depression and affective disorders such as schizophrenia has also been studied (Kanes et al., 2007;O'Donnell & Zhang, 2004). However, these intense efforts toward the development of effective and selective PDE4 inhibitors have not achieved much success yet. For example, rolipram causes emesis at effective anti-inflammatory doses, and roflumilast is limited by nausea and emesis at high doses (Horowski & Sastre-Y-Hernandez, 1985;Spina, 2003). The side effects may be caused by the structural similarity between two isoforms of PDE4, PDE4B (therapeutic effect), and PDE4D (side effect of emesis) (Sharma, Wakode, Lather, Mathur, & Fernandes, 2011). As PDE4 isoforms have unique nonredundant roles, PDE4 isoform-selective inhibitors in terms of PDE4Bselective inhibitors may offer a route for maximizing therapeutic actions while minimizing side effects (Wang et al., 2007).
Therefore, PDE4B as the target for PDE4 isoformselective inhibitors was used in our study, for the objective of finding a new generation of PDE4 inhibitors with increased efficacy and reduced restriction of side effects. Computer-aided drug design methods were performed to find potential PDE4B inhibitors from natural products. We first implemented virtual screening for the discovery of candidate PDE4B-selective inhibitors, which were then evaluated by molecular dynamics (MD) simulation and binding free energy calculation. In addition, the pharmacophore model of potential PDE4B inhibitors was calculated, and its comparison with pharmacophore of existing inhibitors was carried out as well.

Data collection
The structure of human PDE4B was downloaded from RCSB protein data bank (PDB) via accession number 1XMU (Berman et al., 2000). The PDE4D structure that was used to test whether our newly found compounds are more likely to target PDE4B was downloaded from PDB using 3G4L (Burgin et al., 2010). In the 1XMU structure, PDE4B is in complex with roflumilast (Card et al., 2004). Chimera was then used for the protein preparation process, during which water molecules and noncomplexed ions were deleted, hydrogens were added, protonation states for histidine were detected automatically, and charges were assigned to residues (Pettersen et al., 2004). Hydrogens and charges were added to the co-crystallized roflumilast which was used as a reference drug for virtual screening and was saved as another structure file in format of mol2.
The dataset of natural products used for screening was downloaded from ZINC, a free database of commercially available compounds in ready-to-dock formats for virtual screening (Irwin, Sterling, Mysinger, Bolstad, & Coleman, 2012). AnalytiCon Discovery NP is one of the ZINC's natural products catalogs and provided to ZINC by AnalytiCon Discovery under collaborative agreements between them. The catalog has been refined for drug discovery relevance and according to its latest update on 17 February 2013, and there are 11,247 compounds within it. We downloaded this catalog on 5 August 2014 and used it for the discovery of potential selective PDE4B inhibitors through virtual screening.

Virtual screening
In order to acquire dependable results, virtual screening was performed by two tools, DOCK 6 and AutoDock Vina (also wrote Vina) (Lang et al., 2009;Trott & Olson, 2010). They are open-source programs for molecular docking, identifying the low-energy binding modes of a small molecule (or ligand) within the active site of a macromolecule (or receptor) whose structure is known. Taking roflumilast as a reference drug, we implemented virtual screening via DOCK and Vina, respectively, finding out molecules whose binding affinity with PDE4B was better than roflumilast. After that, compounds that ranked top 10 within the two result groups from DOCK and Vina were selected as potential PDE4B inhibitors.
In terms of DOCK 6, the molecular surface of the receptor (PDE4B) was generated with a changeable probe radius of 1.0. The spheres surrounding PDE4B were generated with a maximum and minimum sphere radius of 3.0 and 1.0 Å. A subset of spheres that represents the binding site of PDE4B was identified by the location of the co-crystallized roflumilast, within 10.0 Å RMSD (root mean square deviation) from every atom of it. The box around the active site was created by enclosing 5 Å extra margin in all six directions. Finally, virtual screening by DOCK was performed, and molecules were scored and subsequently sorted.
Regarding AutoDock Vina, the search space at the binding pocket was a rectangle of 40 Å × 32 Å × 40 Å, in the X, Y, and Z dimension. For every compound docked by Vina, three binding modes were generated, with each mode acquiring a score of binding possibility toward PDE4B. After the completion of the virtual screening, the best mode for each of these molecules was used for hits selection.

Enrichment calculation
There are a number of approaches to quantitating the success of a docking protocol for virtual screening. The most often used and simplest to calculate is enrichment at a given percentage of the database screened (Hawkins, Warren, Skillman, & Nicholls, 2008). In order to determine how capable DOCK and Vina are at identifying known active compounds embedded in a random database, 25 actives (i.e. known PDE4B inhibitors) and 500 decoys are docked to PDE4B and the results are rank ordered. Good enrichment is achieved when greater numbers of actives are ranked earlier in the list compared to the decoys (Jiang & Rizzo, 2015).
The 25 existing PDE4B compounds (Supplementary Table S1) as known actives were downloaded from the PDB database after manually selecting crystal structure of PDE4B-inhibitor complexes. The decoy set that comprises 500 compounds was created using DecoyFinder which selects, for a collection of active ligands of a protein target, a set of decoys from a database of compounds (Cereto-Massagué et al., 2012). Decoys are selected if they are similar to active ligands according to five physical descriptors (molecular weight, number of rotational bonds, total hydrogen bond donors, total hydrogen bond acceptors, and the octanol-water partition coefficient) without being chemically similar to any of the active ligands (Cereto-Massagué et al., 2012). Search settings for DecoyFinder are as follows: (1) Active ligand vs. decoy Tanimoto threshold is .75; (2) decoy vs. decoy Tanimoto threshold is .9; (3) hydrogen bond acceptors range is 2; (4) hydrogen bond donors range is 1; (5) log P range is 1.0; (6) molecular weight range is 25; (7) ROTATIONAL bonds range is 1; (8) Number of decoys per active ligand is 20.

MD simulation
All MD simulations for PDE4B-inhibitor complexes were performed using the GROMACS 4.5 package under tip3p water mode and Amber ff99sb force field (Hornak et al., 2006;Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983;Pronk et al., 2013). The system was simulated under the periodic boundary condition with a dodecahedron periodic box also solvated by tip3p water molecules, setting the surface of the protein covered with a water shell of 1.0 nm (Berendsen & Marrink, 1993). The number of sodium and chloride ions needed to balance the system charge was placed in the solvated system with salt concentration of .15 mol/L. The energy minimization was performed to relax the complex system under 1000 kcal/mol/nm using the steepest descent method . Afterward, 100 ps NVT (constant number of particles, volume, and temperature) ensemble stabilized the system at 300 K, and subsequently, 100 ps NPT (constant number of particles, pressure, and temperature) equilibration stabilized the system pressure using coupling reference pressure of 1 bar. Finally, we ran the 10 ns MD simulation with a time step of 2 fs, storing corresponding coordinates every 2 ps. During the simulation, long-range electrostatics was calculated by particle mesh Ewald algorithm (Essmann et al., 1995), and all bonds were constrained by linear constraint solver (LINCS) algorithm (Hess, Bekker, Berendsen, & Fraaije, 1997).

MM-PBSA binding free energy calculation
The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method has been widely used in binding free energy calculation. In this approach, the binding energy for a complex including a receptor and a ligand was calculated by the equation as following (Zoete, Irving, & Michielin, 2010), where ΔG bind is the binding free energy, ΔE gas symbolizes gas-phase energy, ΔG solv represents solvation free energy, and TΔS is the entropy contribution.
The gas-phase contribution denoting the MM part (ΔE MM ) of MM-PBSA is a sum of bonded terms (bond, angle, and torsion energy) and nonbonded terms (van der Waals and electrostatic interactions). The solvation contribution that consists of polar solvation (ΔG polar ) and nonpolar solvation (ΔG nonpolar ) energies is a symbol of the PBSA part. In practice, excluding the entropy term is acceptable provided that similar ligands are docked into the same binding pocket of a receptor (He, Li, & Wu, 2014). Therefore, the binding free energy is evaluated as, The GROMACS tool g_mmpbsa was used in this study for binding energy calculation for protein-ligand complexes (Kumari, Kumar, Consortium, & Lynn, 2014). Of note, only electrostatic and van der Waals components of the gas-phase energy are calculated by g_mmpbsa because the bonded contribution is by definition zero in the singletrajectory approach (Homeyer & Gohlke, 2012). Here, g_mmpbsa calculated binding energy in terms of component ΔE MM , ΔG polar , and ΔG nonpolar for each complex, through eleven snapshots extracted every .2 ns from the 8 to 10 ns MD trajectory. Within the calculation, the nonpolar energy was calculated using solvent-accessible surface area (SASA) model and other parameters for g_mmpbsa were set to default values as described in the developers' publication (Kumari et al., 2014).

Pharmacophore prediction
PharmaGist, a freely available web server for pharmacophore detection, was used for the exploration of pharmacophore model from a given set of molecules (Schneidman-Duhovny, Dror, Inbar, Nussinov, & Wolfson, 2008). The minimum number of features in pharmacophore was set to 3. Weights of aromatic ring, hydrogen bond donor/acceptor, charge (anion/cation), and hydrophobic were 3, 1, 1.5, and .3, respectively. Apart from calculating pharmacophore model of potential PDE4B inhibitors discovered by us, the 25 existing inhibitors' pharmacophore was estimated as well.

Validation of DOCK and Vina
Re-docking a co-crystallized ligand back into its corresponding receptor protein is an efficient way to validate docking methods, and it is also a straightforward way adopted widely (Morris et al., 1998). As an index monitoring how much the docking pose deviates from its original structure, RMSD whose value is less than 2.0 Å ensures the docking accuracy, subsequently proving it is trustworthy to dock other molecules within the binding site of the receptor using this docking protocol .
DOCK re-docked roflumilast back into the binding pocket of PDE4B and got a heavy atom RMSD of 1.13 Å. For roflumilast docked by Vina, the RMSD between its crystallized and re-docked modes was .76 Å. As can be seen from Figure 1, both DOCK and Vina successfully re-docked the crystallized roflumilast to the binding site of PDE4B, with the docked roflumilast almost totally superimposing upon its crystallized conformation. Accordingly, DOCK and Vina are able to explore conformational space sufficiently well to generate correctly docked poses.
After docking the ligand database comprised of 25 actives and 500 decoys to PDE4B by DOCK and Vina, we ranked the docked compounds by score and calculated enrichment for DOCK and Vina, respectively. According to Figure 2, Vina performs markedly better and DOCK is slightly better, demonstrating that both DOCK and Vina could identify active compounds from the ligand database since enrichment curves are better than the random selection curve. This indicates our docking method is accurate and therefore makes the virtual screening carried out by DOCK and Vina reliable.

Potential PDE4B inhibitors
Docked to PDE4B by DOCK, the reference drug roflumilast acquired a DOCK score of −29.15, which was used as a threshold to select hits from the virtual screening result of DOCK. Among the 11,247 compounds of natural products filtered by DOCK, 104 molecules got scores better than the reference drug (Supplementary  Table S2). For these hits, ZINC67912533 (−49.21) with the greatest binding affinity ranked first while ZINC67911372 (−29.16) with a relative better score than roflumilast ranked at the bottom of the table.
By docking the reference drug to PDE4B using Vina, roflumilast got a score of −10. According to this, 684 of 11,247 natural molecules were selected as hits from the Vina virtual screening result (Supplementary Table S2). Among the 684 compounds, ZINC03780340 is regarded as the best hit because it was marked −12.7 by Vina. The molecule ranked last in these hits is ZINC14652075, whose Vina score was −10.1, just a little better than that of the reference roflumilast.
Within hit molecules from DOCK and Vina, respectively, 31 compounds existed in both hit groups simultaneously. As is shown in Supplementary Table S3, DOCK scores of these molecules vary markedly (from −49.21 to −29.17), but their Vina scores seem to be similar   (between −10.1 and −12.1). ZINC67912533 ranks at the first place according to DOCK score on the one hand. On the other hand, Vina score shows ZINC67912780 should be the best. That means it is biased to select candidate PDE4B inhibitors based on hit group from either DOCK or Vina. Therefore, the ranking of these compounds in both groups was taken into account. We compared 10 top hits from the two groups and four molecules were observed to be presented in both of them, thus they were chosen as potent candidate inhibitors against PDE4B. It can be seen from Table 1 that ZINC67912770 has the best DOCK score (−42.2), but its Vina score (−11.6) takes the second place. DOCK scores for ZINC67912780, ZINC72320169, and ZINC28882432 are −38.28, −35.08, and −33.33, and their corresponding Vina scores are −12.1, −10.7, and −11.2, respectively.

Stability of PDE4B-inhibitor complexes
We can see from Figure 3 that RMSD of the PDE4Broflumilast complex experienced a sharp rise during the first 3 ns. After that, the RMSD leveled off at around .2 nm throughout the remained simulation time. However, it did not equilibrate well and seemed to fluctuate slightly during this period of time. Compared to the reference drug, RMSD of the ZINC67912770 system was much better. As displayed in Figure 3(a), RMSD of the complex quickly reached the platform after an upward trend of 2 ns. After that, it stood at about .125 nm and remained stable till the end of the 10 ns simulation. For the ZINC67912780-PDE4B complex, its RMSD equilibrated later but the value was lower than that of roflumilast (Figure 3(b)). The system fluctuated significantly before reaching its equilibrium around .15 nm by approximately 6 ns. Subsequently, it remained stable at this level. According to Figure 3(c), the ZINC72320169 system with a lower RMSD value than the reference drug got equilibrium by about 5 ns. It then remained stable at around .15 nm over the last 5 ns. Even though RMSD of the ZINC28882432-PDE4B complex tended not to be as well as the three other potential PDE4B inhibitors, Figure 3(d) illustrates ZINC28882432 can still be a little better than the reference drug. Its RMSD became stable after a fluctuation of nearly 4 ns and in the next 6 ns, the equilibrium was better, with a lower RMSD value standing at around .17 nm.

Binding modes
It can be seen from Figure 4 that at each check point within the 10 ns MD simulation, these candidate inhibitors bound to PDE4B with contribution from hydrogen bonding interactions. ZINC67912770 formed hydrogen bonds with His274, Asp275, Asn283, Gln284, Met347, and Gln443 when interacting with PDE4B. At the beginning, it only formed one hydrogen bond with Asp275. Two more hydrogen bonds were formed with His274 and Asn283 at 2 ns. The PDE4B-ZINC67912770 interaction formed four hydrogen bonds at 4 ns, one with the conserved residue Gln443, another with Met347, and two others with Asp275. The number of hydrogen bond reduced to three at 6 ns owing to losing interaction with Met347. At 8 ns, ZINC67912770 formed another hydrogen bond with Gln284, resulting in four hydrogen bonds in total. By 10 ns, there were three hydrogen bonds formed, one with Met347 and two with Asp275.
In the ZINC67912780-PDE4B interaction, there were at least three hydrogen bonds formed over the whole simulation. Asp346 and Asp392 were two residues dominating hydrogen bonding interaction with ZINC67912780. Initially, ZINC67912780 formed four hydrogen bonds with Glu300, Ser301, and Asp392. Then, Asp346 and Gln417 participated in the hydrogen bonding interaction at 2 ns. Three hydrogen bonds were formed at 4 ns as only Asp392 and Asp346 were involved in. Met431 was observed to form hydrogen bond at 6 ns, and this binding mode lasted to 8 ns. At 10 ns, a new residue, Tyr233, took part in forming hydrogen bond.
ZINC72320169 interacted with PDE4B and formed hydrogen bonds with Glu304, His307, and Asp392 at the beginning. Then at 2 ns, His274 replaced Glu304 in the hydrogen bonding network. By 4 ns, ZINC72320169 Figure 4. Interaction of each candidate inhibitors with PDE4B in the binding pocket during the 10 ns MD simulation. Potential PDE4B inhibitors are shown as ball and stick, with C colored tan. The light gray surface with 90% transparency is PDE4B. Residues that participate in the formation of hydrogen bonds are presented as sticks with red labels, with C colored yellow. Uniform color codes for other elements: O is red, N is blue, and H is green. lost hydrogen bonding with Asp392 and His307 but formed a hydrogen bond with Gln417. The inhibitor bound to PDE4B with the formation of three hydrogen bonds with His274 and Asp392 at 6 ns. The two residues kept hydrogen bonding interaction with ZINC72320169 by 8 ns. At 10 ns, ZINC72320169 interacted with His274 and His307, forming two hydrogen bonds.
The hydrogen bonding interaction between ZINC28882432 and PDE4B was not as good as the other three potential inhibitors because there was no hydrogen bond observed at 6 and 10 ns. ZINC28882432 formed two hydrogen bonds with His278 and Asp392 at the start of the simulation. Then, it only formed one at 2, 4, 8 ns, with Asp392, Met347 and Tyr233, respectively.

Binding free energy
As demonstrated in Table 2, all four candidate PDE4B inhibitors got a better binding energy than roflumilast. Binding free energy of ZINC67912780 (−317.51 kJ/mol) was the best. ZINC67912770 with a binding energy of −239.44 kJ/mol ranked at the second place. The binding energy of ZINC72320169 was relatively lower, with −215.52 kJ/mol. For ZINC28882432, its binding free energy was much lower than the first three inhibitors. With a binding energy of −165.77 kJ/mol, it was slightly better than roflumilast (−129.05 kJ/mol).
Regarding these four potential inhibitors, their gas-phase contribution, namely the combination of van der Waals energy and electrostatic energy, was favorable for inhibitor binding while the solvation contribution that comes from polar solvation and nonpolar energy was positive and unfavorable. Among the four compounds, the van der Waals energy and nonpolar energy were similar to each other. Nevertheless, electrostatic energy and polar solvation energy of ZINC28882432 was significantly different from the other three, with −38.59 and 136.34 kJ/mol, respectively. The weaker unfavorable polar solvation energy could have benefited ZINC28882432's binding to PDE4B, but the even weaker electrostatic energy consequently leaded to its reduction in binding affinity with PDE4B. Interestingly, unlike the four potential PDE4B inhibitors, the electrostatic energy of roflumilast was a positive value, which was opposite to its binding with PDE4B.

Discussion
Even though our validated virtual screening approach is dependable to carry out the research, the docking-scoreordered list of lead compounds displayed in Table 1 varies from DOCK to Vina. This is not a paradox because DOCK and Vina use different algorithm to assign docking score to docked ligands. DOCK uses physics-based Amber scoring function here which implements molecular mechanics GB/SA (generalized Born with SASA) simulations with the traditional all-atom AMBER force fields and the generalized AMBER force field (Lang et al., 2009). Vina uses a scoring function that combines certain advantages of knowledge-based potentials and empirically weighted scoring methods to do molecular docking (Trott & Olson, 2010). In this study, DOCK and Vina do the scoring process, respectively, so it is not surprising to find an identical lead compound ranked differently between the DOCK result list and Vina result list.
Nowadays, it is increasingly known that among the PDE4 family, PDE4B is a therapeutic target while PDE4D is a target with side effects. As a result, inhibitors that target PDE4D always cause side effects such as emesis and nausea, limiting their use in disease treatment (Sharma et al., 2011). As can be seen in Table 3, that these four newly found inhibitors are more likely to bind to PDE4B than PDE4D, which means they can be PDE4B-selective inhibitors. When docked to PDE4B, these candidate PDE4B inhibitors acquired better docking score than docked to PDE4D. In terms of DOCK score, the difference between PDE4B and PDE4D is noticeable.
The Vina score of ZINC67912770 rose by .2 to −11.4; ZINC67912780 went up from −12.1 to −11.6; ZINC72320169 increased by .3 from −10.7. Therefore, these results illustrate that the four potential PDE4B inhibitors prefer to target PDE4B, without or with low side effects caused by targeting PDE4D.
We can see from Figure 5 that the RMSF (root mean square fluctuation) profile of the four candidate inhibitors in complex with PDE4B displayed similar characteristics to RMSF of the roflumilast-PDE4B system. Two terminal regions and four middle regions of PDE4B were more flexible than other parts. It is not uncommon to find flexible protein terminals from RMSF of proteinligand complexes due to residues at both ends of the protein are seldom involved in the binding region, not participating in the binding interaction consequently. However, there are two remarkable differences presented between roflumilast and the four newly found PDE4B inhibitors. For these four promising PDE4B inhibitors, region Asp340-Arg380 became more flexible, whereas region Leu400-Ile450 was more stable, compared to the reference drug. As the binding modes in Figure 4 shows, these candidate PDE4B inhibitors mainly formed hydrogen bonds with residues from stable regions. Region Leu400-Ile450 has been known to be involved in hydrophobic interaction and hydrogen bonding interaction between PDE4B and inhibitors (Goto et al., 2013;Xu et al., 2000). So such difference suggests the region is more favorable for inhibitors binding to PDE4B. The more unstable region Asp340-Arg380 indicates researchers should not focus on this region when designing novel PDE4B inhibitors. Instead, PDE4B-selective inhibitors may come from by diminishing their interaction with residues in this region.
Secondary structure of region Lys328-Leu461 in Figure 6 clearly shows how the binding of inhibitors affected the protein's structure at the binding site. It can be seen that most of the PDE4B kept their secondary structure unchanged except two regions, Leu373-Arg385 and Thr397-Thr407, in correspondence with findings from the RMSF analysis above. In region Leu373-Arg385, the secondary structure changed between turn and bend frequently. This may be the reason for the more flexible RMSF profile of residues in this region.  The secondary structure for residues within region Thr397-Thr407 changed from turn to 3′-helixe and bend sometimes. According to their RMSF attributes, such change leaded to a favorable conformation for binding of PDE4B inhibitors. Noticeably, unlike other three candidate inhibitors, ZINC67912780 bound to the binding pocket of PDE4B and changed the second structure of residues within region Gly445-Ile450. Where it changed to bend in the first 5 ns but remained unchanged during the time left. Of note, RMSF of the ZINC67912780 complex in this zone was not affected, accompanying with a relative better binding mode and the best binding free energy. Hence, finding inhibitors that can interact with this region should be a good strategy in designing a new generation of PDE4B-selective inhibitors. Among these four potential PDE4B inhibitors, ZINC67912770 and ZINC67912780 are hydroxyphenyl acrylates while ZINC72320169 and ZINC28882432 are tetrahydropyrans. Pharmacophore of the four candidate PDE4B inhibitors is typically composed of six features: one hydrogen bond donor, four hydrogen bond acceptors, and one aromatic ring feature (Figure 7). For the 25 known PDE4B inhibitors, their pharmacophore model consists of three features, two hydrogen bond acceptors and one aromatic ring. Compared to this, pharmacophore of the four promising PDE4B inhibitors possesses another hydrogen bond donor feature and in addition, it has two more hydrogen bond acceptors. Actually, fewer pharmacophore features of these known inhibitors may result from more molecules were used in calculation. Basically, comparison between the two pharmacophore models implies pharmacophore of the four candidate inhibitors provides a clue for further development of potent PDE4B inhibitors.

Conclusion
It is common that the long tedious drug development process can take more than a decade of time and millions of dollars. Worse still, only a small proportion can finish clinical trials and be approved by drug regulatory agency, regardless of a large great number of new therapeutic candidates are being discovered in laboratories every year (Huang et al., 2010). Therefore, the pharmaceutical industry is subject to a widening gap between  productivity and the annual investment. Currently, many pharmaceutical corporations seek to increase short-term profitability via the integration of computational force and molecular biology. The computer-aided drug design that offers a fast track to some limiting factors such as low speed and high cost is complementing traditional drug discovery.
By computational methods integrating virtual screening, MD simulation, binding free energy calculation, and pharmacophore prediction, we found four potential PDE4B-selective inhibitors. Compared to currently used roflumilast, these inhibitors exhibited a stronger binding affinity with PDE4B. They are projected to elevate cellular cAMP level by inhibiting the effect of PDE4B, which hints their therapeutic potentials for the treatment of respiratory disorders, inflammation diseases, neurological and psychiatric disorders (Chen et al., 2011). Even though there is no experimental evidence to confirm the four molecules' bio-activity in inhibiting PDE4B in the current study, our findings are sufficient to demonstrate that they are potent PDE4B inhibitors. To be honesty, we do agree that to what extent the four compounds can be used as approved drugs for humans still requests further in vitro and in vivo testing.

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