A Multi-Target mechanism of Withania somnifera bioactive compounds in autism spectrum disorder (ASD) Treatment: Network pharmacology, molecular docking, and molecular dynamics simulations studies

2.1 Identification and ADME evaluation of potential traditional Chinese medicine compounds

Information about the active compounds of W. somnifera was retrieved from IMPPAT (https://cb.imsc.res.in/imppat/) (Mohanraj et al., 2018), ChEBI (https://www.ebi.ac.uk/chebi/) (Hastings et al., 2016), KNApSAcK Core System database (http://www.knapsackfamily.com/KNApSAcK/) (Y. Nakamura et al., 2014); PCIDB (https://www.genome.jp/db/pcidb/), Dr. Duke's Phytochemical and Ethnobotanical Databases (https://phytochem.nal.usda.gov/phytochem/search) (Duke and Bogenschutz, 1994), and wide-scale review of the literature. The compounds were tested utilizing the ADME evaluation approach, with oral bioavailability (OB) and drug-likeness (DL) being the most essential parameters. For further investigation, bioactive components having an OB of ≥ 30 % and DL ≥ 0.18 were chosen (Xu et al., 2012); (Sadaqat et al., 2023). PubChem (https://pubchem.ncbi.nlm.nih.gov/) (Kim et al., 2016) along with ChemSpider (http://www.chemspider.com/) (Pence & Williams, 2010) were also utilized for information collection on molecular weight and chemical ID (Kaur et al., 2023).

2.2 Prediction and screening of target genes associated with compounds and disease

The STITCH (http://stitch.embl.de/) (Kuhn et al., 2007) and Swiss Target Prediction (http://www.swisstargetprediction.ch/) (Daina et al., 2019) databases were utilized to determine the targets associated with each active chemical. Both databases used “Homo sapiens” as a reference organism. STITCH and Swiss Target Prediction databases were employed to choose targets having a combined score/probability score of ≥ 0.7 (Fatima et al., 2024). A compound-target network was constructed in Cytoscape using these compounds and corresponding targets as input.

The keyword “autism spectrum disorder” was searched to retrieve disease-related targets through the GeneCards database (https://www.genecards.org/) (Safran et al., 2010). Subsequently, a Venn diagram was created using the online tool Evenn (http://www.ehbio.com/test/venn/#/) (Chen et al., 2021) to illustrate the overlapping targets between compound and disease, facilitating the identification of common targets. The resultant common targets were assumed as key targets and processed for further analysis.

2.3 GO enrichment and KEGG Pathway analysis

The DAVID (https://david.ncifcrf.gov/) database was used for gene ontology (GO) analysis that comprised three categories (Dennis et al., 2003). These categories are biological process (BP), cellular component (CC), and molecular function (MF) (Fatima et al., 2024). By lowering the likelihood score to less than 0.05, the top ten GO terms and KEGG pathways associated with ASD were chosen for more in-depth analysis. A bubble chart was constructed using the online program SRplot (https://www.bioinformatics.com.cn/en) to analyze the related functional pathways of the targeted genes.

2.4 Protein-Protein interaction (PPI) analysis

The STRING database (https://string-db.org/) was utilized for the PPI analysis and “Homo sapiens” was taken as the reference (Mering et al., 2003). The output network in TSV format was downloaded and a parameter of combined score ≥ 0.5 was employed (Alshehri & Alshehri, 2023). The network was then loaded into Cytoscape v.3.9.1 where the CytoHubba plugin (Chin et al., 2014) was used. CytoHubba's MCC, MNC, degree, betweenness, and closeness scoring algorithms were utilized to predict hub genes. Common genes from the above algorithms were considered as main genes and processed for further analysis.

2.5 Network construction

To comprehend the mechanism of W. somnifera in ASD, network analysis was carried out. The network diagram was created and displayed using the Cytoscape 3.9.1 program (Kohl et al., 2011). The network's nodes represented bioactive components and target genes, on the other hand edges were used to show the association between the compounds as well as targets. The Cytoscape Network Analyzer (Saito et al., 2012) was then utilized for the calculation of the degree, a topological measure that shows the relevance of a component/target/pathway in the network.

2.6 Analysis of molecular docking

The interaction of the hub gene with bioactive substances was predicted using molecular docking (A. Fatima, Arora, et al., 2023); (A. Fatima, Khanum, et al., 2023). The PDB database (https://www.rcsb.org/) was employed to retrieve and download the 3D structure of the hub gene in PDB format (Kouranov et al., 2006). The PubChem database (https://pubchem.ncbi.nlm.nih.gov/) was utilized to collect the 3D structure of bioactive chemicals in SDF format (Kim et al., 2019); (Jindal et al., 2023). The SDF format was then transferred to the PDBQT format through Open Babel software (O’Boyle et al., 2011). To remove the water molecules and natural ligand groups from the protein structure UCSF Chimera (Pettersen et al., 2004) was employed. The CASTp tool (https://sts.bioe.uic.edu/castp/) was used to forecast binding pockets (Tian et al., 2018). Additional molecular docking investigations were performed via the PyRx software interface using AutoDock Vina (Trott & Olson, 2010). In PyRx the docking grid for IL6 was set at coordinates X: 12.1066, Y: –32.4441, and Z: 0.4771 with dimensions X: 21.1750, Y: 22.1155, and Z: 18.1749 Å from the center, and an exhaustiveness level set to 8. The best-docked complex was chosen with the best postures, higher binding energy, and lowest root mean square deviation (RMSD). The interactions between the top three docking complexes with the higher binding energy were displayed in 2D by utilizing Discovery Studio Visualizer (Headquarters, 2008) and in 3D using UCSF ChimeraX (Pettersen et al., 2021).

2.7 Analysis of molecular dynamic simulation

The AMBER22 program was used to carry out a molecular dynamic simulation (Harris et al., 2022); (Ahmad, 2023). A molecular dynamic simulation (MDs) is an in-silico simulation technique mostly utilized for macromolecule and docked atom analysis (Wani et al., 2021). A macromolecule that is allowed to show dynamic behavior for a specific amount of time in a molecular dynamic simulation pipeline and the atom and molecular trajectories are originated by solving Newton's equations of motion (Gupta, 1992). The AMBER software was used in this research to determine the dynamic behavior of the compounds in complex with the receptor (Mortier et al., 2015). The purpose of this was to read how the compound affinity for the receptor enzyme with time. The complexes were processed via the Antechamber program. In the case of receptor processing the FF19Sb force field was applied while for ligand processing the GAFF2 force field was used. An appropriate quantity of counter ions was introduced into the system to achieve charge neutrality (Guterres et al., 2022). A cubic box size of eight angstroms was supposed to be enough to solvate the complexes and the steepest descent and conjugate gradients algorithms were applied (Nicholls & Honig, 1991). The heating of the complexes was achieved gradually till 310 K followed by system equilibration and a production run of 100 ns. The trajectory analyses were performed using CCPTRAJ (Ahmad et al., 2020).

2.8 MMPB/GBSA analysis

The binding free energies of docked ligands with the enzyme were calculated using MMPB/GBSA analysis (Zaman et al., 2021). This was achieved using a script called MMPBSA.py which is part of the AMBER v22 program (Miller et al., 2012). The script assumed 5000 frames from the trajectories selected at regular intervals. The formula for MMPB/GBSA energy: $$\Delta Gbind=\Delta Gcomplex-\Delta Greceptor-\Delta Gligand$$

The terms in the equation are ΔGligand for free ligand energy, ΔGreceptor for the free energy of the receptor protein, ΔGcomplex for complex free energy, and ΔGbind for total binding free energy. The following formula will be used to find the unique free energies of a complex, protein, and ligand. It is found that the performance of the MM/PBSA and MM/GBSA approaches is similar. The MM/PBSA uses the Poisson-Boltzmann equation to determine the electrostatic energy contribution to the free energy, while the MM/GBSA uses the Generalized Born equation which is believed to be a faster solution for the earlier equation.

3.1 Identification and selection of bioactive compounds of w. Somnifera

Collectively 242 bioactive constituents of W. somnifera have been collected and retrieved from various databases and a review of the literature (Supplementary Table S1). After the screening and removal of duplicates, 13 putative compounds ((+)-Catechin, 24-Methyldesmosterol, Beta-Sitosterol, Campesterol, Fucosterol, Kaempferol, Oleanolic acid, Quercetin, Stigmasterol, Stigmasterone, Withaferin A, Withanolide J, and Withanone) which retained the criteria of DL ≥ 0.18 and OB ≥ 30 % were selected. The details of compounds, their PubChem ID (CID), molecular formula (MF), molecular weight (MW), OB, DL, and 2D structure were enlisted in Table 1.

Table 1 Compounds, their characteristics, and 2D structure.

Compound Name	MF	MW	DL	OB	2D Structure	CID
(+)-Catechin	C15H14O6	290.27	0.64	0.55		9064
24-Methyldesmosterol	C28H46O	398.7	0.76	0.55		193,567
Beta-Sitosterol	C29H50O	414.7	0.78	0.55		222,284
Campesterol	C28H48O	400.7	0.59	0.55		173,183
Fucosterol	C29H48O	412.7	0.85	0.55		5,281,328
Kaempferol	C15H10O6	286.24	0.5	0.55		5,280,863
Oleanolic acid	C30H48O3	456.7	0.37	0.85		10,494
Quercetin	C15H10O7	302.23	0.52	0.55		5,280,343
Stigmasterol	C29H48O	412.7	0.62	0.55		5,280,794
Stigmasterone	C29H46O	410.7	0.5	0.55		14,807,783
Withaferin A	C28H38O6	470.6	0.37	0.55		265,237
Withanolide J	C28H38O6	470.6	0.46	0.55		21,679,022
Withanone	C28H38O6	470.6	0.45	0.55		21,679,027

3.2 Identification and screening of potential targets

From 13 active compounds, a total of 134 putative potential targets were retrieved through STITCH and Swiss Target Prediction databases (Fig. 1A; Supplementary Table S2). On the other hand, a total of 8172 targets linked with autism spectrum disorder were obtained through the GeneCards database (Supplementary Table S3). A Venn diagram was generated to determine the shared targets from compound-related along ASD-related targets (Fig. 1B). As a result, 70 potential anti-ASD genes from W. somnifera were chosen and regarded as key targets (Supplementary Table S4).

Close

Fig. 1

(A) Compound-target network, (B) A Venn diagram drawn between compound and disease-related targets to find common target genes.

Export to PPT

3.3 GO and KEGG signaling Pathway analyses

To comprehend the possible biological function of W. somnifera targets, gene ontology (GO) at three biological levels: biological processes (BP), cellular components (CC), and molecular function (MF) have been done and a total of 234 BP, 34 CC, 68 MF, and 106 KEGG pathways terms were significantly found to be enriched. Additionally, the analysis of KEGG pathways was also carried out. According to BP, most of the genes were involved in the apoptotic process of negative regulation, peptidyl-serine phosphorylation, xenobiotic stimulus–response, and protein autophosphorylation (Fig. 2A). CC indicates that these genes are predominantly present in cytosol, membrane raft, extracellular exosome, and cell membrane (Fig. 2B). MF revealed that these genes are associated with the binding of ATP, DNA, and protein, protein serine/threonine kinase activity, and protein tyrosine kinase activity (Fig. 2C). Furthermore, KEGG pathways analysis pointed toward significant involvement in HIF-1 signaling, MAPK signaling as well and TNF signaling pathways (Fig. 2D).

Close

Fig. 2

GO functional enrichment and KEGG pathways analysis, (A) Biological processes, (B) Cellular components, (C) Molecular functions, and (D) KEGG pathways.

Export to PPT

3.4 Protein-Protein interaction (PPI)

Interactions of proteins with other proteins were determined through the STING database. MCC, MNC, degree, betweenness, and closeness centralities were used. A Venn diagram was generated to predict common genes from all these five methods (Supplementary Table S5; Fig. 3A). As a result, 6 common genes were analyzed. These are IL6, CASP3, AKT1, TNF, SRC, and EGFR (Fig. 3B). In comparison with KEGG pathways one target gene IL6 was utilized for the analysis of molecular docking.

Close

Fig. 3

(A) Calculation of hub genes by MCC, MNC, Degree, Closeness, and Betweenness scoring algorithms of CytoHubba, (B) Venn diagram showing the common genes of 5 different scoring algorithms.

Export to PPT

3.5 Compound-target-pathways network

A network of compound-target pathways was developed with the help of Cytoscape v.3.9.1. Active compounds, hub genes, and pathways were represented by nodes. These nodes are connected through edges (Fig. 4). There were 25 nodes and 58 edges in the network. QSPR models utilizing graph theoretical approaches are vital for predicting molecular properties based on topological indices, offering insights into structure–property relationships efficiently and comprehensively.

Close

Fig. 4

A network of compound-target-pathways. Arrow symbols represent active constituents, circle symbols represent hub genes, and square symbols represent the pathways related to hub genes.

Export to PPT

According to the network ASD is mediated by MAPK, PI3K-Akt, TNF, ErbB, C-type lectin receptor, FoxO, Relaxin, Rap1, IL-17, Toll-like receptor, HIF-1, and JAK-STAT signaling pathways (Table 2).

3.6 Molecular docking

All 13 active constituents were utilized for molecular docking analysis in which some compounds were plant-specific compounds (Withaferin A, Withanolide J, and Withanone). Based on PPI and KEGG pathways analysis one hub gene named IL6 was employed for the analysis of molecular docking. Protein structure (IL6 (PDB ID: 1ALU)) was gained through the Protein Data Bank (PDB). UCSF Chimera v1.16 was used to improve structure, remove non-standard amino acids, and minimize energy use. The possible targets of substances that might reduce the occurrence of ASD were screened out utilizing molecular docking. The binding energy and docking scores were employed as major criteria for the compound selection (Table 3). Withanolide J was stably bound with 1ALU protein with a binding affinity of −6.9 kcal/mol and RMSD of 1.03 Å. Withanolide J side chains formed hydrogen bonds with LEU A:92, VAL A:96, LYS A:120, PRO A:141, ALA A:145, and LEU A:148. Withanone and Withaferin A have the highest binding affinity −6.8 kcal/mol and −6 kcal/mol with RMSD of 0.157 Å and 0.536 Å respectively (Fig. 5). Withanone and Withaferin A made strong bonding with LYS A:120, ILE A:123, PRO A:141, ASN A:144 and GLN A:127, PRO A:141, ASN A:144 residues. All active compounds show strong hydrogen bonding with PRO A:141 except Oleanolic acid, Stigmasterone, and Quercetin. These results suggest that the IL6 target protein is efficiently bound by active components of W. somnifera, which then interact as ASD repressors.

Close

Fig. 5

IL6 docking complexes with their three highest binding affinity compounds.

Export to PPT

3.7 Molecular dynamic simulation

Essentially, the dynamic behavior of macromolecules is confirmed by a molecular dynamic simulation study. The simulations analysis includes root mean square deviation (RMSD) (Sargsyan et al., 2017), radius of gyration (RoG) (Lobanov et al., 2008), and root mean square fluctuation (RMSF) (Pitera, 2014). The foundation for all of these studies was the complexes' carbon alpha atom. The purpose of these investigations was to ascertain whether the ligand binding to the receptors was stable and whether the interactions continued throughout the simulation time. Stable receptor-ligand contact will guarantee that the ligand is correctly presented to the IL6 target protein. The RMSD plot of the systems was uniform at first indicating no discernible structural changes. In Fig. 6A, the mean values of IL6_Withanolide J, and IL6_Withanone were found to be (2.00 Å) and (2.21 Å) respectively, whereas the highest values of the systems' root mean square deviation (RMSD) ranged up to 3 Å. Second, the RMSF was calculated to reveal information on the flexibility of the receptor residues in the presence of the ligand molecule (Fig. 6B). Nearly all system residues were within the region of average stability (<5 Å). The RMSF of the systems for IL6_Withanolide J, and IL6_Withanone showed maximum values of (3.69 Å and 7.25 Å), mean values of (1.01 Å and 1.06 Å), and minimum values of (0.45 Å and 0.44 Å). The higher degree of flexibility seen in some of the residues was shown to be a result of the loop pressure inside the system, the compactness of the system was analyzed against time using the RoG analysis. The way ligands attached to the receptors was unaffected by these variations, though. With IL6_Withanolide J and IL6_Withanone having RoG maximum values of (16.40 Å and 16.85 Å), mean values of (16.11 Å, and 16.34 Å), and minimum values of (15.86 Å and 15.88 Å) (Fig. 6C). RMSD demonstrated that all the systems were compact and did not experience any notable alterations after the simulation. The Beta Factor study yielded the following results: IL6_Withanolide J and IL6_Withanone had mean values of (39.77 Å and 50.09 Å), lowest values of (5.46 Å and 5.13 Å), and maximum values of (358.89 Å and 1385.31 Å) (Fig. 6D).

Close

Fig. 6

Different simulation based statistical assays. (A) RMSD, (B) RMSF, (C) RoG, and (D) Beta Factor plots.

Export to PPT

To learn more about the surface area of the IL6 that interacts with the solvent molecules, solvent-accessible accessible surface area (area) analysis was performed for the ligands. IL6_Withanolide J (9526.49 nm2), and IL6_Withanone (9776.63 nm²) are the average values for the systems. According to (Fig. 7), the lowest values of SASA for IL6_Withanolide J, and IL6_Withanone were (8364.2 nm² and 8891.45 nm²). The maximum values were (10730.5 nm² and 10795.4 nm²), respectively. Plots show the significant variations observed upon ligand binding.

Close

Fig. 7

showing SASA plot for IL6_Withanolide and IL6_Withanone.

Export to PPT

3.8 MMPB/GBSA analysis

The MMPB/GBSA analysis was completed for the selected docked complexes. These methods are believed to be more effective in determining the docked ligand binding affinity with the receptor protein. The very negative net binding energies of all docked and control complexes demonstrated the formation of strong intermolecular systems and stable complexes. It was discovered that the van der Waals force, which keeps ligand docking at the docked site and stabilizes systems, is the most potent force produced by the complexes. The net van der Waals energies of the IL6_Withanolide J, and IL6_Withanone were −36.15 kcal/mol and −13.64 kcal/mol, respectively. Furthermore, it was discovered that the electrostatic energy of every docked complex was incredibly constant. The solvation energy was the least favorable of the calculated energies and contributed negatively to the net energy. In MM-GBSA, the net solvation energy of the IL6_Withanolide J, and IL6_Withanone was 8.20 kcal/mol and 7.68 kcal/mol, respectively. In contrast, in MM/PBSA, the net solvation energies of the IL6_Withanolide J, and IL6_Withanone were 7.54 kcal/mol, and 8.01 kcal/mol, respectively. Additional details on the energy terms and values are provided in Table 4.