Interleukin-10 as Covid-19 biomarker targeting KSK and its analogues: Integrated network pharmacology

COVID-19 caused by the SARS-CoV-2 virus is widespread in all regions, and it disturbs host immune system functioning leading to extreme inflammatory reaction and hyperactivation of the immune response. Kabasura Kudineer (KSK) is preventive medicine against viral infections and a potent immune booster for inflammation-related diseases. We hypothesize that KSK and KSK similar plant compounds, might prevent or control the COVID-19 infection in the human body. 1,207 KSK and KSK similar compounds were listed and screened via the Swiss ADME tool and PAINS Remover; 303 compounds were filtered including active and similar drug compounds. The targets were retrieved from similar drugs of the active compounds of KSK. Finally, 573 genes were listed after several screening steps. Next, network analysis was performed to finalize the potential target gene: construction of protein-protein interaction of 573 genes using STRING, identifying top hub genes in Cytoscape plug-ins (MCODE and cytoHubba). These ten hub genes play a crucial role in the inflammatory response. Target-miRNA interaction was also constructed using the miRNet tool to interpret miRNAs of the target genes and their functions. Functional annotation was done via DAVID to gain a complete insight into the mechanism of the enriched pathways and other diseases related to the given target genes. In Molecular Docking analysis, IL10 attained top rank in Target-miRNA interaction and also the gene formed prominent exchanges with an excellent binding score (> = -8.0) against 19 compounds. Among them, Guggulsterone has an acute affinity score of -8.8 for IL10 and exhibits anti-inflammatory and immunomodulatory properties. Molecular Dynamics simulation study also performed for IL10 and the interacting ligand compounds using GROMACS. Finally, Guggulsterone will be recommended to enhance immunity against several inflammatory diseases, including COVID19.

miRNAs of those hub genes were searched through miRNet, a web-based platform for functional identification. MicroRNAs play an essential role in manipulating gene expression [12]. Enrichment analysis is also more beneficial for finding the Gene Ontology, pathway and other diseases related to the target hub genes. A molecular docking study using AutoDock Vina was implemented by calculating the docking score to verify the interaction between compounds and targets. This analysis also detailed how these novel drugs exert their therapeutic activity. Compound-Target-Pathway Network was constructed in Cytoscape to highlight the potential mechanism of action of the compounds. Finally, MD simulation for IL10 with the ligand molecules Herkinorin, Guggulsterone and mulberrofuran W selected based on the molecular docking results.

Collection and screening of compounds
KSK concoction is composed of 15 plant ingredients. Phytoconstituents of these 15 plants were retrieved from literature reading and IMPPAT online database. IMPPAT is a technological platform for gathering information on traditional Indian medicinal plants and their phytochemical compounds [13]. Also, similar combinations of phytochemical compounds were recouped from TTD based on the high-similarity score of > = 85% [14,15]. Duplicates, repeats, invalids, and no PubChem entry compounds were physically removed from the total compounds collection. Next, an important step is screening to identify the valid compounds based on the given criteria: Lipinski violations, PAINS alerts, Bioavailability score, Catechol A, Solubility class and Synthetic accessibility. PAINS Remover server and Swiss ADME tool were used for the screening process. Lastly, valid compounds were collected for further analysis [16].

Finding target genes
Target genes were retrieved in two ways: a). targets of similar compounds from TTD, and b). targets of COVID-19 from the NCBI Gene database (https://www.ncbi.nlm.nih.gov/gene). Both similar compound genes and COVID-19 genes were collected for network analysis and functional annotation.

Construction of protein-protein interaction network
The online STRING database was used for the construction of the PPI network. Listed target genes were uploaded, and Homo sapiens was nominated as a standard organism [17]. Nodes represent the genes, and edges represent the interaction between the nodes.

Cluster analysis and detecting hub genes
Cytoscape is a free, open-source platform to visualize complex networks, molecular interactions, and several plug-ins to perform diverse analyses [18]. MCODE is a Cytoscape plug-in to find highly interconnected nodes formed as clusters in a network. Top hub genes were identified in clusters. Hub genes were ranked via the MCC algorithm in cytoHubba [19]. miRNet is also used for exploring the functional annotation of genes; miRNA Family with their functions, diseases they are related to, associated clusters and transcription factors. A complex network has been simplified with the help of the Degree filter.

Key enrichment and pathway analyses
Functional Enrichment analysis was done via DAVID to understand the biological information of the genes [21]. GO terms: Biological Process, Cellular Components, Molecular

Integrated network pharmacology
Functions; Pathways Analysis: KEGG Pathway, Reactome; and Diseases related to the given target genes information were retrieved. Pathway and disease relation expression was shown in a bubble plot using ImageGP (http://www.ehbio.com/ImageGP/). GO enrichments results were established as a 2D bar graph using an MS Excel file.

Molecular docking analysis
The three-dimensional crystal structures for the target hub genes were downloaded from PDB [22]. The active sites of the protein structures were predicted by the CastP server (http://sts. bioe.uic.edu/). The missing residues in the target proteins were added using SwissPDB Viewer [23]. Water molecules were removed, and Kollman Charge hydrogen atoms were added for protein preparation in AutoDock [24,25]. 3D structure of ligand molecules was downloaded from the PubChem database. Some ligands do not have a 3D structure in the PubChem database; their 2D files were downloaded and converted to the 3D structure using Chem Sketch [26]. MMFF94 force field parameter was used to optimize the compounds using OpenBabel [27] in Ubuntu. Molecular Docking analysis was done using the AutoDock Vina tool [28]. Protein-ligand 3D interactions were analyzed using PLIP software [29]. 3D and 2D hydrogen bond interactions of the protein-ligand complexes were visualized in PyMOL [30] (DeLano, 2002) and LigPlot+ [31]. The main aim of this molecular docking is to predict the interaction of the protein and ligand molecules and how they function as a complex. The result can be shown mainly based on the binding affinity of the protein-ligand complex. More negative the acute affinity score, the higher the critical relationship between the protein-ligand complexes.

Compound-target interaction network construction
After the Molecular Docking analysis was done, the interaction between the target proteins, ligand molecules, and significant pathways was predicted. This interaction was drawn as a Compound-Target-Pathway network in the Cytoscape tool. Network Analyzer is a Cytoscape plug-in used to compute the topological analysis for the constructed network. It gives information on the number of nodes and edges, neighbors, network diameter, radius, density, etc. Lastly, the calculated result was shown as summary table and node-specific tables [32].

Molecular dynamics simulation
The three dimensional protein structure of IL10 was downloaded from PDB and visualized in PyMOL. Molecular docking results for IL10 showed interactions with top four ligands: Herkinorin, (R,S)-homoaromaline hydrochloride, Guggulsterone and mulberrofuran W. Further MD simulation was performed for the protein-ligand complexes of IL10 and the above four ligands. MD simulation study was performed with GROMACS 2022.2 installed in Ubuntu on Windows. GROMACS is a software suite for molecular dynamics simulation [33]. The missing atoms were checked and filled using energy minimization with Spdb Viewer. Protein topology was created with 'charmm36-jul2021' force field. 'pdb2gmx' in GROMACS reads the.pdb file as input and select the specific force field depending upon the type of working system and add H atoms. It also generates GROMACS coordinate file and topology file. The force field needs to be compatible with the system. CGenFF online server was used for generating ligand topology files [34]. From both protein and ligand topologies,.gro files were combined for creating protein-ligand complex topology and saved as new file. Cubic grid box was generated and solvent system was added. This system was neutralized by adding SOD (positive) and CLA (negative). Energy minimization is a vital step to remove the steric clashes and improper geometries in the system. Steepest descent algorithm is used for the energy minimization. There are two different phases in equilibrating the system: (i) nvt-defines the constant number of particle volume and temperature (to stabilize the temperature); (ii) npt-defines the particle pressure and temperature of all constants. Final step is to run the simulation for the system using 'mdrun'. The runtime for the simulation was 100ns for the system. There are major analyses to be done after MD simulation: (i) RMSD calculation-used to measure the scalar distance for the protein-ligand throughout the trajectory; (ii) RMSF calculation-used for characterizing local fluctuation in the model.

Systematic extraction of active compounds
There are 321 active compounds present in the 15 plant ingredients of KSK. These 321 active compounds were collected from the literature and IMPPAT database. PubChem CIDs were retrieved for 1245 KSK compounds and KSK similar compounds [35] (S1 File). 1, 207 compounds were selected after removing the duplications. Table 1A clearly explains the selection criteria for compounds.
Total compounds were screened to find the appropriate compounds to be taken for further steps. According to the GHS classification report, compounds with Catechol A property were removed because of their toxic nature [36] in the PAINS Remover server [37]. SwissADME was used to screen the compounds based on Lipinski's violations (above 1 violation compounds excluded), PAINS alerts (0 alert compounds), Bioavailability score (0.55 minimum) and synthetic accessibility (1 to 7). Every drug molecule should obey the Lipinski rule and Bioavailability score [38]. Here, 303 compounds were obeyed out of 1,207 compounds (S2 File).

Target genes collection
False positives were removed and 573 unique genes were established from the whole genes (S3 File). A total number of collected target genes was highlighted in Table 1B.

Network analysis of target genes
The PPI network was constructed for the obtained target genes in the STRING database ( Fig  2). Overall, 11,115 edges with 573 nodes interacted with each other. The constructed PPI network was visualized in Cytoscape 3.8.2 version. Cluster analysis was done in Cytoscape using the MCODE plug-in. 16 clusters were revealed from the STRING network (Table 2). After cluster generation, top hub genes were detected, showing highly interconnected regions. Hence, they were considered the significant nodes done by the cytoHubba plug-in in Cytoscape. TNF, IL6, IL10, IL1B, IL4, IFNG, IL17A, CSF2, CD4 and CXCL8 were the top 10 hub genes. These genes were ranked by the MCC method. MCC method is the most effective  Table 2. MCODE result shows 16 clusters score, number of nodes and edges, and node IDs.

Functional annotation of the target genes
To gain further insight into the functions of target genes, functional annotation provides complete information on Gene Ontology Enrichment analysis (which includes Biological Processes, Molecular Functions and Cellular Components), Pathway analysis (KEGG pathway, Reactome) and Other Diseases (DisGeNET) that target genes are related with, and this can be achieved through DAVID. A cut-off criterion for a p-value of <0.05 was taken to indicate the significance. The top significantly enriched GO terms regarding Gene Ontology analysis are shown in Fig 5 and S4 File.

Analyzing protein-ligand interaction by molecular docking
The top 10 hub genes were taken for docking analysis further. The 3-dimensional crystal structures of the hub genes were downloaded using the PDB database. The hub genes and their respective PDB IDs are as follows; TNF (1TNF), IL6 (1ALU), IL10 (1ILK), IL1B (1ILB), IL4 (1BBN), IFNG (1FG9), IL17A (2VXS), CSF2 (1CSG), CXCL8 (1IKL) and CD4 (1CDH). Proteins were prepared using AutoDock software. The missing residues were filled by the process of energy minimization using Swiss-PdbViewer (SPDBV). The active site residues were predicted using the CASTp server. The ligands were prepared using Open Babel. Finally, docking was done with 303 compounds in AutoDock Vina. The binding affinity scores having more negative values were considered a good score, and the protein-ligand complex has a better binding force. Finally, 9 compounds having a better binding affinity score (above -8.7) were taken for interaction analysis. Also, the complexes with zero RMSD value (Distance from best mode) were only taken (Table 4). Non-covalent interaction and Hydrogen bonding were analyzed in PLIP software (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index) and PyMOL software respectively. Non-bonded interactions of Protein-ligand complex were highlighted in Fig 7. 3D and 2D Hydrogen bond interactions were visualized in PyMOL and LigPlot+, respectively.
a. Hydrophobic interactions. Hydrophobic interactions are the strongest among the noncovalent exchanges, and they exhibit non-polar nature unlike hydrophilic. The top 9 compounds formed hydrophobic interactions with the target showing higher affinity in this study. In this part, Phenylalanine (PHE) residue comes under hydrophobic classification, and at times, aromatic side-chain might help for stacking interactions; six out of nine compounds formed hydrophobic interaction with PHE: Tubocurarine and Homoaromoline in IFNG_A, Guggulsterone, (R, S)-homoaromaline hydrochloride, Herkinorin and mulberrofuran W in IL10. Tyrosine residue, which implies polar characteristics, is made from phenylalanine; seven out of nine compounds formed hydrophobic interaction with TYR: Tubocurarine, Homoaromoline, Triptocalline A and 24-hydroxyursolic acid in IFNG_A, Bis(6-hydroxybenzo[b]furan-2-yl) methanone in IL4, Guggulsterone and mulberrofuran W in IL10. Alanine has a nonpolar characteristic which comes under minor hydrophobic classification and can get buried to form many types of interactions; five out of nine compounds formed hydrophobic interaction with ALA residue: Tubocurarine, Homoaromoline and Triptocalline A in IFNG_A, (R, S)-homoaromaline hydrochloride and mulberrofuran W in IL10. Isoleucine contains two non-H atoms linked with c-beta carbon leading to the form of rigid protein structure, and they also come under hydrophobic classification; four out of nine compounds formed hydrophobic interaction with ILE residue: Tubocurarine, Triptocalline A and 24-hydroxyursolic acid in IFNG_A, Guggulsterone in IL10. Both Leucine and Valine are very hydrophobic and have six compounds (Tubocurarine and Triptocalline A in IFNG_A and Guggulsterone, (R, S)-homoaromaline hydrochloride, Herkinorin and mulberrofuran W in IL10), and three compounds ((R, S)-homoaromaline hydrochloride, Herkinorin and mulberrofuran W in IL10) interacted with individual residues in the respective target. Lysine (Lys) is partially hydrophobic, and it interacts with two compounds out of nine: Homoaromoline in IFNG_A and Bis  GLN1, ASP2, PRO3, TYR4, LYS6, GLU7, GLU9, ASN10, LYS12, LYS13, ASN16,  GLY18, THR27, LEU30, GLY31, LEU33, LYS34, TRP36, LYS37, GLU39, ARG42,  LYS43, GLN46, LYS68, GLU71, THR72, GLU75, ASP76, VAL79, LYS80, PHE81 c. Hydrogen bond interactions. Hydrogen bonding is essential in facilitating the molecular interaction between molecules, potentiating diverse cellular functions [39]. Hydrogen bonds form between two molecules by HBond Donor and HBond Acceptor, thus stabilizing the interaction between the complexes formed. In the case of Lysine, Homoaromoline and 24-hydroxyursolic acid created HBond with LYS residue of IFNG_A and Bis(6-hydroxybenzo [b]furan-2-yl) methanone formed HBond with LYS residue of IL4; In case of Aspartic acid, 24-hydroxyursolic acid formed HBond with ASP residue of IFNG_A; In case of Cysteine, Bis (6-hydroxybenzo[b]furan-2-yl) methanone formed HBond with CYS residue of IFNG_A [40].
d. Pi-cation interactions. This is one of the most robust interactions among the non-covalent interactions. Arginine is aromatic and can form pi-cation interaction with small molecules. The compound Bis(6-hydroxybenzo[b]furan-2-yl)methanone created pi-cation interaction with ARG residue of IL4.

Interaction analysis for compound-target network
Network construction was done for top hub genes with their pathways and top interacting compounds. Fig 8 shows the interaction between them based on the data obtained from the above results.
Also, the topological analysis for the given network was done using Network Analyzer. It computes the topological parameters of the network. Two different results were retrieved as Summary statistics and Node Specific Statistics. The Summary Statistics table has been given in S6 File and Node Specific Statistics has been given in Table 5. IL10 gene has a high degree value of 21 with betweenness and closeness centrality.

MD Simulation
For the given input PDB file the command 'pdb2gmx' select the specific force field (charmm36) and fills the missing atoms. It generated three different format files automatically: (1).gro file-It consists of defined force fields; (2).top file-It consists of non-bonded and bonded parameters; (3).itp file-It contains information of restrain position of heavy atoms. The generated '.top' file will be updated in further simulation steps. After protein topology, four ligands topology files were created using CGenFF server: (1).itp file-Ligand topology file; (2).prm file-Ligand dihedral parameters file; (3) ini.pdb file-Ligand structure file; (4).top file. For both protein and ligand topologies,.gro files will be generated as output and these files were combined for protein-ligand complex topology generation. In cubic box generation, Protein will be centered in cubic box, that is, after running the respective command the protein has been embedded in the box. The water molecules (solvent system) were added. To neutralize the system the model was checked whether it was positively or negatively charged. SOD atoms and CLA atoms were selected based upon the type of charge the system has been distributed in protein. The command 'neutral' will automatically recognize and neutralize the molecule by replacing its H 2 O molecule. After completely neutralizing the system, the input file 'minim.mdp' was given and energy minimization was done using 'em' command. Initially   Fig 8. Compound-target-pathway network constructed in cytoscape. *Ellipse-shaped nodes represent ten hub genes. Rectangle-shaped nodes represent the pathways involved by the genes. Triangle-shaped nodes represent the top nine compounds, and the dotted line represents the interaction between the nine compounds with the genes IFNG, IL10 and IL4.
https://doi.org/10.1371/journal.pone.0282263.g008 equilibration was done to the whole system before going to actual simulation step; else it will disintegrate the system and may give false results. The input file 'nvt.mdp' was given and temperature was stabilized (for 100ps) using 'nvt' command. The input file 'npt.mdp' was given to achieve complete pressure constant using 'npt' command. Finally, MD simulation was done by changing runtime (nsteps = 100ns) in 'md.mdp' file.
Calculating the RMSD of backbone atoms is the widely used method for checking the stability of the protein with respect to the initial structure. Generally RMSD values used to compare the docked conformation with the referenced conformation where RMSD value between 1-3Å is considered.
In this study, the interactions of ligands with IL10 were taken for MD simulation. Out of four ligands interacted with IL10, for (R,S)-homoaromaline hydrochloride compound had no scientific evidences regarding its functions, hence it is not considered for MD simulation. explains the RMSD results for the remaining three compounds (Guggulsterone, Herkinorin and mulberrofuran W) and the RMSD of backbone atoms of IL10. From the results obtained, the RMSD value reached 10Å for backbone conformation; the RMSD value is less than 1Å for ligand conformation. The RMSD value has very large deviation in the beginning in all three backbone conformation. In interaction with Guggulsterone, after certain period of time (from 50 ns) the backbone RMSD reaches the equilibrate state and fluctuations. In interaction with Herkinorin and mulberrofuran W, backbone RMSD had not reached equilibration state and fluctuation and was in constant dynamic state. Also the RMSD results of ligands shows that, Guggulsterone is considerably stable with respect to IL10 structure when compared to Herkinorin and mulberrofuran W stability. Despite of the above results, the RMSD of the backbone atoms for all the three results showed major deviation in the protein. Hence the protein showed major conformational change in the structure and the ligands might not be stable with the protein binding pocket. The local fluctuations in the protein RMSF results were shown in the Fig  10. The most fluctuated regions of the protein during the simulation were indicated in the resulted plots for each of the ligand molecules. The fluctuations were observed in both the ends of the protein than in the middle regions.
From the RMSF results, the interactions with all the three ligand molecules showed more fluctuations in the protein. These results showed that the protein secondary structure elements might be less rigid in complex state with the respective ligand molecules.
IFNG-A plays a significant role in anti-viral, anti-microbial and anti-tumor activity in humans. Typically, T-cells or Nature Killer cells produce the type II interferon for activating the immune system, thus enhancing the antigen [41,42]. PubChem evidence suggests that Tubocurarine, despite having a top, binding score with the target, has toxic nature and hence it is neglected in our result. Homoaromoline compound is an alkaloid and presented in the roots of Thalictrum lucidum L. It has anti-microbial activity against Mycobacterium smegmatis [43]. 24-hydroxyursolic acid has chemo-preventive nature and is presented in the Diospyros kaki (Persimmon) plant. It stimulates AMP-activated protein kinase (AMPK), thus inhibiting the cyclooxygenase activity (COX-2) [44].
IL4 plays a part of the role in B-cell activation. Also, it stimulates the expression of MHC II molecules on quiescent B-cells, thus increasing the expression of immunoglobulin molecules such as IgE and IgG1 [45]. Bis (6-hydroxybenzo[b]furan-2-yl)methanone is a synthetic compound, but no complete information is available on the compound [46].
IL10 plays a role in an anti-inflammatory function where it binds with the heterotetrameric receptor. The heterotetrameric receptors are IL10RA and IL10RB, which activates the JAK1 and STAT2 [47]. It restricts the antigen-presenting cells, thus reducing the expression of MHC-II and co-stimulatory molecules owing to inhibiting the T cell activation [48]. Guggulsterone compound properties suggest anti-inflammatory function, anti-hyperlipidemic, antioxidant, immunomodulatory and craving stimulus activity. It is a bioactive compound present in plants such as Commiphora mukul and Commiphora wighti. Our study suggests that in the oral, TPSA, Blood-Brain Barrier (BBB) penetration and absorption properties are interconnected. Clinical studies supported that Guggulsterone reduces the cholesterol level in the body. It has a tremendous potential immunomodulatory effect, yet in vitro studies are insufficient to validate the results regarding Guggulsterone. But our in-silico study strongly supports that Guggulsterone is the potential compound for the IL10 gene [49]. Literature evidence is not available for (R, S)-homoaromaline hydrochloride. Up-to-date, no research works were documented on (R, S)-homoaromaline hydrochloride compound. Herkinorin is the first semisynthetic analog of Salvinorin A. Salvinorin A present in the Salvia divinorum (Lamiaceae) plant. Herkinorin has anti-nociceptive properties and has been evaluated in a study elucidating inflammatory aches in rats [50]. Mulberrofuran W is a bioactive compound present in the Morus plant. It has anti-hyperglycemic activity via inhibiting the tyrosine phosphatase and anti-carcinogenic activity via inhibiting the hypoxia-activated HIF-1 α [51].
From the functional annotation results, other diseases related to the genes have been directly or indirectly linked with the COVID-19 infection. IL10 have found to be associated with conditions such as Leishmaniasis (Visceral, Cutaneous, and Urban cutaneous), Colitis, Arthritis (Experimental, Adjuvant-induced and Collagen-induced), Inflammation, Bright Disease, Glomerulonephritis, Allergic reaction, Hypersensitivity and Rheumatoid Arthritis have been highlighted in this study. IFN-gamma production is the protective immune response against the COVID-19 infection and the resultant Natural Killer cells and CD8+ T-cell activation. At the same time, IFNG is the Th1 immune response against Leimaniasis; a case study also reported co-infection of Leimaniasis with COVID-19 in an immunocompromised patient [52]. A case study regarding the GI system proposed that Ulcerative colitis has been diagnosed in a patient with diarrhea and abdominal pain following post-treatment of COVID-19 pneumonia. Data suggests that COVID-19 might influence the GI system by ACE2 protein expression [53]). One of the most common symptoms occurring in COVID-19 patients is arthralgia. However, COVID-19 might attack the musculoskeletal system of the host in its infection or post-infection stage, leading to inflammatory arthritis, but rheumatoid manifestations are yet to be confirmed [54]. Undeniably, smokers seem to be more susceptible to developing cancer and lung diseases and are more vulnerable to COVID-19 infection and COVID-19 pneumonia because smoking can increase ACE2 protein expression, leading to the entrance of the SARS--CoV-2 virus into the host more quickly than usual [55]. The conformational changes were more in the protein backbone after protein-ligand complex interaction. Though the results of MD simulation for IL10 protein and ligand molecule complexes might not showed favorable results, the functions of the protein cannot be interpreted only based on the computational analyses.

Conclusion
Overall, our study identified many of the similar compounds showing interactions between target genes; but the Guggulsterone has reached top position among all the interactions and comparing with the previous literatures. However the MD simulation results were not promising regarding protein-ligand complexes, we cannot rely only on the computational results and might do experimental studies in future to validate the functions of IL10 and Guggulsterone in COVID severe cases. Based on this whole analysis, we identified IL10 plays a main role in antiinflammatory function and which reacts to Guggulsterone in SARS-CoV-2. The same IL10 attained the top position in miRNA network study and integrated network analysis. In specific, we recommend the compound Guggulsterone for the target IL10 against COVID-19 infection. Utterly, this is the first in silico study conducted with Guggulsterone against IL10 in COVID-19 infection. So, we hope that in future in vitro studies might confirm our study.