Structure-based discovery and bio-evaluation of a cyclopenta[4,5]thieno[2,3-d]pyrimidin-4-one as a phosphodiesterase 10A inhibitor.

Phosphodiesterase10A (PDE10A) is a potential therapeutic target for the treatment of several neurodegenerative disorders. Thus, extensive efforts of medicinal chemists have been directed toward developing potent PDE10A inhibitors with minimal side effects. However, PDE10A inhibitors are not approved as a treatment for neurodegenerative disorders, possibly due to the lack of research in this area. Therefore, the discovery of novel and diverse scaffolds targeting PDE10A is required. In this study, we described the identification of a new PDE10A inhibitor by structure-based virtual screening combining pharmacophore modelling, molecular docking, molecular dynamics simulations, and biological evaluation. Zinc42657360 with a cyclopenta[4,5]thieno[2,3-d]pyrimidin-4-one scaffold from the zinc database exhibited a significant inhibitory activity of 1.60 μM against PDE10A. The modelling studies demonstrated that Zinc42657360 is involved in three hydrogen bonds with ASN226, THR187 and ASP228, and two aromatic interactions with TYR78 and PHE283, besides the common interactions with the P-clamp residues PHE283 and ILE246. The novel scaffold of Zinc42657360 can be used for the rational design of PDE10A inhibitors with improved affinity. This journal is © The Royal Society of Chemistry.

Introduction

The cyclic nucleotides, 3′,5′-cyclic adenosine monophosphate (cAMP) and 3′,5′-cyclic guanosine monophosphate (cGMP), are well characterised and evolutionarily conserved second messengers that regulate a multitude of cellular functions involved in neuronal signal transduction and synaptic transmission.[1-3] Neuronal regulation of the magnitude and duration of cAMP/cGMP elevation is important to achieve specificity and accuracy of signal transduction. The alterations in cyclic nucleotide availability contribute to changes in neuronal cell functions in the central nervous system (CNS), hence precipitating, maintaining or triggering cognitive, motor or psychiatric disturbances. The availability of cyclic nucleotides is controlled by their rate of synthesis through activation of adenylate and guanylate cyclase and their rate of degradation by phosphodiesterases (PDEs) that hydrolyse cAMP and cGMP to inactive nucleotide monophosphates.[1,4]

PDE10A, is a dual substrate enzyme that catalyses the hydrolysis of cAMP & cGMP. It is expressed in both striatonigral direct (dopamine D1 receptor-expressing) and striatopallidal indirect (dopamine D2 receptor-expressing) pathway medium spiny neurons (MSNs) in the striatum.[5-7] The affinity of PDE10A for cAMP is higher than cGMP by approximately 20-fold.[8] In the direct pathway neurons, PDE10A inhibition activates cAMP/protein kinase A (PKA) signalling leading to the potentiation of D1-receptor signalling, while in the indirect pathway neurons, PDE10A inhibition activates cAMP/PKA signalling by simultaneous potentiation of adenosine A2A receptor signalling and inhibition of D2-receptor signalling.[9] A study of the neuronal type-specific regulation of dopamine- and cAMP-regulated neuronal phosphoprotein (DARPP-32) phosphorylation at Thr34 using neostriatal slices showed that PDE10A inhibitor raised DARPP-32 phosphorylation by 6-fold in the indirect pathway neurons, while it raised DARPP-32 phosphorylation by <2-fold in the direct pathway neurons, thus indicating the predominant effect of PDE10A inhibition in the indirect pathway neurons. Interestingly, the indirect pathway-specific effect of a PDE10A inhibitor is observed with antipsychotic medications, which raise DARPP-32 phosphorylation in dopamine D2-receptor-expressing striatal neurons in mice.[6] This effect is responsible for the improvement of the positive symptoms in schizophrenia.[7] However, the PDE10A inhibition might be associated with extrapyramidal side effects, similar to the side effects observed with the D2-receptor antagonists, which may explain why PDE10A inhibitors have not reached the market yet as antipsychotic treatment.[10]

The PDE10A inhibitors are being clinically investigated as a treatment for Huntington's disease, obsessive-compulsive disorder, schizophrenia and Tourette's syndrome due to the high expression of PDE10A in the indirect pathway neurons.[9,11,12] In addition, it has been suggested that the inhibition of PDE10A provides a novel approach for the treatment of Alzheimer's disease by overcoming the detrimental effects of the amyloid-β peptide on the cAMP-response element-binding protein (CREB) pathway, which is a key control point for long-term memory formation.[13,14] Recently, the extensive efforts of the medicinal chemists have directed toward developing potent PDE10A inhibitors with minimal side effects. However, PDE10A inhibitors are not approved as treatment for neurodegenerative disorders, possibly due to the lack of research in this area. Therefore, the discovery of novel and diverse scaffolds targeting PDE10A is still required. Over the last two decades, the computational methods of drug design and discovery have influenced the overall process of drug discovery.[15] The early application of computational chemistry as a means to understand the molecular basis of PDEs inhibition goes back to the 1980s. The first reported study attempted to explain the relationship between specific physicochemical properties and potency of known inhibitors, and thereby define pharmacophore for PDE inhibitors.[16-18] Since then, a number of PDE10A inhibitor has been discovered and optimised by means of structure-based methods and ligand-based methods.[18] In the present study, a potent inhibitor of PDE10A has been identified using structure-based methods that include virtual screening, molecular docking and molecular dynamics (MD) simulations. The identified inhibitor was able to inhibit the PDE10A activity in the micromolar range in vitro.

Experimental

Computational methods

Selection of appropriate protein structure for the development of pharmacophore model

The first step in the development of a good receptor-based pharmacophore model is the selection of the most appropriate protein structure and active site. Initially, few crystal structures of the target protein, PDE10A, in complex with their co-crystallised ligands were selected from the Protein Data Bank (PDB).[19,20] The selection of the protein structures was based on the following criteria (1) resolution must be high. (2) The active site completeness. (3) The side chain completeness. (4) And the active site must contain a co-crystallised ligand. Secondly, one of the selected complexes was assigned as the reference complex where all the other protein complexes were superimposed on it. Finally, the active site of the reference protein was examined to determine whether the superimposed ligands can fit completely into the selected site of the reference protein without any steric clashes that might affect the binding. Maestro (Release 2019-4, Schrödinger software) and Discovery Studio (DS) Visualizer were used for performing the binding site alignment.[21,22]

Receptor-based pharmacophore modelling

Based on the previous step, one PDE10A in complex with its co-crystallised ligand was chosen to develop the receptor-based or structure-based pharmacophore model by finding the pharmacophoric features in the active site that are important for the binding of the ligand. LigandScout 4.4 has been used for the generation of a pharmacophore model that comprises all the important chemical features.[23] The interactions between the target and the co-crystallised ligand were used to identify the aromatic ring (AR), hydrophobic area (HA), hydrogen bond acceptor (HBA), hydrogen bond donor (HBD) and ionisable area features. These features represent the pharmacophoric features of the receptor-based pharmacophore.

Receiver operating characteristic curve

The receiver operating characteristic (ROC) curve has been used to evaluate the capability of the developed pharmacophore model in discriminating the active compounds and decoys. The ROC curve is generated by plotting the Sensitivity against 1-Specificity, where the sensitivity indicates the inclusion of the true positives, whereas; specificity indicates the exclusion of the true negatives from the validation result.[24] The ROC lies between 0 and 1, where a value near one is considered favourable. From the reported PDE10A inhibitors in the literature, 11 inhibitors were selected to constitute the test set whereas, 1000 decoys (without reported PDE10A inhibitory activity) were obtained from PubChem. The developed pharmacophore model was validated by screening against the training set and decoy using LigandScout 4.4. The reliability of the pharmacophore model was verified by calculating the ROC value, in addition to other statistical values like the goodness of hit (GH) and enrichment factor (EF).[25]

Multi-step virtual screening

Pharmacophore-based screening

More than six million lead-like compounds were downloaded from the zinc database to be used in the virtual screening. The first step of the screening was the preparation of the database by converting the lead-like chemical compounds into the IDB file format. This step was performed because LigandScout uses the IDB format for storing the compound libraries. Then, the validated pharmacophore model of PDE10A was used as a query for the screening of the database and retrieving hits that match the corresponding pharmacophoric features. The retained hits were further subjected to drug-likeness analysis using Lipinski's rule of five and other filters, i.e. polar surface area <60 Å, blood/brain partition coefficient 0–1.2, and CNS active drug 0–2. Lipinski's rule of five has been modified according to the CNS active drug's criteria. In which, compounds with molecular weight <450 Da, log P 2–5, hydrogen bond donor <3, hydrogen bond acceptor <10, and rotatable bonds <10, were retained. Then, a pharmacophore fit score (FPS) was applied to retrieve compounds with FPS ≥48 (the selected value based on the FPS of the training set compounds). Finally, the hit compounds that passed all the screening filters were submitted to the PAINS-Remover server to remove the Pan Assay Interference Compounds (PAINS) from the retrieved hits.[26]

Molecular docking

Protein preparation

The crystal structure of PDE10A in complex with its co-crystallised ligand was prepared by the Protein Preparation Wizard module in Maestro (Schrödinger software).[21,27] The protein structure was refined by assigning the bond orders, adding the missing hydrogens, removing the water molecules beyond 5 Å from the co-crystallised ligand and water molecules with less than 3 hydrogen bonds to the co-crystallised ligand. This followed by optimising the protonation state using PROPKA at pH: 7.0 and finally minimising the protein structure to relieve the stearic clashes using optimised potentials for liquid simulations (OPLS_2005) force field.

Molecular docking using AutoDock Vina

AutoDock Vina 1.1.2 was used to perform molecular docking. The program uses a sophisticated gradient optimisation algorithm and a hybrid scoring function (empirical + knowledge-based) based on the X-Score function.[28,29] At first, the prepared PDE10A structure was imported into AutoDock 4.0, where the graphical user interface program AutoDockTools (ADT) 1.5.6 was used to convert and save the protein as pdbqt format. Then, all the ligands obtained from the virtual screening were imported separately into AutoDock 4.0 where ADT was used to prepare the ligands for docking by setting the rotatable bonds and allow all the torsions to rotate for the ligands. Finally, the ligands were saved as pdbqt format. The PDE10A inhibitor, TAK063, was used as a standard for comparing the docking results. The location of the grid box was assigned based on the coordinate information of the active site and the atoms of the co-crystallised ligand; thus, the grid box encompasses the active site completely. The configuration file was prepared, which; contain information about the target receptor, the ligand, the dimensions (X: −20.304, Y: 12.903, Z: −35.36) and size (40 × 40 × 40) of the grid box and where to send the docking results. The file was saved as conf.txt. By preparing the protein, ligand, grid box and conf. file, molecular docking was carried out using the command “vina.exe --config conf.txt” in AutoDock Vina 1.1.2.[28] The analysis of the docking results and the protein–ligand interactions was performed using DS Visualizer.[22]

Molecular dynamics simulations

The docked protein–ligand complex with the lowest binding energy was subjected to molecular dynamics (MD) simulations using PMEMD.CUDA in AMBER 14 package.[30] The general AMBER force field (GAFF) parameters were assigned for the ligand, and the restricted electrostatic potential (RESP) was calculated under Gaussian g09. The parameters and topologies for the protein–ligand complex were obtained by tleap, and Amber forcefield (ff14SB) was assigned to the complex. The complex system was solvated in a cubic box of TIP3P water with a margin distance of 10 Å in each direction and was neutralised by the addition of Na+ ions. All simulations were carried out under periodic boundary conditions, and the long-range electrostatics were treated using the particle mesh Ewald (PME) method. MD simulations were carried out in four steps (1) relaxation and energy minimisation with a series of steepest descent and conjugated gradient up to 1000 cycles to remove the steric hindrance in the model and use the initial local minimum structure prior running MD simulation; (2) heating the system to 310.15 K using Langevin thermostat at constant volume; (3) equilibration by fixing and equilibrating the protein atoms for 1000 ps at a constant pressure of 1 atm; and (4) finally, a production simulation run for 150 ns using the NPT ensemble under a constant temperature of 310.15 K and pressure of 1 atm. The time step was set to 2 fs while the trajectory was recorded every 1 ps.

Molecular dynamics trajectories analysis

The generated trajectories from MD simulations were analysed through the cpptraj module of AmberTools 14. The root mean square deviation (RMSD) for the backbone atoms of the complex relative to the starting structure was analysed to quantify the conformational changes in the protein structure. The root mean square fluctuation (RMSF) of the amino acid residues was analysed to determine the flexibility of the residues, especially the residues located in the binding pocket. Additionally, the hydrogen bond analysis was performed to evaluate the strength of the hydrogen bonds formed between the protein and the ligand.

Binding free energy calculation and residue decomposition using MM-GBSA

The binding free energy of the complexes was calculated from the stable MD trajectories of the last 20 ns using, the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) method implemented in AMBER 14. In this method, the binding free energy, was calculated as the difference between the free energy of the bound and unbound state of the protein and the ligand.ΔΔΔWhere ΔGbind is the total binding free energy upon protein–ligand binding; Gcomple, Gprotein and Gligand are the free energies of the complex, protein and ligand, respectively. ΔEMM (gas-phase energy) is the change of the molecular mechanics (MM) potential energy upon binding which involves van der Waals (VDW), electrostatic (ELE) and internal energies (INT). The INT results from angle, bond and dihedral terms in the MM force field (this term always amounts to zero in the MM/GBSA calculations based on the single trajectory of a complex). ΔGsolv is the solvation free energy which is the dissolution free energy penalty upon the ligand binding. It is calculated from the polar and nonpolar components of the solvation free energy. −TΔS is the change of the entropy upon binding which was not considered in this study since we are comparing states of similar entropy i.e. ligands binding to the same protein.[31-35] Additionally, to obtain information on the main residues involved in the binding of the protein–ligand complex, the free energy of binding was decomposed to the energy contributed by a single residue. It is calculated by integrating the interactions of each residue over all residues in the system. Graphical representation was generated by OriginPro 9.1 and Excel 2016 software.

Biological assay

Phosphodiesterase10A inhibition assay

To evaluate the potency of the inhibitor with the highest affinity for PDE10A, the compound was tested in vitro to determine the PDE10A-IC50 value. The compounds, Zinc42657360, TAK063 and 3-isobutyl-1-methylxanthine (IBMX), were procured from Mcule, Inc. (Palo Alto, CA, USA). The compounds' purity was determined to be ≥90% by liquid chromatography-mass spectrometry (LCMS) and nuclear magnetic resonance (NMR). The PDE inhibitory activity of the tested compound was evaluated using the PDE-Glo phosphodiesterase assay according to the manufacturer's portocol (Promega Corporation, Madison, WI, USA). 7.5 μl of 1× PDE-Glo Reaction Buffer containing 16.5 ng of recombinant human PDE10A (SignalChem, Richmond, Canada) was pipetted into 96-well plate wells (SPL Life sciences, Gyeonggi-do, Korea). The tested compound was dissolved in Dimethyl Sulfoxide, and a serial dilution of the inhibitor was performed using 1× PDE-Glo Reaction Buffer. 5 μl of the diluted inhibitor and 12.5 μl of the cAMP solution was added to each well. After 12 min incubation at 30 °C, 12.5 μl of PDE-Glo Termination Buffer and 12.5 μl of PDE-Glo Detection Solution were added, and the plate was further incubated for 22 min at 25 °C. Finally, 50 μl of Kinase-Glo Reagent was pipetted into each well and after 10 min of incubation at 25 °C, the luminescence was measured using a microplate reader (FLUOstar Omega, BMG LABTECH, Ortenberg, Germany). The data analysis was performed with GraphPad Prism 8.4.3 software (GraphPad Software Inc.) using a sigmoidal dose-normalised response (variable slope) equation. All data points were the average of two replicates.

Results and discussion

The structure of PDE10A (PDB ID: 5UWF) in complex with the inhibitor (16d); was selected for the development of the pharmacophore model based on the binding site alignment.[36] In which, the alignment of four protein complexes (PDB ID: 3HQW, 4DDL, 4P0N, and 4HF4) on 5UWF showed complete fitness of the four ligands into the active site of 5UWF (Fig. 1, Table 1).[37-40]

Fig. 1

Location of the binding site of PDE10A. (A) The four proteins, 3HQW (green), 4DDL (purple), 4HF4 (yellow) and 4P0N (grey) are superimposed on the reference protein 5UWF (red). (B) Position of the five co-crystallised ligands in the binding site of 5UWF. 5UWF co-crystallised ligand (grey), 3HQW co-crystallised ligand (green), 4DDL co-crystallised ligand (purple), 4HF4 co-crystallised ligand (yellow) and 4P0N co-crystallised ligand (red).

The structural parameters of the selected complexes of PDE10A

PDB ID	Active site's residues	Resolution (Å)	Co-crystallised ligand
5UWF	TYR524, HIS525, LEU675, VAL678, ILE692, TYR693, PHE696, ILE711, MET713, MET714, GLN726, PHE729	1.8
3HQW	VAL668, ILE682, TRY683, PHE686, PRO702, MET703, LYS708, GLU711, VAL712, GLY715, GLN716, PHE719, ALA722, VAL723	1.7
4DDL	SER563, LEU625, LEU665, VAL668, ILE682, PHE686, ILE701, MET703, GLN716, PHE719	2.0
4P0N	SER563, LEU625, LEU665, VAL668, ILE682, TYR683, PHE686, PRO702, MET703, GLU711, GLY715, GLN716	2.0
4HF4	TYR514, HIS515, LEU665, SER667, VAL668, ILE682, TYR683, PHE686, PRO702, MAT703, LYS708, GLU711, VAL712, GLY715, GLN716, PHE719	2.0

As a result, the receptor-based pharmacophore model was developed based on the hydrophilic and lipophilic interaction points available in the active site of 5UWF. The crucial interactions for binding of 16d to 5UWF were considered during the development of the pharmacophore model. Consequently, the final model comprised six pharmacophoric features; one HBA oriented towards GLN726, one AR represented by the thiophene, and four HA represented by the phenol, sulphide and the fluorine groups (Fig. 2).

Fig. 2

Receptor-based pharmacophore model. (A) The interactions between the amino acid residues of 5UWF and the pharmacophore groups of its co-crystallised ligand 16d. (B) Pharmacophoric features of PDE10A pharmacophore model.

Then, the capability of the receptor-based pharmacophore model in classifying the compounds correctly as active or inactive was evaluated using the ROC curve (Fig. 3). Based on the results, the pharmacophore model was able to identify all the active compounds in the dataset and predict nine inactive to be active. The sensitivity of the pharmacophore was found to be one, and the specificity was 0.99. The obtained ROC, GH and EF values were 1, 0.65, and 50.55, respectively. All the calculated values pointed to the appropriateness of the pharmacophore model to be used as a query in the virtual screening to identify PDE10A inhibitors.

Fig. 3

ROC curve for PDE10A-pharmacophor.

Multi-step virtual screening was conducted in six phases (1) pharmacophore-based screening; (2) drug-likeness analysis; (3) pharmacophore fit score; (4) PAINS filter; (5) molecular docking, and; (6) molecular dynamics simulations. The schematic workflow of the multi-step virtual screening protocol employed in this study is shown in Fig. 4. Pharmacophore-based screening of zinc database was performed using the validated receptor-based pharmacophore model as a 3D query to identify potential PDE10A inhibitors with matched pharmacophoric features. As a result, a total of 7541 hits were retrieved in the first phase. Then, several filters were applied to retain only compounds with drug-likeness properties. Consequently, 516 hits were retrieved, in which the physicochemical properties of those hits fit the criteria of the CNS active drug. In this criteria, the compound must be relatively small, moderately lipophilic, has a small polar surface area <60 Å, blood/brain partition coefficient in the range of 0–1.2 and the value of CNS in the range of 0–2, to be considered as a CNS active drug. Compounds that fit the mentioned criteria are most probably able to cross the BBB and have CNS activity. Following the drug-likeness analysis, the 516 hits were subjected to the pharmacophore fit score filter. In this filter, compounds with fit values ≥48 to the pharmacophoric features were retained. As a consequence, only 14 screened hits were obtained in this step (Table S1†). Subsequently, the 14 hits were subjected to the fourth filter in which they were submitted to the PAINS-Remover server to remove the false positive compounds. This filter screens the compounds based on a number of sub-structural features that assist in the identification of problematic screening compounds, which appear as frequent hitters (promiscuous compounds) in many screening libraries. The compounds can be considered false positive for several reasons, i.e. interference in the binding interactions by the formation of aggregates; interference in the assay signalling; and protein-reactive compounds.[26] Accordingly, all the 14 hits passed this filter and were considered in docking studies.

Fig. 4

Schematic workflow of the multi-step virtual screening protocol employed in the identification of PDE10A inhibitor.

Molecular docking studies are performed to predict the affinity of the ligand to the target receptor, in addition to the preferred pose and conformation in the complex. In this study, AutoDock Vina 1.1.2 has been used to carry out the molecular docking of 14 compounds. This step was considered the fifth filter in the screening process, where all the compounds were filtered based on their affinities to PDE10A. Before molecular docking, the protocol was validated by re-docking of the co-crystallised ligand, JY4, into the binding site of PDE10A. This step was performed to evaluate the ability of the docking protocol to produce the bioactive conformation. As a consequence, the docked pose with the lowest binding energy adapted a binding mode similar to that of the co-crystallised ligand. Among the docked ligands, the standard compound, TAK063, and Zinc42657360 showed the lowest binding energy of −9.1 kcal mol−1 and −9.2 kcal mol−1, respectively (Tables 2 and S2†). Therefore, the underlying binding interactions of these two compounds with PDE10A were further analysed to explore the structural features that contribute to the PDE10A affinity. The architecture of the binding site of PDE10A achieves all the required criteria for a drug-able binding site. It consists of three pockets, (1) M pocket represents the area in the enzyme that contains the important metal ions (Zn+2 and Mg+2) for catalysing the cAMP and cGMP hydrolysis. (2) S pocket, signifies the solvent filled side pocket which has polar residues and, (3) Q pocket which is known as an inhibitor pocket and further divided into hydrophobic clamp (P-clamp) and conserved purine-selective glutamine. The P-clamp is a rigid and small <300 Å hydrophobic cavity that is ∼14 Å width, ∼13 Å depth and ∼8 Å height. It contains a conserved aromatic PHE located at the roof of the P-clamp; and two hydrophobic residues, ILE and PHE, located on the floor of the binding site. Whereas the purine-selective glutamine pocket contains an invariant substrate-recognising GLN residue which is critical for substrate or ligand recognition. Several studies have shown that the PDE10A inhibitors share two features, a planar ring structure held within the hydrophobic residues of the P-clamp (PHE283 at the roof of the P-clamp and ILE246 and PHE250 on the floor of the binding site); and a hydrogen bonding interaction with the invariant glutamine residue (GLN280). The cyclic nucleotides are recognised by the enzyme upon the formation of two hydrogen bonds for cAMP and one hydrogen bond for cGMP with the Q pocket. Additionally, the phosphate moiety forms a complex with the metal in the M pocket, promoting the hydrolysis of cAMP/cAMP. Thus, the PDE10A inhibitor should be able to occupy either the Q/M pocket or both to block the entry and hydrolysis of cAMP/cGMP.[41,42]

Molecular docking binding affinities to PDE10A

Ligand	PDE10A binding energy (kcal mole−1) PDB ID: 6MSA	2D-structure
Co-crystallised ligand (JY4)	−8.0
TAK063 (standard)	−9.1
Zinc01397213	−7.1
Zinc02156284	−6.1
Zinc42657360	−9.2
Zinc43638301	−6.5
Zinc47464611	−8.8
Zinc71439134	−6.3
Zinc71759377	−7.6
Zinc72553806	−8.3
Zinc72878277	−7.2
Zinc79055898	−7.6
Zinc82446000	−5.9
Zinc82779572	−7.2
Zinc82779574	−7.9
Zinc82779590	−7.7

The docking study has revealed that both TAK063, and Zinc42657360 occupied the P-clamp of PDE10A and interacted with the main residues, ILE246, PHE250 and PHE283. In regard to Zinc42657360, the cyclopentathiophene moiety interacts with the two residues that form the P-clamp, ILE246 and PHE283. Whereas, for TAK063, the two pyrazole groups interact with PHE250 and PHE283, respectively, and the fluorobenzene group interacts with ILE246 and pHE250 (Fig. 5). Aromatic interactions have seen between the cyclopentathiophene and the pyrimidine groups of Zinc42657360 and the residues PHE283 and TYR78, respectively; and between the pyridazine group of TAK063 and the aromatic residue HID79. Further, three hydrogen bond interactions were observed only in the PDE10A-Zinc42657360 complex. The hydrogen bonds were formed between the hydroxyl group of the ligand and the residues THR187 and ASP228 and between the fluorine group and the residue ASN126 (Fig. 6). Another compound has shown a high affinity to PDE10A in comparison to the rest of the ligands. This compound is Zinc47464611 with a binding energy of −8.8 kcal mol−1. The interactions observed between PDE10A and Zinc47464611 were similar to that observed between PDE10A and the two ligands, Zinc42657360 and TAK063. In which, Zinc47464611 also occupy the P-clamp and interact with ILE246, PHE250 and PHE283. Moreover, the trifluoromethyl benzene group displays aromatic interaction with HID79 (Fig. 7).

Fig. 5

The binding interactions of PDE10A with TAK063. Hydrophobic interactions are presented by pink and Pi–Pi stacking interactions by magenta dotted lines.

Fig. 6

The binding interactions of PDE10A with Zinc42657360. Hydrogen bond interactions are presented by green, hydrophobic interactions by pink, Pi–Pi stacking interactions by magenta dotted lines.

Fig. 7

The binding interactions of PDE10A with Zinc47464611. Hydrophobic interactions are presented by pink and Pi–Pi stacking interactions by magenta dotted lines.

According to the binding mode of the three ligands within the active site of PDE10A, we observe that they all formed aromatic interactions or Pi–Pi stacked interactions due to the presence of one or more aromatic rings in their structures. In addition, the three ligands occupied the P-clamp and interacted with the main residues within the clamp. This occupation plays an important role in binding and stabilising the ligands in the active site of PDE10A. However, Zinc42657360 has shown the highest affinity to PDE10A among the three ligands, which could be attributed to the three hydrogen bond interactions formed between the ligand and the active site residues. Thus, the aromatic interactions and the hydrophobic interactions, as well as the hydrogen bonding, showed that TAK-063, Zinc42657360 and Zinc47464611 have an orientation within the active site that would inhibit the entry of cAMP/cGMP into the catalytic domain. These results are consistent with the previous studies that reported the structural elements and key interactions of PDE10A inhibitors like papaverine and MP-10 to achieve high potency and selectivity. These interactions included hydrophobic clamp occupation and aromatic interactions in addition to hydrogen bond interactions with the residues in Q pocket.[37,42]

The three compounds bound similarly to PDE10A as illustrated in Fig. 5–7. However, the difference in the values of the binding energy for the three compounds cannot be illustrated by molecular docking alone. This could be attributed to the inaccuracy in the calculation of the binding energies of protein–ligand complexes using molecular mechanics-based force field, and the induced-fit effects which are not included in the docking protocol. Therefore, four independent MD simulations were performed for unbound-PDE10A and PDE10A in complex with TAK063, Zinc42657360 and Zinc47464611, respectively to predict the dynamic binding behaviour between the protein and the three ligands in the aqueous environment. The MD results of the various simulations are presented and analysed comparatively in order to better characterise the structure and flexibility of the protein and the ligand's mode of binding. At first, UCSF Chimera 1.13.1 was used for visualising the MD trajectories of each protein–ligand complex that produced in the production stages before proceeding with the analysis of the results; to ensure all the ligands present in the binding site of the target receptor for the entire simulation time. As a consequence, Zinc47464611 was found unbound to its target receptor and away from the binding site at 127 ns (Fig. 8), which; indicates the instability of the ligand within the active site of PDE10A. Therefore, Zinc47464611 was excluded from the rest of the analysis.

Fig. 8

PDE10A-Zinc47464611 complex. (A) The ligand presents within the binding site at 126 ns. (B) The ligand is unbound to the receptor at 127 ns.

In order to analyse the changes in the relative position of the backbone atoms of the two complexes and validate the stability of the complexes, the time evolution of the RMSD of the backbone atoms' coordinates was monitored during the 150 ns with respect to the initial crystallographic structures. As seen in Fig. 9, the RMSD values were oscillating steadily up to 150 ns simulation for the unbound-PDE10A and the two PDE10A-complexes, which indicate the minor conformational changes in the structure of the protein upon ligands' binding. In regard to the unbound-PDE10A, the RMSD value showed a slight increment to 2.79 Å at 35 ns; in PDE10A-standard (TAK063) complex, the RMSD value increased to 3.2 Å at 48 ns; and in PDE10A-Zinc42657360 complex, the RMSD value increased to 3.0 Å at 75 ns. However, the two ligands remained bound within the binding site of PDE10A, and their RMSD values remained constant throughout the entire simulation time at a range of 1.6–3.2 Å.

Fig. 9

Time evolution along 150 ns MD simulation trajectories of the root mean squared deviation computed for the backbone atoms. Black: Unbound-PDE10A, Grey: PDE10A-TAK063, Green: PDE10A-Zinc42657360.

In general, more fluctuations were observed in PDE10A-Zinc42657360 complex than that in the unbound-PDE10A and PDE10A-TAK063 complex (Fig. 10). The values of the RMSF, which reflect the individual residues flexibility during the MD simulations, showed minimum fluctuation in the three systems; whereas the RMSF values of the residues HID79 (2.0 Å, 0.7 Å and 0.9 Å, respectively), ASN80 (2.0 Å, 0.7 Å and 0.9 Å, respectively), ASP228 (1.5 Å, 0.4 Å and 0.4 Å, respectively), ASN244 (3.0 Å, 1.1 Å and 1.6 Å, respectively), ARG270 (3.1 Å, 1.6 Å and 2.1, Å respectively), ASP271 (3.5 Å, 2.03 Å and 2.4 Å, respectively) and LYS272 (3.3 Å, 1.8 Å and 2.3 Å, respectively) varied significantly. The highest fluctuation in the unbound-PDE10A was observed in the residues GLY261 and GLY322 with RMSF score of 2.0 Å. While, in PDE10A-Zinc42657360 complex, the highest fluctuations were observed in GLY261(RMFS score 3.5 Å) and ASP271 (RMSF score 3.5 Å); and in PDE10A-TAK063 complex, the highest fluctuations were observed in THR205 (RMSF score 2.9 Å) and GLY206 (RMSF score 2.8 Å). These results indicated the rigidity of the residues within the binding pocket once the inhibitor bound to PDE10A, in which binding of these two inhibitors decreases the flexibility of these residues (TYR78, LEU189, LEU229, ILE246, PHE250, MET267, GLN280, PHE283) hence, making the complexes more stable.

Fig. 10

Root mean square fluctuation of the residues in the PDE10A–ligand complexes. Black: Unbound-PDE10A, Grey: PDE10A-TAK063, Green: PDE10A-Zinc42657360.

The number of hydrogen bonds and occupancies were also calculated to determine the strength of the intermolecular hydrogen bond interactions between the receptor and the ligands. The hydrogen bond is a direct interaction between the donor and acceptor where the interaction strength depends on the distance and electronegativity of the acceptor. By comparing the hydrogen bond interactions between the proteins and the ligands obtained by molecular docking and MD simulations, we notice that the results were inconsistent for both TAK063 and Zinc42657360. This might be attributed to the rigidity of the docking process, where the protein is treated as a rigid structure and the ligand as a flexible structure. Whereas in MD, both protein and ligand are treated as flexible structures. Thus, the ligand might not form a hydrogen bond with a specific conformation in molecular docking while forming hydrogen bonds with the different conformations of the same protein in MD simulations. The MD results have shown that Zinc42657360 was involved in four hydrogen bond interactions with PDE10A. These hydrogen bonds were formed between THR187 and the hydroxyl group of the ligand with 72.6% occupancy, in which THR187 acted as hydrogen bond acceptor and between ASN126 and fluorine group of the ligand with 5.6% occupancy, where ASN126 acted as hydrogen bond donor. These hydrogen bonds were displayed by the molecular docking results as well. Additionally, LEU189 formed a hydrogen bond with the oxygen group of the ligand with 34.1% occupancy, and ASN126 formed another hydrogen bond with the second fluorine group with 0.1% occupancy. The hydrogen bond between ASN126 and fluorine of the ligand was very week that it stood for a very short time during the MD simulation. In regard to PDE10A-TAK063, the two hydrogen bonds were formed between the residues GLN280 and LEU189 and the amine groups of the ligand with occupancies of 23.7% and 0.6%, respectively. Both residues were acted as hydrogen bond donors (Fig. 11 and 12).

Fig. 11

Hydrogen bond occupancies between PDE10A and TAK063 during MD simulations.

Fig. 12

Hydrogen bond occupancies between PDE10A and Zinc42657360 during MD simulations.

The calculations of the binding free energy for PDE10A complexes have been performed using the MM-GBSA method, which is usually utilised to predict the binding affinity that could assist in identifying the compound with inhibitory potential. The calculation of the free energy difference between two states, bound and unbound, of protein and ligand enables us to estimate the binding affinity by calculating the average results of interaction energies (ΔGgas) and the solvation free energy (ΔGsolv) of the ligand, protein and complex in which these energies contributed favourably to the protein–ligand binding. Both complexes displayed high energy of binding indicated by ΔGtotal value of −21.6 and −24.0 kcal mol−1 for PDE10A-TAK063 complex and PDE10A-Zinc42657360 complex, respectively. The major contributors to the binding free energy were VDW and ELE energies as calculated using molecular mechanics (MM) force field. The VDW interaction is considered the favouring ligand binding energy where large ligands with more atoms are prone to have higher VDW interactions than small ligands. Whereas, the ELE interaction is an important force in the primary approach of the ligand and receptor to each other. These types of interactions are crucial in the stability of the protein–ligand complex. Moreover, the nonpolar component of the solvation energy (ESURF) were almost similar in both complexes and favourable for binding as well. However, they have less contribution to the binding energy due to the relatively small negative values. In contrast, the polar solvation free energy (EGB) was shown to be unfavourable for binding, as indicated by the positive values (Table 3).

Binding free energy and its components' values for PDE10A–ligand complexes from two MD simulationsa

Method	Contribution	PDE10-TAK063	PDE10A-Zinc42657360
MM (kcal mol−1)	VDW	−46.7 ± 2.0	−34.4 ± 3.4
	EEL	−16.4 ± 3.1	−46.7 ± 5.0
	ΔGgas	−63.2 ± 3.6	−81.2 ± 5.2
GBSA (kcal mol−1)	EGB	45.9 ± 3.1	60.5 ± 3.8
	ESURF	−4.3 ± 0.1	−3.3 ± 0.1
	ΔGsolv	41.5 ± 3.0	57.1 ± 3.8
	ΔG total	−21.6 ± 2.4	−24.0 ± 3.7

*MM: molecular mechanics energies, VDW: van der Waals interactions, EEL: electrostatic interactions, ESURF: non-polar contribution to solvation, EGB: polar contribution of solvation, ΔGtotal: total binding free energy.

In order to identify the important amino acid residues in PDE10A that interact with the ligands, the analysis of the per-residue free energy decomposition was performed to identify the fundamental basis of the protein–ligand interaction. The energy contribution of a single residue in the binding of the receptor with the ligand is divided into three parts VDW, ELE and solvation. Each residue exhibits a negative or positive influence on the receptor–ligand binding. Accordingly, the residues are considered to be effective contributors in the stabilisation energy if their relative energies are <−1 kcal mol−1. The decomposition of energy shows that the per-residue interaction energy varies in the range 0.02 to −3.1 kcal mol−1 for PDE10A complexes. LEU189 contributed significantly to the binding in both complexes in which the decomposition energy was <−3 kcal mol−1. Further, TAK063 bound strongly to PDE10A through residues ILE246, PHE250, GLN280 and PHE283 with favourable energy <−1 kcal mol−1 (Fig. 13). In contrast, the energy contributed by GLN280 in the PDE10A-Zinc42657360 complex was 0.01 kcal mol−1. This energy was unfavourable for binding; hence it could play the opposite role in the protein–ligand binding and result in weakening the stability of the complex.

Fig. 13

Binding free energy decomposition of the significant amino acids residues of PDE10A complexes. Yellow: PDE10A-TAK063, Grey: PDE10A-Zinc42657360.

Phosphodiesterase10A inhibition assay

The effect of Zinc42657360 on the PDE10A enzyme activity was investigated using PDE-Glo phosphodiesterase assay. TAK063 and IBMX were used as standard inhibitors. At first, the inhibitory activity for the three compounds was tested at a fixed concentration of 10 μM. Both TAK063 and Zinc42657360 inhibited PDE10A significantly at this concentration, in contrast to IBMX, which did not show any inhibitory activity against PDE10A. Thus, a serial dilution of different concentrations was prepared for TAK063 and Zinc42657360 and screened against PDE10A to identify the concentration that inhibits 50% of enzyme activity. As a result, the obtained IC50 value for TAK063 was 0.34 nM, which is close to the reported value in the literature (0.30 nM). Whereas; the obtained IC50 value for Zinc42657360 was 1.60 μM (Fig. 14).[43]

Fig. 14

PDE10A activity (%) in the presence of different concentrations of the tested inhibitors.

A number of studies have emphasised the strong relationship between scaffold type and activity toward PDE10A, in which scaffolds possessing less nitrogen (N) have a lower PDE10A inhibitor potency.[44,45] Zinc42657360 has two nitrogen groups which are part of the pyrimidine ring. This ring is important for inhibiting the activity of PDE10A by interacting with specific residues in the active site, as shown in the molecular docking study. Accordingly, Zinc42657360 can be used as a starting point for the identification of selective and potent PDE10A inhibitor by optimising the structure of the compound and designing potent analogues.

Conclusion

An integrated virtual screening campaign was used to discover a new PDE10A inhibitor for the treatment of several neurodegenerative disorders through pharmacophore model screening integrated with comprehensive molecular docking. The MD simulations were performed to probe the binding mode and stability of the predicted protein–ligand complexes. Fourteen compounds from the zinc database were subjected to molecular docking in which two compounds displayed the highest affinity to PDE10A, thereby selected for MD simulations. Out of the two compounds, Zinc42657360 was subjected to PDE-Glo phosphodiesterase assay, where it exhibited a significant inhibitory activity of 1.60 μM against PDE10A. The above findings suggested that Zinc42657360 can be explored for further development of new inhibitors of PDE10A. It can be used as a starting point for the discovery of a new series of potent analogues with improved affinity. The strategy used in this study makes a good concession between computational cost and improvement of the hit ratio of drug discovery, thus, may have a broad application in further rational drug design.

Author contributions

Mayasah Al-Nema: conceptualization, methodology, investigation, writing-original draft. Anand Gaurav: conceptualization, methodology, validation, writing – review & editing, supervision, project administration, funding acquisition. Vannajan Sanghiran Lee: writing – review & editing, supervision. Baskaran Gunasekaran: writing – review & editing, supervision. Ming Tatt Lee: writing – review & editing, supervision. Patrick Okechukwu: writing – review & editing, supervision. Piyarat Nimmanpipug: editing & funding.

Conflicts of interest

There are no conflicts to declare.