Insights into the potential of withanolides as Phosphodiesterase-4 (PDE4D) inhibitors

Abstract Medicinal herbs have been used as traditional medicines for centuries. The molecular mechanism of action of their bioactive molecules against various diseases or therapeutic targets is still being explored. Here, the active compounds (withanolides) of a well-known Indian medicinal herb, Ashwagandha (Withania somnifera), have been studied for their most potential therapeutic targets and their mechanism of action using ligand-based screening and receptor-based approaches. Ligand-based screening predicted the six top therapeutic targets, namely, Protein kinase C alpha (PRKCA), Protein kinase C delta (PRKCD), Protein kinase C epsilon (PRKCE), Androgenic Receptor (AR), Cycloxygenase-2 (PTGS-2) and Phosphodiesterase-4D (PDE4D). Further, when these predictions were validated using receptor-based studies, i.e. molecular docking, molecular dynamics simulation and free energy calculations, it was found that PDE4D was the most potent target for four withanolides, namely, Withaferin-A, 17-Hydroxywithaferin-A, 27-Hydroxywithanone and Withanolide-R. These compounds had a better binding affinity and similar interactions as that of an already known inhibitor (Zardaverine) of PDE4D. These results warrant further in-vitro and in-vivo investigations to examine their therapeutic potential as an inhibitor of PDE4D. Communicated by Ramaswamy H. Sarma


Introduction
Adenosine-3 0 ,5 0 -cyclic phosphate (cAMP) and guanosine-3 0 , 5 0 -cyclic phosphate (cGMP) function as an intracellular secondary messenger for cellular signaling. Their role in the inflammatory process, cell cycle progression and various cancers are well reported (Levy et al., 2011;Raker et al., 2016). Phosphodiesterases (PDEs) are enzymes that hydrolyse the cyclic AMP (cAMP). Based on the amino-terminal domains and their functions, the PDEs are divided into 11 groups (PDE1-11). Out of them, PDE4, 7, and 8 are specific to cAMP (Azevedo et al., 2014). Further, PDE4 are subdivided into various isoforms, namely, PDE4A, PDE4B, PDE4C and PDE4D (Azevedo et al., 2014;Li et al., 2018). PDE4D is one of the crucial isoforms responsible for the degradation of cAMP without changing the cGMP concentration (Massimi et al., 2019;Ricciarelli & Fedele, 2015). In recent studies, the role of PDE4D in various cancers has aslo been reported (Hsien Lai et al., 2020;Powers et al., 2015). A lower level of cAMP helps in the survival and growth of cancer cells (Fajardo et al., 2014). These studies have also revealed that an increase in the expression of the PDE4D level was found in lung, gastric, breast and ovarian cancer cell lines, and targeting PDE4D could be the potential therapeutic approach for multiple cancers (Lin et al., 2013). It has been shown that in A549 lung cancer cells, transforming growth factor-b (TGF-b) stimulation increased the PDE4D expression and activity, which spikes the epithelial to mesenchymal transition. It was concluded that PDE4D functions as a tumor-promoting factor and a unique targetable enzyme of cancer cells (Kolosionek et al., 2009). Most of the inhibitors for the PDE4D, such as rolipram, cilomilast, Zardaverine and NVP-ABE171, are used to treat respiratory and inflammatory disorders such as asthma (Dyke & Montana, 2002;Lipworth, 2005). However, in recent in-vivo and in-vitro studies, it has been reported that PDE4D inhibition by Zardaverine and Rolipram could be helpful in various cancer inhibition as well (Goldhoff et al., 2008;Sun et al., 2014;Yeo et al., 2014). Further, the role of cAMP and cGMP in brain development has also been reported in the literature (Delhaye & Bardoni, 2021;Son et al., 1998). Some PDEs like PDE2A, PDE3, PDE4D and PDE10A inhibitors have been reported to treat neurodevelopmental diseases such as autism and intellectual disability (Delhaye & Bardoni, 2021;Richter et al., 2013;Zamarbide et al., 2019). Apart from these, the role of PDE enzymes in regulating cAMP in spermatozoa is also well known (Balbach et al., 2018). The level of the cAMP plays a direct role in the mortality of the sperms (Freitas et al., 2017). In some previous studies, pentoxifylline (PTX) has been shown to inhibit PDE4A, PDE4D, and PDE10A, to enhance sperm function (Shen et al., 1991). However, the weaker binding affinity of the PTX has been a major limitation; hence, more potent drugs and therapeutics need to be developed (Balbach et al., 2018).
Various bioactive compounds from medicinal herbs have been studied in the past as PDE inhibitors. In one of the studies, Moracin M, a phenolic compound from Morus alba L, was found to be a potent inhibitor of PDE4D2 and PDE4B2 through molecular docking and MD simulation (Chen et al., 2012). In another experimental study, where phenolic compounds extracted from the herb Selaginella pulvinate were found to be potent against PDE4 having IC50 in the range of 0.11-5.13 mM (Liu et al., 2014). Similarly, three natural resveratrol analogs were found to have inhibitory action against PDE4D with IC50 < 100 mM (Zhao et al., 2013). Likewise, the compounds from leaf skin of Aloe barbadensis Miller and Gaultheria yunnanensis are also reported to be effective against PDE4D (Cai et al., 2016;Zhong et al., 2015). All these studies have been performed using the concept of the one protein target, where only one target protein has been focused and investigated. However, the major limitations of these studies were the off-target effects of the chosen compounds. In this study, we have chosen the Indian Ayurvedic herb Ashwagandha (Withania somnifera) that has been extensively used in traditional home medicine systems. It has been reported to have prophylactic and therapeutic activities Kalra et al., 2021;Munagala et al., 2011). We attempted to screen the top protein targets of withanolides using the inverse screening method and then using the structure-based molecular dynamics study; we tried to find the best target and its crucial interactions with withanolides to minimize the off-target effects. We have screened the well characterized withanolides against all the possible targets through ligand shape-based inverse screening using the Swiss Target Prediction server. The top common protein targets of withanolides were then selected through network analysis. Further, all the top targets were investigated for binding interactions, energy and dynamics with selected withanolides at the structural level. It was found that out of all the top selected protein targets, only PDE4D was best interacting with withanolides. Therefore, we examined the potential of natural compounds (withanolides) found from Ashwagandha against PDE4D. We took Zardaverine as the positive control for this study, whose IC50 value (0.8 mM) for PDE4D has already been reported (Schudt et al., 1991). The molecular docking study showed that among selected withanolides, Withaferin-A, 17-Hydroxywithaferin-A, 27-Hydroxywithanone, and Withanolide-R had higher docking energy than the Zardaverine. Further, the Molecular Dynamics (MD) Simulation results, along with Molecular Mechanics Generalized Born Surface Area (MM/GBSA) free binding energy calculations, indicated that withanolides had an approximately equal binding affinity towards PDE4D in comparison to Zardaverine.

Ligands retrieval, filtering, and inverse virtual screening
A total of 26 active compounds of Ashwagandha were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/). Further, the SMILES notation was used to calculate the ADMET properties of the compounds using the Swiss ADME web server (http://www.swissadme.ch/). The compounds which were violating any of the Lipinski rules were filtered out. Finally, 20 compounds were selected for the descriptorbased inverse virtual screening of therapeutic targets using the Swiss Target Prediction server (http://www.swisstargetprediction. ch/) (accessed in November 2020). The list of the compounds is provided in Table S1. The Swiss Target predictions are based on the similarity principle of the known ligands. The server was updated in 2019 by training the machine learning model with additional datasets. The updated tool has been reported to correctly identify at least one correct human target in the top 15 predictions for >70% of external compounds (Daina et al., 2019). Therefore, the top 15 target information was collected for each of the 20 Ashwagandha compounds that have been listed in Table S2.

Network analysis of the selected compounds-targets
The top 15 selected targets of each of the 20 compounds were used for network analysis to evaluate the best common targets and their interacting compounds using the Cytoscape software (version 3.7.2) (Shannon et al., 2003). The whole network was treated as directed, and Cystoscope's network analyzer plugin was used to analyze the network. Node, edge size and color were set as 'low values to small sizes' and 'low values to bright colors'.

Preparation of the protein-ligand complex
The 3D structure of Protein kinase C alpha (PRKCA), Cyclooxygenase-2 (PTGS2), PDE4D, and Androgenic receptor (AR), resolved using X-ray diffraction, were retrieved from Protein Data Bank (PDB) having IDs-3IW4, 5IKR, 1MKD and 1E3G, respectively (Lee et al., 2002;Matias et al., 2000;Orlando & Malkowski, 2016;Wagner et al., 2009). The crystal structures for Protein Kinase C Delta (PRKCD) and Protein Kinase C Epsilon (PRKCE) were not available in the PDB, and hence they were modelled using Swiss Homology web server using PRKCh as a template. The GMean score of the modeled structures was 0.35 that depicted a well-modeled structure. Further, the modelled structures were prepared using a protein preparation wizard and simulated for 50 ns. The structure from the last frame of the converged simulated trajectory was used for further studies. The minimized and simulated modeled structures were then validated using the Ramachandran plot using the PROCHECK server (Laskowski et al., 1996). The structures were then prepared using the protein preparation wizard of the Schrodinger suite (Madhavi Sastry et al., 2013;Schr€ odinger, 2020). The preparation steps involved the addition of the missing disulfide bonds, removal of water molecules, the addition of missing hydrogen atoms, filling of missing amino acids side chains and optimization of hydrogen bonds. After preparing the proteins, the 20 selected Ashwagandha compounds were prepared using the Ligprep module of the Schrodinger suite (Madhavi Sastry et al., 2013; Schr€ odinger, 2020).

Grid generation and molecular docking
The prepared protein structures were used for generating the grid for molecular docking. For PRKC isoforms (PRKCA, PRKCD and PRKCE), the grid (10 Å 3 ) was generated by selecting identical amino acid residues in all three proteins, as they had identical catalytic residues interacting with the positive control Sotrastaurin. The grid was generated taking the centroid of the interacting residues Phe358, Lys368, Met417, Glu418, Tyr419, Val420, Asp424, Asp467, Asn468 and Asp481 [PDB ID 3IW4] ). In the case of PDE4D, the PDB structure (1MKD) had a bound inhibitor (Zardaverine) and hence the coordinate of the inhibitor was used to generate the grid. Similarly, the bound inhibitor was chosen for grid generation in PTGS-2 (PDB ID: 5IKR) and AR (PDB ID: 1E3G). The prepared ligands were docked at the binding site by choosing the generated grid using the glide Extra precision (XP) docking algorithm with a flexible ligand sampling option. All other options were kept as default (Schr€ odinger, 2020).

MD simulations and analysis of the trajectories
The docked complexes were further taken for MD simulations to investigate the stability and crucial interactions between the protein-ligand complexes using the Desmond MD tool integrated with Maestro Schrodinger software (Schr€ odinger, 2020). Each system was solvated with the TIP3P water model in an orthorhombic periodic boundary box (Jorgensen et al., 1983). To prevent the interaction of the protein complex with its own periodic image, the distance between the complex and the box wall was kept at 10 Å. The system was then neutralized by adding the appropriate number of Na þ /Cl À ions depending on the protein-ligand complex using OPLS3e forcefield (Roos et al., 2019). The energy of the prepared systems was then minimized by running 100 ps low-temperature (10 K) Brownian motion MD simulation (NVT ensemble) to remove steric clashes and move the system away from an unfavourable high-energy conformation. Further, the minimized systems were equilibrated in five steps in NVT and NPT ensembles using the 'relax model system before simulation' option in the Desmond Schrodinger suite. The equilibrated systems were then subjected to 100 ns unrestrained MD simulations in NPT ensemble with 300 K temperature maintained by Nose-Hoover chain thermostat, constant pressure of 1 atm maintained by Martyna-Tobias-Kelinbarostat and an integration time step of 2 fs with a recording interval of 100 ps. The simulated trajectories were then analyzed for the Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF) and the number of hydrogen bonds between the complexes throughout the simulations using the Simulation event analysis tool of the Desmond suite (Schr€ odinger, 2020). Using Desmond's Simulation interaction diagram tool, the fraction of the interactions between the protein and ligand was calculated. Further, to compare the essential dynamics of the complexes, the Principal Component Analysis (PCA) of the simulated trajectories was computed using the Schrodinger script ($SCHRODINGER/run trj_essential_dynamics.py). Finally, to assess the binding affinity of the ligands toward the receptors, MM/GBSA binding free energy was calculated using the Prime module of the Schrodinger, the details of which are described elsewhere (Genheden & Ryde, 2015;Malik et al., 2021;Zhang et al., 2017).

Results and discussion
3.1. Ligand-based inverse screening predicts protein kinase C isoforms, AR, PTGS-2 and PDE4D as therapeutic targets of withanolides In classical structure-based drug discovery, a library of compounds is screened against a single therapeutic target of interest, and in this process, the off-target effects of the screened ligands are ignored.
Here, the natural compounds withanolides were used for ligand-based inverse screening of the best potential therapeutic targets. For this study, the structures of 26 withanolides were retrieved from the PubChem database. These 26 withanolides were chosen based on their structural availability and previously known bioactivity. Out of the 26 withanolides, 17 withanolides were chosen based on the recent study that reported the activity of withanolides against COVID-19 (Coronavirus Disease-19) (Khanal et al., 2021) and the rest of them were mined through individual literature searches. The selected withanolides were screened for their ADMET properties based on the Lipinski rule of five (Lipinski, 2000). Any compound violating any of the rules was screened out from the study. After the screening, 20 compounds out of 26 were found to follow all the Lipinski rules. These 20 compounds were then used for the ligand-based inverse screening of the best potential therapeutic targets using the SwissTarget server. For all the 20 compounds, their top 15 targets predicted by the server were retrieved. To find the best probable targets for the withanolides, the network for the Plant-Compounds-Targets was built using Cytoscape (version 3.7.2). The total number of nodes were 79, with a network density of 0.1, network diameter 4 and characteristic path length 2.8. Further, the interaction of withanolide Q with CDC7 (Cell division cycle 7) had the highest edge betweenness of 156. PRKCE, AR, PRKCA had the average shortest path length of 1.75. Then the top six protein targets of withanolides were identified based on the degree. The degree showed the number of connected nodes with the individual nodes. The top target proteins had more than 15 interactions with the withanolides in the network. It was found that AR, PRKCA, and PRKCE had the degree 19, followed by PTGS2, PDE4D and PRKCD, which had 17. Hence, these top 6 proteins are suggested to be hubs with which most withanolides interact (Figure 1). These predictions by SwissTarget Server were based on a trained machine learning algorithm that predicts therapeutic targets for any new ligand based on the previously known targets of similar ligands (Daina et al., 2019). The predicted targets were taken for further validation using molecular docking and classical MD simulations. The major protein group that was predicted to be a target for the withanolides was Protein kinase C (PRKC) that belonged to the family of serine/threonine kinases. PRKC has 10 isoforms that are categorized into 3 groups based on the requirement of cofactors for their activation (Steinberg, 2008). PRKCA belonged to the conventional group requiring Ca þ2 and DAG for its activation, while PRKCD and PRKCE were the novel isoforms requiring only DAG (G eczy et al., 2012;Steinberg, 2008). These PRKC isoforms play a crucial role in early T-cell activation and inhibition of their activity . Sotrastaurin (AEB071) has been widely reported to inhibit the phosphorylation of PRKC isoforms , and withanolides have also been shown to inhibit the autophosphorylation in kinases (Malik et al., 2021;Santhi & Aishwarya, 2011). Hence, in this study, a pan-kinase inhibitor, Sotrastaurin, has been used as a positive control against PRKC isoforms to compare the withanolides. For PRKCA, it was found that none of the top docked withanolides; namely, Withanolide-Q (À6.97 kcal/mol), 27-Hydroxwithanone (À6.95 kcal/mol), and 27-Hydroxywithanolide-B (À6.35 kcal/mol) with PRKCA was not even close in the docking energy in comparison to the already reported pan-PRKC inhibitor, Sotrastaurin (À13.84 kcal/mol). Hence, in the structure-based analysis, it was found that withanolides were not potent against PRKCA as compared to the known inhibitor Sotrastaurin. Then, the withanolides were investigated against modelled PRKCD and PRKCE. When these simulated modelled structures were validated through the Ramachandran plot, both PRKCD and PRKCE showed that more than 98% of the residues were in the most favourable and allowed region and that their simulated RMSD plot further showed a converged trajectory within 20 ns of simulation ( Figure S1). In the case of PRKCD, it was found that top docked withanolides; namely, Withanolide-A (À5.44 kcal/ mol), 17-Hydroxywithaferin-A (À4.99 kcal/mol) and 27-Hydroxywithanone (À5.09 kcal/mol) had less but similar docking as that of Sotrastaurin (À6.55 kcal/mol) against the PRKCD. However, when these docked complexes were subjected to MD simulation and MM/GBSA free binding energy calculation, it was found that the number of interactions and binding affinity of these natural compounds were too low compared with the positive control Sotrastaurin ( Figure S2). A similar case was observed for withanolides against PRKCE. The best-docked withanolides, namely 27-Hydroxywithanolide-B (À4.85 kcal/mol), Withanolide-P (À5.40 kcal/mol) and Withanolide-E (À5.15 kcal/ mol), had a similar docking score as Sotrastaurin (À5.25 kcal/ mol). However, when these docked compounds were subjected to 100 ns of MD simulation, it was found that Withanolide-E could not bind stably with PRKCE and came out of pocket within 30 ns of simulation. While other compounds could bind stably for 100 ns, the interaction and binding free energy indicated that they were not as potent as positive control Sotrastaurin ( Figure S3). In summary, the simulation and thermodynamic studies against PRKC isoforms suggested that even if these molecules could bind to the PRKCD and PRKCE, they are not potent and selective enough against these targets compared to already known inhibitors. After PRKC isoforms, AR was predicted to be one of the crucial targets in inverse virtual screening. Androgen's action depends on the AR, a liganddependent nuclear transcription factor. AR has been reported to be overexpressed in various cancers. Here, Metribolone was used as a positive control against AR, as it is an anabolic steroid that was reported to be effective against advanced breast cancer, but further research was abandoned on it due to hepatotoxicity (Matias et al., 2000). When the binding and interactions of the withanolides against AR were investigated, it was found that none of the top docked withanolides, namely, 17-Hydroxywithanone (À3.80 kcal/mol), withanolide-E (À2.76 kcal/ mol), and withanolide-A (À2.57 kcal/mol) had even half the docking score as compared to the already known inhibitor Metribolone (À12.58 kcal/mol). Therefore, these complexes were not subjected to further MD simulation and binding free energy analyses. The next predicted target was PTGS-2, which plays a primary role in synthesizing prostaglandins from the Arachidonic acids. The prostaglandins produced via PTGS-2 have a role in the inflammatory effects. The inhibitor of PTGS-2 effectively reduces the pain and inflammation, mainly in diseases like rheumatoid arthritis (Ricciotti & FitzGerald, 2011). Therefore, the binding and interactions of withanolides against PTGS-2 were investigated. It was found that none of the withanolides could bind at the ligand-binding site, possibly due to the smaller catalytic binding pocket of PTGS-2 ( Figure S4). Finally, the withanolides were studied for interaction with PDE4D. The docking results showed that the top docked withanolides had a better docking score than the positive control Zardaverine. MD simulations results showed that the binding interactions, binding-affinity, and dynamics were similar to the positive control Zardaverine throughout the simulations. Among all the top targets selected through inverse screening and network analysis, withanolides showed the best binding affinity and structural interactions with the active site residues of PDE4D only, compared to positive control. Hence, the MD simulation analysis of the withanolides against PDE4D has been described in detail in the following section.

Structure-based study shows the inhibitory potential of withanolides against PDE4D
Firstly, the crystal structure of PDE4D (PDB ID: 1MKD) bound with Zardaverine (chosen as positive control) was retrieved and prepared for the molecular docking and MD simulation studies. The bound inhibitor Zardaverine was re-docked to assess the best binding pose, interactions and validation of the docking algorithm. It was found that the docked pose complex was similar to the crystallized pose, and no difference was found in the RMSD (<0.1 Å) value aligned with the crystallized pose ( Figure S5). In the best docking pose, Zardaverine made two hydrogen bonds with Asn418 and Gln466 (Figure 2A). The docking score in the best binding pose was À8.26 kcal/mol. Further, to compare and investigate the binding interactions and energy, all the 20 screened withanolides were docked at the same active site. It was found that 4 out of the 20 withanolides had a higher docking score than Zardaverine. The 17-Hydroxywithaferine-A (À9.48 kcal/mol) had the highest docking score in the best binding pose, and it was making identical hydrogen bonding with Asn418 and Gln 466 similar to Zardaverine ( Figure 2B). 27-Hydroxywithaone (À9.41 kcal/mol) had the second-highest docking score and it was making a hydrogen bond with Gln466 in the best binding pose ( Figure 2C). These were followed by Withanolide-R (À8.61 kcal/mol) and Withaferin-A (À8.31 kcal/mol), and both were making a hydrogen bond with Gln466 in the best binding pose ( Figure  2D and 2E). All the polar and the non-polar interactions of the best-docked pose and docking scores are reported in Table 1 and their 2D representation in Figure S6.
From the observed interactions, it was evident that these four withanolides could be the inhibitor of PDE4D, and Gln466 was the crucial interactive residue that provided the binding affinity to these ligands at the catalytic site. However, from various previous studies, it is known that the docking score alone cannot be used to describe the binding stability of the ligand at the catalytic pocket. To validate the binding dynamics and the interactions of these ligands, a classical MD simulation of 100 ns was performed. After performing the MD simulation, the RMSD of the complexes was analysed, and it was found that all the complexes got stabilized within the 20 ns of the simulations. Although a bit of fluctuation was seen in the PDE4D-Zardaverine complex around the 50 ns, the fluctuation was less than 1 Å. Overall, the PDE4D-17-Hydroxywithaferin-A (2.18 ± 0.17 A) had the lowest average deviation throughout the simulations, followed by PDE4D-Withaferin-A (2.24 ± 0.17 A), PDE4D-Withanolide-R (2.47 ± 0.24 A), PDE4D-27-Hydroxywithanone (2.60 ± 0.37 A) and PDE4D-Zardaverine (2.72 ± 0.31 A) ( Figure  3A). The ligand RMSD with respect to protein backbone was then analysed, and it was found that all the ligands got stabilized within the first 50 ns of the simulation time, while major fluctuation was seen only in Withaferin-A ( Figure S7). Further, to investigate residue-wise fluctuation, RMSF was calculated and analyzed. No major fluctuation was found in any of the complexes ( Figure 3B). All the MD analysed data has been shown in Table 2.
As the global and local dynamics of the complexes could not be analysed easily, PCA was calculated and plotted in the two principal components using a kernel density plot. The PCA was calculated for all the simulated complexes on the last 60 ns of the simulation trajectory. In all the cases, a total of 10 Principal components were calculated; however first two components were sufficient to tell more than 60% of the motion of the complexes. Therefore, the first two principal components (PC1, PC2) were selected to analyze their projection of trajectories during the simulations. The kernel density plot showed that the PDE4D-27-Hydroxywithanone covered the highest phase space, followed by PDE4D-Zardaverine, PDE4D-Withanolide-R, PDE4D-Withaferin-A. PDE4D-17-Hydroxywithaferin-A had the most stable movement. However, there was no similarity in the direction of the movement in any of the complexes; and, when the average structures of PDE4D-Withanolides from each cluster were superimposed with PDE4D-Zardaverin, it was found that the change in RMSD was less than 1.5 Å, and no major change in the global structure was observed ( Figure S8). As it has been well reported and known that hydrogen bonding significantly improves the ligand-binding affinity towards the receptor, the average number of hydrogen bonds throughout the simulation was then calculated. The hydrogen bond numbers indicated that Withaferin-A (1.60 ± 0.60), 27-Hydroxywithanone (1.77 ± 0.58) and Zardaverine (1.83 ± 0.48) had more than 1 average hydrogen bonding per frame followed by 17-Hydroxywithaferin-A (0.92 ± 1.0) and Withanolide-R (0.14 ± 0.36) ( Figure 3C). To analyze the binding affinity of the studied compounds towards PDE4D, the MM/GBSA binding free energy was then calculated. The MM/ GBSA binding free energy showed that Withanolide-R (À60.51 ± 0.40 kcal/mol) had the highest binding affinity, which is mainly contributed by Van der Waal and lipophilic energy, followed by Zardaverine (À51.29 ± 0.29 kcal/mol), here apart from Van der Waal, and lipophilic, electrostatic energy contributed significantly in comparison to others, followed by Withaferin-A (À48.20 ± 3.11 kcal/mol) and 27-Hydroxywithanone (À55.50 ± 20.51 kcal/mol) ( Figure 3D). The individual average energy contribution due to Electrostatic, H-bond, lipophilic and Van der Waal interactions has been shown in Table S3. Finally, the simulation interactions were analyzed to assess the crucial residues that were making significant interactions (>30% of simulation time). The simulation interaction diagram showed that, in the case of Zardaverine, Thr430 and Gln446 were involved in the hydrogen bonding for more than 0.6 fractions of the simulation time, while Ile433 and Phe469 were the crucial residues  making hydrophobic contacts significantly ( Figure 4A). Similar was the case in all the withanolides, except Withanolide-R. In Withaferin-A, 17-Hydroxywithaferin-A and 27-hydroxywithanone, the residue Gln466 was significantly involved in hydrogen bonding, while the residue Phe469 was in the hydrophobic contact ( Figure 4C-E). However, in the case of the PDE4D-Withanolide-R complex, it was found that Pro453 and Ser465 were making the crucial interactions through hydrogen bonding and water bridge contacts, while Tyr472 was involved in hydrophobic interactions ( Figure 4B). Further, the total number of all the interactions throughout the simulation is shown in Figure S9.
The overall results and analysis showed that four withanolides had similar interactions as well as a binding affinity towards PDE4D like Zardaverine, and hence these lead compounds could be investigated further by in-vitro and in-vivo assays.

Conclusion
In the current study, the most probable targets of the withanolides were screened using ligand-based methods. Inverse screening results were further validated using molecular docking, classical MD simulations and MM/GBSA free energy calculations to find the best target and natural inhibitors. The overall analysis suggested that the PDE4D was the best therapeutic target of the withanolides. Withaferin-A, 17-Hydroxywithaferrin-A, 27-Hydroxywithanone and Withanolide-R were the top hit compounds with similar interactions and binding affinity towards the PDE4D like positive control Zardaverine. It was also observed that Gln466 and Phe469 provided stability to the potent inhibitors at the catalytic pocket through hydrogen bonding and hydrophobic contact. It is clear from the data presented in this study that the activity of PDE4D can be modulated by the natural bioactive compounds (withanolides) from Ashwagandha. The presented results warrant further in vitro and in vivo studies. However, the presented insights into the mechanism of action would be crucial in the novel drug design for PDE4D-related diseases.

Disclosure statement
The authors declare no conflict of interest.

Funding
The author(s) reported there is no funding associated with the work featured in this article.

Data availability statement
The authors confirm that the data supporting the findings of this study are available within the article and/or its supplementary materials.

Authors' contribution
A.R. contributed to experiments and manuscript writing. V.K. contributed to conception, experiments and manuscript writing. D.S. contributed to conception, design, manuscript writing and resources. All authors have read and agreed to the published version of the manuscript.