Computational Approaches in the Development of Phosphodiesterase Inhibitors

Computational drug design tools have become indispensable in the quest for new drugs. There is hardly any drug discovery program where computational methods are not employed, be it structure‐based or ligand‐based methods. Numerous drug targets have been explored for discovery of new drugs using computational methods. In recent times, discovery of newer and selective phosphodiesterase as medications for inflammatory disorders, CNS disorders, and many other diseases has been the focus of many research groups worldwide. Most of these groups have employed computational methods of drug design and discovery at different stages of their research. This chapter reviews the reported application of computational methods used in the discovery and development of phosphodiesterase inhibitors.


Introduction
The early application of computational methods as a means to identify and design phosphodiesterase (PDE) inhibitors goes back to the 1980s where the first reported study tried to explain the relation between specific physicochemical properties and potency of known inhibitors [1]. The result was the identification of pharmacophore for PDE inhibitors [2]. Computational methods of drug design and discovery have impacted the overall process of drug discovery in a big way over the last two decades [3,4]. The cost as well as time of drug discovery has been reduced due to the inclusion of computational methods as well as highthroughput synthesis and screening. Computational methods have undergone development importance of atom-based topological and whole molecular E-state as well as 3D topological indices.
Hologram QSAR, a relatively newly developed QSAR technique, relates biological activity to structural fragments [15]. HQSAR eliminates the need for generation of 3D structures, putative binding conformations, and molecular alignments. The process of generating HQSAR models involve the fragmentation of the data set structures and then hashing into array bins. Molecular hologram fingerprints are then generated. Holograms are constructed by cutting the fingerprint into strings at various hologram length parameters. After the generation of descriptors, partial least square (PLS) methodology is used to find the possible correlation between dependent variable (activity) and independent variable (descriptors generated by HQSAR structural features) [15]. The best HQSAR model for cinnolines as PDE10A inhibitors was found to be statistically significant (q 2 = 0.664, R 2 pred = 0.513), and it highlighted some important structural features.
Pharmacophore mapping was also employed in this case. The pharmacophore hypothesis showed the importance of hydrogen bond acceptors and ring aromatic and hydrophobic features for higher activity. The compounds were mapped to the pharmacophore as per their activity, and the pharmacophore provided an efficient means of aligning the compounds for comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) studies.
Molecular dynamics (MDs) simulations were also performed for two ligand-receptor complexes, and the findings of MD simulations were consistent with the interpretations obtained from ligand-based analyses methods described above. The role of hydrogen bond acceptors and aromatic and hydrophobic features in determining PDE10A inhibition potential was unraveled. The study clearly puts forth guideline for the design and development of cinnolines as potent and selective PDE10A inhibitors and demonstrates the role of ligand-based drug design methods combined with structure-based drug design (SBDD) methods.

PDE4 inhibitors
PDE4 has long been considered as target for design of antiinflammatory molecules, and numerous PDE4 inhibitors have been developed [16][17][18][19][20][21]. Quinolines as a third-generation of PDE4 inhibitors show an increased antiinflammatory effect without being dose limited by side effects compared to the first-and second-generation of PDE4 inhibitors such as rolipram and roflumilast. Recently, newer quinoline derivatives were developed as selective PDE4B inhibitors using ligand-based pharmacophore and atom-based 3D-QSAR modeling along with structure-based docking and ADME methods [21].
These studies, in combination, led to new understanding about the selection of target molecules from many candidates. In this case, computational methods were not the primary tools used for designing the ligands rather they were used for identifying and confirming the binding site of the ligands on PDE4B. The 3D-QSAR model of PDE4B inhibitors developed for this study proved to be reliable with r 2 value of 0.96 and q 2 value of 0.91. The specific pharmacophore for PDE4B was then mapped and selected for virtual screening, and the potent PDE4B inhibitors were finally confirmed their selectivity ability for PDE4B by docking, ADME analysis, and MD, and then, molecules were developed as selective PDE4B inhibitors (Figure 1) [21].

PDE7 inhibitors
Like PDE4, PDE7, a cAMP-specific phosphodiesterase, is also highly expressed in human immune system such as thymus, lymph nodes, spleen, and blood leukocytes, suggesting PDE7 as a possible target for treating CNS and airway diseases [22][23][24][25]. A recent study by Cichero et al. whereby they applied 3D-QSAR methods for the development of selective PDE 7 inhibitors is worth noting. Herein, 72 derivatives of thieno [3,2-d]pyrimidin-4(3H)-one were selected from literature and were used to develop new selective PDE7 inhibitors [26]. The uniqueness of the study was the application of 3D-QSAR to identify the structural features required for the selectivity for a particular target.
Docking-based 3D-QSAR followed by redocking was performed to identify the most suitable bioactive conformations of derivatives and most comparable binding mode with that of X-ray structure. The conformations of the compounds showed good agreement with the conformations in the known PDE7A-ligand complex. The compounds were further submitted to a CoMFA and CoMSIA analysis. Since inhibitors which showed high selectivity for PDE7, rather than PDE4, are candidates for antiinflammatory molecules, two models (PDE7 selectivity (model A) and   Figure 2). The role of steric, electrostatic, and hydrophobic features in the binding mechanism of ligands to both PDE7 and PDE4 was identified. The study unraveled structural information for the development of new PDE7 inhibitors with high selectivity.

PDE11 inhibitors
Another unique application of ligand-based methods was reported by Cichero et al. earlier in 2013. Since the X-ray structure of PDE11 was not available, ligand-based homology modeling technique was applied to explore the 3D structure of PDE11 [27]. Consistent with the sequence similarity of PDE11 and PDE5, several tadalafil analogs (PDE5 ligands) showed binding affinity to PDE11 too [28][29][30]. Thus, the 3D structure of PDE11 was built on the basis of PDE5-tadalafil complex using homology modeling technique [31]. Besides the conventional homology modeling steps such as insertion and deletion of extra atoms during the energy minimization stages, evolution, including model building and refinement were performed in this specific case. Specifically, residues located in H-loop and M-loop were involved in obtaining PDE11 structure as they are important for substrate recognition [32]. The coordinates of PDE11 were derived from PDE5-tadalafil complex, and the amino acid sequences were obtained from the SWISSPROT database. This was followed by the development of models using MOE software without significant main chain deviations. AMBER94 force field was applied to minimize the structure energy, and Ramachandran plots were applied to assess the final obtained model. Successively, reliability of the derived PDE11 model was further confirmed by molecular docking and MD simulations with selective PDE11 inhibitors and dual PDE5-PDE11 inhibitors. This study provided new information of target structure of PDE11 to design selective PDE11 inhibitors for the treatment of cardiac pathologies [33]. The success of this approach is supposedly due to the synergic interaction between theory and experiment.

Structure-based methods: general aspects and perspective
Structure-based drug design (SBDD), as one of the in silico methods, is almost an integral part of any drug discovery and development project. These methods, in addition to the binding affinity between a specific protein target and a ligand, also provide insight into the interaction between the two. This helps in devising substituent modifications around the ligand scaffold leading to improved binding [34]. Based on the existing knowledge of 3D structures of PDEs, the great potential and success of structure-based computational methods have been visible in the development of newer PDE inhibitors. In addition to the routine screening and interaction studies, many unique applications of structure-based methods have been reported for the development of PDE inhibitors.

Selective PDE4B inhibitors
In one of the recent reports by Jing Li and team, virtual screening of natural product database was carried out to discover novel selective PDE4B inhibitors from natural products database [35]. Structure-based approaches like docking and molecular dynamics simulation were used for the purpose, although pharmacophore-based screening was also employed as initial filter. The screening led to the identification of four potential PDE4B-selective inhibitors (ZINC67912770, ZINC67912780, ZINC72320169, and ZINC28882432; Figure 3). Compared to the reference drug (roflumilast), they scored better during the virtual screening process. DOCK and Vina were used for docking, and results were agreeable. Binding free energy was −317.51, −239.44, −215.52, and −165.77 kJ/mol, which is better than −129.05 kJ/mol of roflumilast. The MD studies also showed that ZINC28882432-PDE4B complex reached stable RMSD faster than roflumilast-PDE4B complex. This study demonstrates the successful application of MD in a virtual screening workflow.

Selective PDE5 inhibitors
The cross-reactivity of PDE 5 inhibitors like sildenafil, vardenafil, and tadalafil with hERG1 is well documented. However, hERG1 is proven to be involved in the regulation of human ventricular myocyte action potential, thus blocking this hERG channel results in the function loss of the PDE5 inhibitors and further leads to serious life-threatening disorders and cardiovascular problems [36,37]. One of the recent studies explored this cross-reactivity using docking. Open and open-inactivated states of hERG1 potassium channel were used as protein structures and binding interaction patterns and affinity between the proteins and PDE5 inhibitors were studied [38]. Three different structure-based docking tools including GOLD, MOE, and AUTODOCK were applied to identify the binding interactions between Sildenafil and hERG1 channel. In open-state conformation of hERG1, both GOLD and MOE showed common interactions and further AUTODOCK and GOLD scoring analysis gave a 2.16 kcal/mol lower energy compared to open-inactivated state. In open-inactivated conformation of hERG1, three docking tools showed similar interaction patterns but MOE docking results were more specific in terms of Hbonding, distance, and key residues. To further confirm the key residues which are responsible for binding affinity to the hERG1 channel, in silico alanine mutagenesis study was performed, and new promising molecules were designed based on the interaction patterns and alanine mutagenesis studies. MD simulations were finally used to confirm that the complex of the new compounds with hERG1 channel is much less stable than that with PDE5. This study showed new approach toward the design of PDE5 inhibitors with lower affinity for hERG1 channels.
In another study dealing with cross-reactivity of PDE5 inhibitors, Kayık and team explored cross-reactivity of PDE5 inhibitors with PDE6 and PDE11 [39]. The major challenge of their study was designing novel PDE5 inhibitors with decreased cross-reactivity with PDE6 and PDE11. For this aim, the similarity-based virtual screening protocol was applied for the "clean drug-like subset of ZINC database that contained more than 20 million small compounds. Moreover, molecular dynamics (MD) simulations of selected hits complexed with PDE5 and off-targets were performed to get insights for structural and dynamical behaviors of the selected molecules as selective PDE5 inhibitors. Since tadalafil blocks hERG1 K channels in concentration-dependent manner, the cardiotoxicity prediction of the hit molecules was also tested. The study revealed important structural information for the design of novel, safe, and selective PDE5 inhibitors by applying SBDD.

Dual PDE1 and PDE5 inhibitors
Dual inhibitors of PDE1 and PDE 5 are known to produce vasodilation and have thus been explored as therapeutic agents for the treatment of different cardiovascular diseases, i.e., hypertension, angina, heart failure, and arteriosclerosis [40,41]. Yamazaki et al. employed a ligand-based virtual screening for the identification of novel lead candidates with potent dual inhibition of PDE1 and 5 [42]. These compounds have application in the treatment of different cardiovascular diseases, i.e., hypertension, angina, heart failure, and arteriosclerosis [42,43].
They applied virtual screening that consists of classification and regression tree (CART) analysis with the utilization of 168 2-center pharmacophore descriptors and 12 macroscopic descriptors which can result in finding of new lead compounds and drug candidates efficiently. The method applied started with the learning step of CART analysis where a prediction model is configured, and the explanatory variables are selected as per the discrimination index of training set. These variables should be independent from each other.
In the following step, they checked the potential energy term and solvation energy, which are the main constituents of the binding energy used in molecular modeling and simulation, to select the most suitable set of explanatory variables. The potential and solvation energy was obtained by summing up electrostatic potential term, van der Waals potential term, solvation free energy term, and hydrogen bond energy term. This situation is represented qualitatively by pharmacophore where it provides the binding features with the position of these features. In chemoinformatic analysis, the pharmacophore descriptor is either a part or all of pharmacophore. The n-center pharmacophore descriptor is the number of features that participate in binding and the distance between them, for example, 2-center, 3-center, etc. In this research, they adapted 2 for n to reduce the number of explanatory variables. Six binding features were used with eight classes of distance between them, this gives 168 two-center pharmacophore descriptor.
In next stage, the virtual screening has been performed to screen the library of commercially available chemical compounds that was supplied by SPECS Inc. for PDE5 inhibitor activity. Based on Lipinski's rule of five, the compounds with undesirable physicochemical and pharmacokinetic properties were filtered out, followed by selection of compounds with desirable inhibition activities for PDE5 by CART model. The next step involved the selection of structurally diverse 100 compounds. Finally, 19 drug-like compounds out of 100 were selected and obtained from SPECS Inc. and assayed in vitro to test their inhibitory activity against PDE1 and PDE5 [42]. The results showed that the virtual screening along with the CART analysis have a high prediction capability for biological activity of new chemical compounds.

Selective PDE1 inhibitors
One of the least explored PDEs as a drug target singly has been PDE1, although it has been a combined target with PDE5 for vasodilator drugs [11,44]. However, some of the studies have shown that it is a promising target for the treatment of cognitive impairments [45]. One such report by Li et al. recently presented novel PDE1 inhibitors for the treatment of cognitive impairments. They applied both ligand-and structure-based drug design methods to discover these novel PDE1 inhibitors. Achieving high selectivity among PDE enzymes is challenging, since all PDE enzymes share a high degree of sequence homology in the catalytic domains [46][47][48][49][50]. They started by analyzing the published PDE inhibitors with various scaffolds and found out that pyrazolo [3,4-d]pyrimidinones, which was reported by Xia and coworkers, could be the first step for designing a new, potent, and selective PDE1 inhibitors [51]. At the beginning, these compounds were developed as dual inhibitors for both PDE1 and PDE5 [46][47][48][49].
According to the available X-ray crystal structures of human PDE enzymes, a hydrogen bond is formed between the N-7 nitrogen of purine ring in guanine and adjacent amino acid residue in catalytic site. In pyrazolo [3,4-d]pyrimidinone compounds, the hydrogen bond network is disrupted due to the shift of nitrogen in pyrozole ring from position 3 to position 2. Li et al. designed polycyclic 3-aminopyrazolo [3,4-d]pyrimidinone scaffold 2 to reestablish the hydrogen bond through the substitution of amino group at carbon 3 of pyrazole ring [50]. By applying a combination of ligand-based and structure-based drug design methods, they designed numerous novel scaffolds as PDE1 inhibitors. Their work resulted in the discovery of a clinical candidate (ITI-214) which has excellent selectivity toward PDE1 (Figure 4). ITI-214 is now in clinical trial phase I and is being tested for its ability in the treatment of cognitive deficits associated with neurodegenerative and neuropsychiatric disorders and other CNS and non-CNS diseases [45,50].

PDE9A inhibitors
Another PDE which has emerged as a promising drug target recently is PDE9A [52]. PDE9A is now considered as an important therapeutic target for the treatment of diabetes and Alzheimer's disease (AD) [53]. Most of the inhibitors of PDE9A until recently were based on the pyrazolopyrimidinone scaffold, thus there was need to identify novel scaffolds possessing this activity [54]. Zhe Li and team used a combinatorial method including pharmacophores, molecular docking, molecular dynamics simulations, binding free energy calculations, and bioassay to discover novel PDE9A inhibitors with new scaffolds (Figure 3) [55]. SPECS database containing about 200,000 compounds was screened using their combinatorial approach. The combination of ligand-and structure-based methods in combinatorial fashion was done with the aim to reduce computational cost as structure-based methods alone are computationally too expensive for virtual screening. The results were encouraging as has been the case whenever such combinatorial approach has been employed. Fifteen hits out of 29 molecules (a hit rate of 52%) with five novel scaffolds were identified to be PDE9A inhibitors with inhibitory affinities no more than 50 mM to enrich the structural diversity, different from the pyrazolopyrimidinone-derived family. The high hit ratio of 52% for this virtual screening method indicated that the combinatorial method is a good compromise between computational cost and accuracy. Binding pattern analyses indicate that those hits with nonpyrazolopyrimidinone scaffolds can bind the same active site pocket of PDE9A as classical PDE9A inhibitors. The five novel scaffolds discovered in this study can be used for the rational design of PDE9A inhibitors with higher affinities (Figure 5).

PDE2 inhibitors
PDE2 is a key enzyme that hydrolyzes both cAMP and cGMP. It has been suggested that selective PDE2 inhibitors could be a promising therapy for some of the CNS disorders, such as Alzheimer's disease, memory deficit, and depression since PDE2 modulates neuronal signaling Figure 5. Virtual screening strategy to discover novel PDE9 inhibitor scaffolds.
Quantitative Structure-activity Relationship involved in concentration, learning and memory, and emotion [56]. Bo Yang and coworkers examined the binding structures and free energies for PDE2 and benzo [1,4]diazepin-2-one derivatives (PDE2 inhibitors) by combining the molecular docking, molecular dynamics (MD), calculations of binding free energy, and binding energy decompositions [57,58]. Molecular docking was performed followed by energy minimization. The docked structures were analyzed followed by the selection of the best pose for each ligand. The next step was MD stimulation. The binding free energy (ΔG bind ) of PDE2A-ligand complex was determined based on the MD trajectory for each complex [59]. The binding energy decomposition was estimated by using molecular mechanism/generalized Born surface area (MM/GBSA) method.
The results put forth important information regarding the PDE2A-ligand binding patterns including the intermolecular interactions, hydrophobic interactions, and hydrogen bonding.
The estimated PDE2A-inhibitor binding patterns and the agreement between theoretical and experimental outcomes provided a firm base for further design of new, selective, and more potent PDE2A inhibitors [57].
In recent years, several potent inhibitors for PDE2 were developed but none reached the market due to either lack of selectivity or poor pharmacokinetic properties. This led work to focus on the optimization of pharmacokinetic properties of PDE2 inhibitors. Zhang et al. tried to discover selective potent PDE2A inhibitors with improved pharmacokinetic properties [60][61][62][63][64][65][66][67]. In their study, they described the identification of novel PDE2A inhibitors by structure-based virtual screening. They combined pharmacophore model-based screening and molecular docking along with MD simulations and bioassay to find new selective compounds with significant improvement in inhibition activity (Figure 6).
Beginning with the selection of small-molecule database SPECS (comprising almost 200,000 small molecules) for virtual screening, this database was filtered by Lipinski's rule of five to constitute databaset0 which is the initial data set. The crystal structures of PDE2A in a bound  state were used to generate the 3D pharmacophore model. This model was applied for minimizing the size of database through efficient screening of dataset0 to obtain dataset1. Panassay interfering compound substructures (PAINS) screening was used to obtain dataset2 by eliminating the false positive compounds that might interfere with other detecting methods and result in false positives.
In the next step, the dataset2 compounds were submitted for molecular docking to predict the preliminary optimal binding modes and binding energies between PDE2A and ligand. Those compounds with better scores than the reference and proper binding modes constituted dataset3. In the last step, MD stimulations and molecular mechanism/Poison-Boltzman surface area (MM-PBSA) method were applied for more accurate prediction of the binding modes and binding energies. AMBER 10.0 was utilized to stimulate the binding of PDE2A with the compounds in dataset3, whereas MM-PBSA was used for the calculation of binding free energy. At the end, 30 molecules with optimal binding patterns and top binding energies were selected to make up the final dataset. Nine hits out of 30 molecules (a hit rate of 30%) were identified with less than 50 μM affinity for PDE2A.
The result of this study was the discovery of new compound LHB-8 with IC 50 = 570 nM (Figure 4). This compound poses novel scaffold benzo[cd]indol-2(1H)-one among PDE2A inhibitors, which can be used as the novel scaffold for designing of potent inhibitors of PDE2A [67].

PDE3 inhibitors
One of the early studies reporting the application of structure-based drug design on PDE3 inhibitors was done by Fossa et al. [2]. They combined homology modeling techniques with docking to gain knowledge about the molecular requirements of selective inhibition of PDE3 and identified the important amino acids residues for substrate and inhibitor binding. They built a homology model of PDE3A catalytic site based on coupling of PDE4B2B crystal structure and PDE3A mutagenesis data. The amino acids sequences of PDEs were obtained and multiple sequence alignments of PDE isozymes were performed [68][69][70]. Then, phylogenetic relationship was used to calculate the amino acid conservation based on the obtained alignments followed by prediction of secondary structure [71,72]. Fold recognition and sequence to structure alignment was also done during the course [2,73]. 3D models were generated by the application of restraint-based homology modeling methods [74]. At the end, the Amber force field was applied for energy minimization, and the generalized Born solvation model as implemented in AMBER was used for the calculation of molecular mechanic. The structural evaluation of the model proved that it is suitable for docking studies [75].
Another study was proposed by Kim et al. where they used a virtual screening approach to discover novel PDE3 inhibitors for obesity treatment. They started with the analysis of the structural features of the known PDE3 inhibitors and screening of virtual library with 30,000 diverse compounds, followed by docking study based on the 3D structure of PDE3B. In this work, Kim et al. built 3D structures of PDE3B-ligand complex by utilizing a cocrystal structure. Docking was carried out, and 80 compounds with low energy conformation were identified by docking. These compounds were examined for their adipocyte lipolysis activity and finally four structurally unrelated compounds were identified. Among these leads, the most potent compound has IC 50 = 14.8 nM for PDE3A activity and 88.4 nM for PDE3B activity [76].

Future perspective
Apart from those discussed in the previous section, there are many more examples of successful application of in silico methods in the discovery and design of PDE inhibitors. The applications are very diverse ranging from lead discovery to optimization. However, the major problem with various PDE subtype inhibitors, designed or discovered, is the lack of selectivity. Although with the application of computational methods inhibitors with significant selectivity for PDE subtypes have been designed lately, still none of them have been able to make it to the market. The advancement in computational hardware and MD simulation methods will be the major boost for solving the selectivity problem in future. These advancements will allow more deeper studies on a tinier time scale allowing the computational chemist to capture interactions and structural features required for high degree of selectivity. These developments bring hope that the unfulfilled potential of various PDE inhibitors will be realized in near future.