An insight from computational approach to explore novel, high-affinity phosphodiesterase 10A inhibitors for neurological disorders

Abstract The enzyme Phosphodiesterase 10A (PDE10A) plays a regulatory role in the cAMP/protein kinase A (PKA) signaling pathway by means of hydrolyzing cAMP and cGMP. PDE10A emerges as a relevant pharmacological drug target for neurological conditions such as psychosis, schizophrenia, Parkinson's, Huntington’s disease, and other memory-related disorders. In the current study, we subjected a set of 1,2,3-triazoles to be explored as PDE10A inhibitors using diverse computational approaches, including molecular docking, classical molecular dynamics (MD) simulations, Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculations, steered MD, and umbrella sampling simulations. Molecular docking of cocrystallized ligands papaverine and PFJ, along with a set of in-house synthesized molecules, suggested that molecule 3i haded the highest binding affinity, followed by 3h and 3j. Furthermore, the structural stability studies using MD and MM-PBSA indicated that the 3h and 3j formed stable complexes with PDE10A. The binding free energy of −240.642 kJ/mol and −201.406 kJ/mol was observed for 3h and 3j, respectively. However, the cocrystallized ligands papaverine and PFJ exhibited comparitively higher binding free energy values of −202.030 kJ/mol and −138.764 kJ/mol, respectively. Additionally, steered MD and umbrella sampling simulations provided conclusive evidence that the molecules 3h and 3j could be exploited as promising candidates to target PDE10A. Communicated by Ramaswamy H. Sarma


Introduction
Phosphodiesterases (PDEs) are a group of highly conserved isoenzymes regarded as hydrolases and plays a key role in regulating the intracellular level of second messengers: cyclic guanosine monophosphate (cGMP) and cyclic adenosine monophosphate (cAMP) (Amin et al., 2021), and the role of these messengers is well reported in several physiological processes like neuronal signaling, cell-cycle regulation, muscle relaxation, visual transduction, cell proliferation and differentiation, inflammation, gene expression, apoptosis, energy homeostasis and cell growth (Boswell-Smith et al., 2006;Peng et al., 2018).PDEs are involved in diverse biological functioning in the human body and could be considered as a pharmacologically relevant target against related disorders.
The human genome has 11 families of PDE genes (PDE1-11) (Xi et al., 2022), comprising 21 genes that generate approximately 100 proteins due to different transcription start sites, multiple promoters, and alternate splicing (Francis et al., 2011).The PDE families are highly regulated and structurally related but differ in their three-dimensional structures, intracellular localization, cell expression, catalytic properties, and responses to specific modulators and inhibitors (Jansen et al., 2016).All the PDE families based on substrate specificity are divided into three groups.One group hydrolyzes cAMP constitutes PDE4, PDE7, and PDE8, whereas the other group consist of PDE5, PDE6, and PDE9, which hydrolyzes cGMP.The third group members include PDE1, PDE2, PDE3, PDE10, and PDE11, possess dual-specificity, having the capability to hydrolyze both cAMP and cGMP with varying affinities (Conti & Beavo, 2007).
PDE10A is a dual substrate enzyme capable of accomodating and hydrolyzing both cAMP and cGMP with Km values of 56 nmol/L and 4.4 mmol/L, respectively (Yang et al., 2020).The expression of PDE10A is not restricted up to the central nervous system, and it is also expressed in other organs like the kidneys, testis, lungs, thyroid, colon, ovary, and thyroid.A higher PDE10A expression was observed in adipose tissues, brain, colon, and kidneys, while moderate in the heart, hippocampus, thalamus, lymph node, testis, thyroid, and frontal cortex (Azevedo et al., 2014;Bhardwaj & Purohit, 2021).
The binding site of all the PDEs consists of a pocket and has a conserved GLN residue at position 726 (Manallack et al., 2005).The comparison and profiling of binding sites of known structures of PDEs led to determining the unique features of PDE10A.The PDE10A contains GLY at position 725, which is adjacent to conserved GLN, and the same is not the case among other PDE families (Verhoest et al., 2009).Additionally, the M-loop residues which form the back of this lipophilic pocket connects helix 14 and 15, and this loop is longer in PDE10A due to insertion of LYS/ARG-709 than other PDEs except for PDE6 and PDE5; as a result, the pocket in PDE10A is more profound than other PDEs (Chappie et al., 2012).Finally, to anchor an inhibitor inside the pocket, PDE10A has TYR-693, making a hydrogen bond.Apart from PDE10A, only PDE2 has this potential amino acid residue making hydrogen bond, but a leucine residue of PDE2 restricts access to this pocket (Verhoest et al., 2009).All these salient features of the substrate recognition site of PDE10A make it a potential drug target candidate for structure-based drug designing and discovery for countering related disorders.It must be taken into consideration that the reported structures of PDE10A in the protein data bank did not follow the consistent pattern for numbering the residues of the active site.Wyeth-Elbion-group referred to the conserved GLN as GLN-716, while in some other studies it was also referred to as GLN-726 (Kehler & Nielsen, 2011).
PDE10A dysregulation is thought to be associated with Huntington's disease as the inhibition of PDE10A (using TP-10, PDE10A-selective inhibitor, or by PDE10 knockout mice) altered the signaling pathway associated with Huntington's disease and provided beneficial effects in Huntington's disease models (Erro et al., 2021;Kleiman et al., 2011).Also, the administration of PDE10A inhibitors MP-10 and papaverine produces beneficial effects by improving glutamatergic and dopaminergic dysfunctions thought to be associated with schizophrenia (Grauer et al., 2009;Siuciak et al., 2006).But in the case of papaverine, a cross reactivity was observed with PDE3 family members which hinder its use to be a selective inhibitor for PDE10A (Helal et al., 2011;H€ ofgen et al., 2010).Other inhibitors from Lundbeck (AF11167), Pfizer (PF-02545920), and Takeda (TAK-063) were developed to act as potential antipsychotic agents, but none of them during Phase-II clinical trials is able to possess a robust antipsychotic efficacy (Menniti et al., 2020).Despite its association with several neurological disorders, it was also found to be a suitable drug target in several cancers like colon, breast, and lung cancers, but inhibitor designing was focused on neuronal disorders, and the functioning of PDE10A in tumorigenesis is yet to be explored (Russo et al., 2014;Zhu et al., 2015Zhu et al., , 2017)).The research area in still wide open; here, we screen a set of inhouse synthesized 1,2,3 triazoles using classical and enhanced computational methods to identify novel, potent and selective PDE10A inhibitors.In this study, we identified that 3h and 3j bind strongly with the active site of PDE10A without any significant change in the conformation of the protein.Hence, these molecules could be exploited as the novel, potent and selective inhibitors for PDE10A and further subjected to in-vitro and invivo studies.

Molecular docking
Molecular docking is considered one of the highly reliable computational tools to find the orientation of ligand molecules into the target protein's active site and understand the precise mechanism of selectivity/binding of substrate or inhibitor (Sack et al., 2016).The CDOCKER module of Discovery studio was utilized to perform molecular docking of PDE10A with standard and in-house developed molecules.CDOCKER is a grid-based docking tool that employs the CHARMm force field to dock molecules into the target protein's active site.The results were expressed as negative CDOCKER energy or score.The CDOCKER energy combines van der Waals, electrostatic, and hydrogen bond interactions between ligand molecules and target protein (Wu et al., 2003).For molecular docking, the already present inhibitor PFJ was removed from the active site of PDE10A, and the site of binding of PFJ was considered as a ligand-binding site for docking.A sphere of 10 Å radius centered at X ¼ 12.358, Y ¼ 19.950 and Z ¼ 43.041 was set to facilitate the interaction of ligands with the PDE10A.During docking, each ligand underwent 1000 dynamics steps, with ten rotating ligand orientations with respect to initial random conformations.Therefore, the in-house synthesized 1,2,3 triazoles and standard molecules, papaverine and PFJ, were docked into the active site of PDE10A.

Molecular dynamics (MD) simulations
Protein dynamics at the native state and interactions at atomic resolution are essential and can be monitored using a powerful computational method named Molecular Dynamics (MD) Simulations.Here, MD simulations were conducted on PDE10A complexed with papaverine, PFJ and selected complexes from molecular docking with the GROMACS 5.0.6 (Van Der Spoel et al., 2005).The force field utilized to simulate the PDE10A complexes was GROMOS96 43a1.All the simulations were performed on a high-end, in-house cluster machine.The gmxgrep command was employed to extract the respected ligands from the protein-ligand complexes of PDE10A from molecular docking and further subjected to the generation of topologies for both target protein and ligand molecules.The topology and other parameters for ligands including standard molecules were generated using the GlycoBioChem PRODRG server (Sch€ uttelkopf & Van Aalten, 2004).Further, to prepare the complexes for MD simulations, the generated topologies of the ligands from PRODRG and topologies of target protein prepared using pdb2gmx module of gromacs were rejoined.The selected complexes were placed in a cubic box having dimensions of 10 Å filled with water molecules.The gmxeditconf module was utilized to set box parameters, and subsequently, solvation was done with Simple Point Charge (spc216) water model using gmxsolvate module.The gmx genion script was then employed to maintain electroneutrality in the simulation box by adding appropriate sodium or chloride ions.
Now the system needs to attain minimum energy, and it was done using the steepest descent energy minimization.Further, NVT and NPT ensembles were employed in subsequent steps for 1000ps each at 310 K to equilibrate the system.Soon after the equilibration was done, the wellprepared system without any restrains was subjected to an MD production run of 100 ns with an integration step of 2fs using the leap frog integration algorithm.Throughout the simulation, to maintain a temperature of 310k particle-mesh Ewald method was employed (de Souza & Ornstein, 1999).A constant pressure of 1 bar is also maintained with a coupling constant of 1.0 ps by using the Parrinello Rahman pressure coupling method (Parrinello & Rahman, 1981).
For calculating long-range columbic and electrostatic interactions, particle mesh Ewald method was implemented, while for constraining the hydrogen bond length and other bonded interactions, SHAKE algorithm was implemented (Bhardwaj & Purohit, 2021).Finally, the RMSD, RMSF, and hydrogen bonds were analyzed using GROMACS utilities gmx rms, gmx rmsf, and gmx hbond, respectively, and further plot generation was done using Grace toolkit.

MM-PBSA (molecular mechanic/Poisson-Boltzmann surface area) calculation
Binding free energy for a protein and ligand complex was calculated from resultant MD trajectories using the MM-PBSA method.The equation which applies during MM-PBSA is as follow: Here, G (Complex) represents the total free energy of the protein ligand complex and G (Ligand) and G (Receptor) signify the free energies of ligand and protein, respectively.The solventaccessible surface area (SASA), polar solvation, Van der Waals and, electrostatic energies in combination represent the end state binding free energy.The selected protein ligand complexes were evaluated using the g_mmpbsa (Kumari et al., 2014) script to calculate binding free energies.

Entropy calculation by quasiharmonic method
In order to calculate the relative and absolute entropies, Schlitter suggested a method based on the Quasiharmonic (QH) approach.These entropies were obtained using MD simulation trajectories to evaluate the covariance matrix of Cartesian coordinates.According to Schlitter's paradigm, the computation of absolute entropy identifies the maximum possible entropy in quantum mechanics (Schlitter, 1993).

Steered molecular dynamics simulation
A slightly different method from conventional MD simulations is Steered molecular dynamics (SMD) simulations.A constant time-dependent pulling force is applied to unbind ligand/peptide from the target protein (Lesitha Jeeva Kumari et al., 2017).GROMACS package was utilized to perform SMD simulations.Solvation of the system is followed by adding ions to neutralize the system.Now the system was subjected to energy minimization to remove wrong contacts.Simultaneously a 500ps NPT equilibration was carried out.The target protein remained immobile during SMD simulation, whereas all the position restraints were removed from the ligand molecule.Finally, SMD simulation was carried out by pulling the ligand away from the target protein along the Z-axis over 500ps with a spring constant of 250 kJ/mol/nm 2 .Throughout the simulation, 0.01 Å/ps pull rate was maintained.
Umbrella sampling was performed on the chosen trajectories obtained from the pulling simulations.The center of mass (COM) distance obtained from the pulling simulation was used as the reaction co-ordinates for umbrella sampling simulations.An asymmetrical distribution pattern was followed for sampling windows with 0.10 nm window spacing for up to 3 nm COM separation and 0.2 nm spacing for the rest of the COM separation.As a result, total of 31 windows were collected, and further, each sampling window was subjected to 10 ns of MD simulation.Finally, the weighted histogram method was employed to analyze umbrella sampling outcomes (Li et al., 2016).

Results and discussion
PDE10A, a dual substrate phosphodiesterase with a preference for cAMP degradation, is an emerging drug target for neurological diseases due to its high expression level in the mammalian brain (Abdel-Magid, 2018;J€ ager et al., 2012).Hence, the hunt for identifying more selective and potential PDE10A inhibitors is still on.Here, we screened in-house synthesized molecules as potential candidates to inhibit PDE10A through molecular docking, followed by the classical MD, SMD, and umbrella sampling simulations.

Molecular docking
Molecular docking, a computational approach, provides significant information about the active site residues of target protein and interactions of molecules with the protein (Yadav et al., 2017).In this study, seven in-house synthesized molecules, 3f, 3g, 3h, 3j, 3i, 5d, and 5e were docked with the active site of PDE10A to identify potential inhibitors.Also, two cocrystallized ligands, papaverine and PFJ were taken as reference compounds to compare the potential of selected molecules against PDE10A.
The cocrystallized ligands papaverine and PFJ had a binding affinity of À 153.12 kJ/mol and À 179.17 kJ/mol, respectively.All the selected molecules showed better binding affinity than the standard molecules, ranging from À 153.32 kJ/mol to À 205.69 kJ/mol (Table 1).Papaverine formed hydrogen bond with amino acid residue GLN-716 and van der Waals interactions with the GLY-715, TYR-683, ALA-722, TYR-720, SER-667, LEU-665 and PHE-686 (Figure 1).The other interactions like Pi-Pi stacking, Pi-alkyl, Pi-sigma, Pi-sulphur and alkyl were also observed for papaverine complexed with PDE10A.The Pi-alkyl, Pi-Pi T shaped, Pi-sigma and Pi-Sulphur interactions comes in the broad category of non-covalent interactions and mainly involves charge transfer between the ligand molecule and amino acid residues of protein, thus responsible for intercalation of ligand molecules inside the active site of protein (Arthur & Uzairu, 2019).
PFJ, which is a cocrystallized ligand, also formed hydrogen bond with the GLN-716 and all other interactions depicted in Figure 1.Out of the seven selected molecules, 3i exhibited a better binding affinity of À 205.69 kJ/mol, which is the best among the selected molecules and better than papaverine and PFJ.The molecule 3i formed hydrogen bonds with two amino acid residues (SER-563 and ASN-562), while several residues formed van der Waals interactions (ALA-626, ASP-626, HIS-515, HIS-519, TYR-514, ASP-664, SER-667, VAL-668, GLN-716, TYR-683, MET-704, and ILE-701), PHE-686 and PHE-719 formed Pi-Pi Tshaped and Pi-stacking interactions respectively (Figure 1).The MET-703, which is an active site residue, formed Pi-sulphur interactions.The two other molecules, 3h, and 3j showed the binding affinity of À 189.98 kJ/mol and À 194.74 kJ/mol, respectively, with PDE10A.The molecule 3h formed hydrogen bonds with SER-563 and ASN-562, while 3j formed hydrogen bond with TYR-683.Almost similar kind of other interactions were observed for 3j with active site residues as that of 3h, except an additional Halogen interaction was made by GLY-718 in the case of 3j (Figure 1).
The molecules 3f, 3g, and 5d also showed better binding affinity than the cocrystallized molecules.The binding affinity of À 176.86 kJ/mol, À 177.11 kJ/mol, and À 177.73 kJ/mol was observed in the case of 3f, 3g, and 5d, respectively.The molecule 3f formed the hydrogen bond with the TYR-683, while no conventional hydrogen bond was observed for the other two molecules (3g, and 5d).The molecule 5e exhibited a binding affinity of À 153.32 kJ/mol, which is almost similar to papaverine and less than the other cocrystallized molecule PFJ.An invariant substrate-recognizing residue that is conserved in all PDEs is GLN-716 (Manallack et al., 2005), and it was evident from interaction analysis of triazole molecules that they are making interactions with this residue.Apart from it to confer selectivity to the inhibitor molecules, TYR-683 plays important role as it resides inside the selectivity pocket of PDE10A (Bhardwaj & Purohit, 2021).The 3j formed hydrogen bond with this residue and it also interacted with other triazole molecules with some other kind of interactions which is important to achieve selectivity.The interaction analysis also suggested that the papaverine and PFJ formed only one or two hydrogen bonds and mainly, formed van der Waals and other interactions.These other interactions including hydrogen bonds mainly account for the binding affinity of these molecular complexes.
As all the in-house synthesized molecules exhibited better binding scores than the standard molecules except 5e, and even good interaction pattern were observed for all the molecules with PDE10A; hence all the molecular complexes, including cocrystallized ligands, were subjected to MD simulations to access the stability of these complexes.

MD simulation analysis
The one of the significant limitations of molecular docking is its inability to distinguish between the conformational changes of the receptor or ligand; hence MD simulations has a part to play as it provides in-depth structural information about the protein-ligand complexes.It offers a complete information about the strength, overall stability, and flexibility of the protein-ligand/drug-receptor complexes (Modestia et al., 2019;Sharma et al., 2021).
Here, all the docked molecules (3f, 3g, 3h, 3i, 3j, 5d, 5e) and two standard molecules (PFJ and Papaverine) were subjected to MD simulations for a time period of 100 ns.After that, various time-dependent MD simulation analyses were carried out.

Root mean square deviation (RMSD)
Root mean square deviation (RMSD) is an important parameter that determines the structural stability of protein-ligand complexes (Kumar Bhardwaj et al., 2022).Resultant trajectories after MD simulations were analyzed for molecules complexed with PDE10A to compute the RMSD values.It was observed that the RMSD for 3h, 3j, PFJ, and papaverine were stabilized and equilibrated at 30ns for a fraction of time.Afterwards, the RMSD for papaverine was suddenly increased and showed compratively higher deviations throughout the simulation time period.Also, the other standard molecules PFJ showed the deviation in RMSD at around 55-60ns and at 75ns of the time period.The RMSD for two triazole derivatives, (3h and 3j) remained consistent for almost the whole 100ns of time period once they stabilized at 30ns, while a minute fluctuation was observed for 3j at around 80ns (Figure S2).The average RMSD values for PFJ, papaverine, 3h and 3j were 0.366 nm, 0.425 nm, 0.376 nm, and 0.385 nm, respectively.The RMSD values for other molecules, i.e. 3f, 3g, 5d, 5e, and 3i, converge at around 40ns and remain consistent up to 60ns.After that, slight fluctuations were observed for each molecule (Figure S3).The overall RMSD analysis suggested that the 3h complexed with PDE10A was more stable than the other molecules, but no drastic changes were observed even for other molecules.The RMSD obtained from the MD trajectories was comparable with the previous studies (Al-Nema et al., 2022;Bhardwaj & Purohit, 2021), which also suggested the robustness of our methodology in finding suitable inhibitors for PDE10A.

Root mean square fluctuations (RMSF)
As the RMSD provides information about the C-a atom, further, root mean square fluctuations (RMSF) were carried out to analyse the effect of amino acid residues on the functionality and stability of the protein complexes (Kumar et al., 2022).The movement of amino acids for all the complexes, including standard molecules were determined for a period of 100ns.The critical residues of the active site (highlighted using black dotted lines in Figure 2) showed slight fluctuations upon binding of 3h and 3j while starting residues at around 701 amino acid exhibited more considerable changes for PFJ and papaverine.The binding of 3h and 3j to the active site of PDE10A enhanced the rigidity of amino acid residues as compared to PFJ and papaverine which results in more stable complexes.As the molecular complexes of3h and 3j with PDE10A, remained stable and slightly fluctuated throughout the simulation period, this stability might be the possible explanation for the enhanced estimated binding affinity (Table 1) of these molecules as that of the standard papaverine and PFJ.The RMSF for other hit molecules was depicted in Figure S4.As seen from the figure, all the molecular complexes showed higher amino acid fluctuations, and it was also supported by RMSD analysis where significant changes in RMSD values were observed for these molecules.

Hydrogen bond analysis
For a protein-ligand complex to remain stable, it must form a hydrogen bond between the ligand molecules and the target protein.Hydrogen bonds are used in medication design to increase specificity because of their severe distance and geometric restrictions (D.Chen et al., 2016).Hydrogen bonds provide the foundation for many processes used to study protein-ligand complexes, from docking and scoring to in-depth analyses of binding activities and catalytic activity (Nittinger et al., 2017).In order to obtain pertinent results for drug discovery research, computational techniques were utilized in several studies to access the number of hydrogen bonds.
To further validate our results, we also performed the hydrogen bond analysis of all the selected complexes (3f, 3g, 3h, 3j, 5d, 5e and 3i), including standard molecular complexes (papaverine and PFJ).The papaverine formed only one hydrogen bond up to 40ns of the time period, and after that, none of the hydrogen bond formed by papaverine with PDE10A throughout the entire simulation of 100ns.In contrast, the other standard molecule, PFJ, comprised a maximum of three hydrogen bonds, and for starting 10ns of simulations, it showed two hydrogen bonds.After that, up to 75ns, very few hydrogen bonds were observed (Figure S5).PFJ formed maximum hydrogen bonds from 80ns to 100ns of time.The two hit molecules, 3h and 3j, followed more consistent pattern in hydrogen bond formation.These two molecules formed a maximum of four hydrogen bonds, and in the case of 3h, hydrogen bonds were consistently formed throughout the simulation period, while 3j consistency was seen up to 65ns of the time period.
The hydrogen bond pattern for other hit molecules was depicted in supplementary Figure S6.As seen from the figure, the 3i formed maximum of eight hydrogen bonds, 3f, and 5d fromed maximum of six hydrogen bonds, while 5e formed a maximum of seven hydrogen bonds.Moreover, it is crystal clear from the figure that none of the molecules form consistent hydrogen bonds throughout the simulation time period.The fluctuation of RMSD and RMSF for these hit molecules might be the reason for inconsistent hydrogen bond pattern.

Secondary structure element analysis
The secondary structure after ligand binding to PDE10A was evaluated, the selected hit molecular complexes and standard molecules were subjected to Dictionary of Secondary Structure of Proteins (DSSP) analysis.The PDE10A complex with cocrystalized ligand at around 695 to 716 residue position, includes important residues of the active site of PDE10A attained a helix conformation at around 65ns and it remained as helix upto 80ns of simulation time period while, turn, bend and coil conformations were observed throughout the simulation (Figure S7a).The papaverine complex exhibit coils, bends, turns and few traces of 3-helix during simulations at same residue positions (695 to 716) (Figure S7b).For molecule 3h complexed with PDE10A 3-helix conformation remains dominated at around 695 to 716 residue position and other structural conformations like bend, coil and turn were also seen (Figure S7c) in contrast for 3j the 3-helix conformation was seen after 80ns and remained upto 100ns (Figure S7d).So it was observed that the binding of 3h with PDE10A allowed it to attain more helix conformations as compared to cocrystallized ligand and other molecules.As we know the helix structure provides stability to protein complexes hence the 3h formed much stable complex with PDE10A and due to augment in stability the molecule 3h remained firm inside the active site and had stronger interactions with the active site residues of target protein as compared to papaverine, PFJ and other molecules.
The secondary structure changes were also observed at other positions within the protein structure which are far away from the active site.The starting residues from 453 for PDE10A after binding of co-crystallized ligand PFJ attained a helix conformation up to 25ns.After that, helix completely lost and coil conformation dominated with few traces of bend conformations (Figure S7a).The same residues for other standard molecule papaverine showed helix configuration up to 55ns; after that, a few traces of bend confirmations followed by coils for the rest of the simulation time (Figure S7b).In the case of 3h complex, these starting residues attained helix and coil conformations upto 50ns and after that coil and turn conformations remained dominated, while for 3j, the helix conformations were predominantly found.For amino acid residues 595-605, beta-sheet conformation was observed at around 45-65ns and 75-100ns of the time for PFJ complex.In the case of papaverine at similar residue position bend conformations with few traces of coil and turn conformations was observed.
The molecular complex of 3j at about 595 to 605 amino acid residues attained turn and bend conformations while a minor fraction of beta-sheet was observed at around 15-20ns of the time.The 3h molecular complex at similar positions attained diverse conformations as coil, turn and bends were predominantly found.After 80ns, few traces of 5-helix and around 95ns to 100ns beta-sheet interspersed in between turns were observed.
DSSP profiles for other hit molecular complexes were illustrated in supplementary Figure S8.Collectively helix conformation was dominated throughout the simulation time period with little fluctuations and no significant conformational changes in PDE10A was observed due to binding of selected molecules as compared to the native ligand PFJ.As a result the PDE10A complexes with hit and standard molecules remain stable throughout the simulation period.Further to strengthen the evidence from above analysis, all the hit molecular complexes were subjected to MM-PBSA and free energy landscape analysis (FEL).

MM-PBSA analysis
MM-PBSA is one of the popular approaches for accessing the binding free energy of small ligand molecules and biological macromolecules in a dynamic environment (Rifai et al., 2019).The binding free energy of 3f, 3g, 3h, 3j, 5d, 5e, 3i, papaverine, and PFJ complexed with PDE10A was calculated for 100ns of MD simulation trajectories.The polar solvation energy, SASA, van der Waal energy, and electrostatic energy make up the end-stage binding free energy calculated in MM-PBSA.
Papaverine and PFJ, the standard molecular complexes used in the study, exhibited the binding free energy of À 202.030þ/À 16.732 kJ/mol and À 138.764 þ/À 20.939 kJ/mol, respectively (Figure 3) (Table 2).The Electrostatic energy, Van der Waal energy, and SASA energy predominantly contributed to the final, binding free energy of papaverine complex.In contrast, only Van der Waal energy and the SASA energy contribute to the end-stage binding free energy of the PFJ complex.Out of the selected hit molecules 3h complex exhibited binding free energy of À 240.642þ/À 17.051 kJ/mol while the molecular complex of3j exhibited À 201.406þ/À 35.658 kJ/mol binding free energy (Figure 3) (Table 2).Electrostatic energy, Van der Waal energy, and the SASA energy contributed to the final, binding free energy for both the molecular complexes.Both the molecules exhibited lower binding free energy (higher affinity) than the molecule #3 suggested as a potential candidate to inhibit PDE10A in a previous study (Bhardwaj & Purohit, 2021).Table 2 also depicts the binding free energy for the rest of the selected hit molecular complexes.Graph for the rest of hit molecules was depicted in Figure S9.
From selected hit molecules, the 3h complex showed the least binding free energy and is even better than both the standard molecules (papaverine and PFJ), followed by the 3j complex, which performs better than PFJ and have almost similar binding free energy as that of papaverine.The significant difference observed in the binding free energy of the 3h and 3j complex is due to the electrostatic energy (Table 2).
The selectivity pocket residue, TYR-683 (Verhoest et al., 2009), also showed the contribution energy of À 4.3095 kJ/ mol, À 1.7032 kJ/mol, À 0.7025 kJ/mol, and À 10.631 kJ/mol in the case of PFJ, papaverine, 3h, and 3j, respectively.If we talk about the selectivity, the PDE10A has GLY-715 adjacent to conserved GLN-716 and this GLY-715 had favorable contribution energy for molecular complexes of 3h and 3j as compared to standard molecules.The contribution energy for 3f, 3g, 3i, 5d, and 5e was depicted in Figure S10.
Above analyses clearly depicted that both the top hit molecules 3h and 3j accommodate well inside the active site and 3j also had significant interactions with the selectivity pocket residue while 3h interacted strongly with the key residues of the active site, (MET-703 and PHE-719).Supplementary Figure S10 depicts the graph for contribution energy of residues for other selected molecules.Additionally, using a QH technique, the configurational entropy for the selected protein-ligand complexes including standard molecules was estimated.Entropic effects are essential for protein folding, ligand binding, molecular interaction, and stability.All the molecular complexes of PDE10A with 3f, 3g, 3h, 3i, 3j, 5d and 5e had comparable entropic values as that of papaverine and PFJ (Table S1).

FEL analysis
The principal components named PC1 and PC2 were taken into account to generate the free energy landscape plots for PDE10A complexes with standard and hit molecules.GROMACS scripts called g_anaeig, g_sham, and g_covar were utilized to perform these tasks (Majumder & Mandal, 2022).
The FEL tells about the global minimum energy conformation of a receptor-ligand complex.A single cluster of conformation represent a stable and strong protein ligand interactions and if multiple clusters of energy minima was observed then it means that the interactions between a protein and ligand is weak or unstable.The contour 2D and 3D plots of papaverine, PFJ, 3h, and 3j complexes have been depicted in Figure 5.A distinct pattern of FEL was observed for every complex.From the Figure 5b, it was clear that the papaverine complex showed a single noticeable energy minima, and also three energy basins were observed.The other inhibitor PFJ had demonstrated at least two broad global energy minima (Figure 5a).The 3h complex exhibited single deep and wide global energy minima (Figure 5c), while the 3j complex had observed single narrow energy minima, and a few scattered energy minima were also observed (Figure 5d).It suggested that the hit molecules 3h and 3j form energetically favorable, stable and strong complexes with PDE10A.
All these analyses showed that the molecules 3h and 3j, among the other hit molecules, perform much better during the conventional MD simulations.Both the molecular complexes remain stable in the binding site of PDE10A as that of standard molecules papaverine and PFJ, and even among two, 3h had been the potential one as suggested through MM-PBSA analysis.Both the molecules can be exploited as potential PDE10A inhibitors.In order to provide further strength to these outcomes, the 3h and 3j with PDE10A, along with standard papaverine and PFJ, were subjected to SMD and umbrella sampling simulations.

Unbinding mechanism of proteinligand complexes
The binding and unbinding mechanisms of biological macromolecules with small molecules or peptides in structural biology are monitored through an advanced version of MD simulations, regarded as SMD simulations (Rodriguez et al., 2015).SMD is similar to atomic force microscopy, where an external pulling force is applied to pull the ligand away from the target protein's active site (Singam et al., 2014).It can be noticed that the binding affinity of ligand molecules shares a direct correlation with the pulling force.SMD offers more effective ranking procedures while reducing computational expense, and previously utilized in various computationalbased drug screening studies (Le et al., 2010;Mai & Li, 2011; V. K. Singh et al., 2017).Hence the selected ligands 3h, 3j, and two standard molecules, papaverine and PFJ, complexed with PDE10A subjected to SMD analysis.SMD simulation results are depicted in Figure S11.It was inferred from the results that the time required for 3h, 3j, papaverine, and PFJ to move out of the PDE10A ranges between �100 and 350ps and force ranges from �190 to 325kJ/mol/nm.The molecule 3h, which had the highest binding free energy, required a larger force to pull away from the complexed state, followed by PFJ, 3j, and papaverine.At the start of pulling simulations, a rise in pulling force was seen for all the selected complexes, and pull force reached the maximum of �300kJ/mol/nm in the case of 3h, and also it takes more time to pull it away from the binding site as clearly inferred from the graph (Figure S11).It strongly indicated that out of the four selected complexes, including standard molecules, the 3h is confined strongly within the active site of PDE10A, which emphasizes it as a potential PDE10A inhibitor.Further number of contacts formed in PDE10A molecular complexes during pulling simulation was analyzed to strengthen the candidature of 3h as a potential inhibitor for PDE10A (Figure S12).At the initial stages of simulations, the maximum number of contacts made by papaverine, and PFJ was �820 and �1010, respectively, while initially, 3h made �1310 and 3j showed �1160 contacts with PDE10A.As the pulling simulation progressed, the number of contacts increased up to �180ps, and a maximum of �1550 was observed in between for 3h, while the same is not the case with the 3j, PFJ, and papaverine.It further validate the SMD results as 3h look confined much better within the active site of PDE10A, and as a result larger force is required to pull it away from the active site.

Umbrella sampling
A technique named umbrella sampling, an extended version of SMD simulations was utilized to monitor the binding or unbinding of protein-ligand complexes through calculating the potential of mean force (PMF) (R. Singh et al., 2021;Tung et al., 2019).The enhanced sampling method had been utilized in several studies, providing an important insight into the free energies of inhibitors bound to respective targets (J.Chen et al., 2017aChen et al., , 2017b)).The PMF graph plotted for two best hits, and standard molecules were depicted in Figure 6, and furthermore, individual free energy values were calculated manually from the plotted graph (Table 3).
It was evident from the graph that the energy minima for papaverine and PFJ with PDE10A takes place at about center of mass (COM) distance of 1.11 nm and 1.08 nm, respectively, while for best hit molecules, 3h and 3j fall at about COM distance of 1.5 nm and 1.01 nm respectively.The energy minima at larger COM distance for 3h compared to standard molecules hints about its retention for a larger period of time within the active site of PDE10A and the possibility of having stronger interactions.Even the other hit molecule, 3j, had almost similar COM distance as that of the papaverine and PFJ.Therefore, the possibility of 3j being a potent inhibitor for PDE10A should not be excluded.
The other parameter from umbrella sampling simulations is the binding free energy (difference between highest and lowest PMFvalues).The binding free energy for papaverine and PFJ with PDE10A, was À 70.60 kJ/mol and À 55.35 kJ/mol.For 3h, the binding free energy of À 81.65 kJ/mol was observed, while for 3j, a much lower value of binding free energy,(À 101.69 kJ/mol) was observed.
The molecule 3h results agreed with previous MM-PBSA analysis as it performed much better than the standard molecules, while slight deviation was observed for 3j molecule in umbrella sampling simulations.Hence the possibility of 3j to inhibit the PDE10A cannot be excluded, and both the molecules (3h and 3j), could be assessed as potential lead molecules to inhibit PDE10A in related disorders.

Drug-likeliness properties
The ADMET properties plays a crucial role because a large percentage of compounds failed clinical trials due to subpar ADMET properties (Lagorce et al., 2017).All the triazole derivatives included in the study were subjected to ADMET analysis and compared with the standard papaverine and PFJ.The complete information about ADMET parameters is depicted in Table S2.The Blood Brain Barrier (BBB) is an essential aspect of molecules designed for neurological applications.The potential lead molecules 3h and 3j were predicted to cross BBB and be non-toxic, as indicated by the AMES toxicity test.Compounds with less than 30% value for human intestinal absorption are poorly absorbed (Guleria et al., 2021).All the triazole molecules had good potential to be absorbed in the intestine, similar to standard molecules.Another critical aspect of the drug-like molecule is its metabolism when administered within the body.Cytochrome P450 is responsible for the metabolism of drugs, and any drug inhibiting cytochrome P450, cannot be metabolized in the body.The predicted values for selected molecules in the study suggested that none of the molecules had the potential to inhibit CYP2D6/CYP3A4, which are important isoforms of Cytochrome P450.

Conclusion
In the current investigation, extensive computational methods including molecular docking, conventional MD simulations, MM-PBSA, steered molecular dynamics, and Umbrella sampling simulations were applied to identify novel PDE10A inhibitors.In-house synthesized molecules were docked with PDE10A using the CDOCKER module of the discovery studio.The selected molecules exhibited good binding affinity during molecular docking as that of standard/cocrystallized ligand molecules, so for optimization of the inhibitory potential of the molecules, they were subjected to MD simulations and MM-PBSA analysis.The comparative analysis with standard molecules suggested that the molecular complexes of 3h and 3j with PDE10A had stable trajectories during RMSD analysis.The MM-PBSA analysis revealed that the 3h had the least binding free energy followed by 3j.Furthermore, the pulling simulations suggested that the unbinding of 3h required to overcome a higher energy barrier than the rest of the molecules.In contrast, the umbrella sampling simulations revealed better binding free energy for 3j followed by 3h.Overall the results obtained from the study revealed that the molecules 3h and 3j could be exploited as potential lead molecules in in-vitro and in-vivo experiments against PDE10A.

Figure 2
represents the RMSF values of 3h, 3j, PFJ, and papaverine complexes.It is clear from the figure that the amino acid residue fluctuations were more in the case of standard molecules PFJ and papaverine than the 3h and 3j.The RMSD fluctuations for papaverine and PFJ might be the reason for their higher RMSF values throughout the simulation period.The amino acid residues at position 603, 604 and 640-643 showed higher fluctuations in the case of 3h, but none of the residues constitute the active site of PDE10A.The molecule 3j, showed few amino acid fluctuations throughout the simulation as shown by other molecules.

Figure 3 .
Figure 3. Graphical plot representation of MM-PBSA calculations in the form of average total energy for the binding of papaverine, PFJ, 3h and 3j.

Figure 4 .
Figure 4. Contribution energy of different amino acid residues towards PDE10A complexes.

Figure 6 .
Figure 6.Potential mean force (PMF) plot as a function of reaction co-ordinates for Papaverine, PFJ, 3h and 3j complexed with PDE10.

Table 1 .
CDOCKER interaction energy of standard inhibitor (papaverine and PFJ) and selected molecules.

Table 2 .
Binding free energy of complexes with standard ligands and selected molecules.

Table 3 .
Binding free energy from umbrella sampling.