Pharmacoinformatics-based investigation of bioactive compounds of Rasam (South Indian recipe) against human cancer

Spice-rich recipes are referred to as “functional foods” because they include a variety of bioactive chemicals that have health-promoting properties, in addition to their nutritional value. Using pharmacoinformatics-based analysis, we explored the relevance of bioactive chemicals found in Rasam (a South Indian cuisine) against oxidative stress-induced human malignancies. The Rasam is composed of twelve main ingredients, each of which contains a variety of bioactive chemicals. Sixty-six bioactive compounds were found from these ingredients, and their structures were downloaded from Pubchem. To find the right target via graph theoretical analysis (mitogen-activated protein kinase 6 (MAPK6)) and decipher their signaling route, a network was built. Sixty-six bioactive compounds were used for in silico molecular docking study against MAPK6 and compared with known MAPK6 inhibitor drug (PD-173955). The top four compounds were chosen for further study based on their docking scores and binding energies. In silico analysis predicted ADMET and physicochemical properties of the selected compounds and were used to assess their drug-likeness. Molecular dynamics (MD) simulation modelling methodology was also used to analyse the effectiveness and safety profile of selected bioactive chemicals based on the docking score, as well as to assess the stability of the MAPK6-ligand complex. Surprisingly, the discovered docking scores against MAPK6 revealed that the selected bioactive chemicals exhibit varying binding ability ranges between − 3.5 and − 10.6 kcal mol−1. MD simulation validated the stability of four chemicals at the MAPK6 binding pockets, including Assafoetidinol A (ASA), Naringin (NAR), Rutin (RUT), and Tomatine (TOM). According to the results obtained, fifty of the sixty-six compounds showed higher binding energy (− 6.1 to − 10.6 kcal mol−1), and four of these compounds may be used as lead compounds to protect cells against oxidative stress-induced human malignancies.


Prediction of in silico pharmacokinetic and physicochemical properties.
In silico prediction of ADME and physicochemical properties of the selected bioactive chemicals plays a major role in determining its integrity and efficiency. Selected bioactive chemicals into account, properties like molecular weight, molar refractivity, solubility, bioavailability, bioavailability radar plot, egg-boiled model, brain penetration, and human gastro intestinal absorption properties of the active bio-compounds have been determined using SwissADME (http:// www. swiss adme. ch/) webserver. The SwissADME webserver is a free tool that can predict the pharmacokinetic and drug-likeness properties of the test bioactive compounds 37, 38 . Toxicity prediction. Toxicity was predicted by determining the safety profile of the intended bioactive chemicals, which must have deadly effects on people and cause organ damage. As a result, the toxicity of the chosen bioactive chemicals was assessed using pkCSM-pharmacokinetics web-based server (http:// biosig. unime lb. edu. au/ pkcsm/ predi ction) 39 .
Molecular dynamics simulation. The molecular dynamic simulation was evaluated to determine the binding stability, conformation and interaction modes between the selected bioactive compounds (ligands) and receptor (MAPK6). The selected ligand-MAPK6 complex files were subjected to molecular dynamics studies using GROMACS 2019.2 software [40][41][42] . The selected ligands topology was downloaded from PRODRG server 43 .
The system preparation of all the complexes were as described earlier 24 . For molecular dynamic simulation, first vacuum was minimized using the steepest descent algorithm for 5000 steps. The complex structure was solvated in a cubic periodic box of 0.5 nm with a simple point charge (SPC) water model. The complex system was subsequently maintained with an appropriate salt concentration of 0.15 M by adding a suitable amount of Na + and Cl − counter ions. Each complex was allowed a simulation time of 50 ns from the NPT (Isothermal-Isobaric, constant number of particles, pressure, and temperature) equilibration was subjected in NPT ensemble for final run. The trajectory analysis of root means square deviation (RMSD) and root mean square fluctuation (RMSF) was performed in the GROMACS simulation package through the online server "WebGRO for Macromolecular Simulations (https:// simlab. uams. edu/)".

Molecular mechanics Poisson-Boltzmann surface area (MMPBSA) calculation. The MMPBSA
method was used to calculate the protein-ligand binding free energy of each complex. The free energy of binding was determined using the g_mmpbsa tool developed for GROMACS 24, 44 . Consent to participate. All authors agree to participate. All  www.nature.com/scientificreports/

Results
Graph theoretical network analysis. A graph was constructed for the MAPK6/ERK3 signaling pathway network. Numerous entities of genes, proteins (nodes) and their interactions (edges) in the current research work was displayed in Fig. 1. Based on centrality criteria such as degree, proximity, eccentricity, eigen vector, and radiality, the network has 92 nodes and 153 edges. The measured values of Betweennes (3483. 24) degree (2), closeness (0.004), eccentricity (0.125), eigen vector (0.0538), radiality (10.08), and stress (20,850) have shown the threshold value of all measures as well as significant node in the network (Tables 1, 2). The protein mitogen activated protein kinase 6 (MAPK6)/extracellular signal-regulated kinase 3 (ERK3) was identified as a potential drug target based on the centrality measure and its threshold values.
Retrieval of bioactive compounds and preparation. The accessible bioactive components of the selected spices (tamarind, red pepper, black pepper, cumin seed, fenugreek, asafoetida, garlic, tomato, coriander, curry leaves, sesame oil, and mustard) were searched using IMPPAT database. From the database, a list of sixtysix bioactive compounds were obtained from the twelve spices along with one standard MAPK6 inhibitor listed in Table 3.
Binding site identification. Our evaluation of the crystal structure of MAPK6 (PDB: 7AQB) revealed the existence of 11 binding pockets, according to binding site analyses. The protein's recovered binding site residue are shown in the Fig. 2. Molecular docking investigations were also conducted using the obtained complex structure of the binding sites. Grid generation in molecular docking resulted in more reliable ligand posture scoring. As a result, we created a receptor grid for the selected MAPK6 protein based on the previously acquired binding site residues to achieve more precise scoring of our ligand poses. A receptor grid with a box dimensions of X = 38.6666, Y = 62.5914, and Z = 31.9740 in angstrom (Å) was created and utilized further for molecular docking experiments.

Molecular docking.
The optimum intermolecular interaction between the target protein and bioactive chemicals were investigated using molecular docking analysis. To analyse their binding capability, a specific number of bioactive chemicals (sixty-six) and known MAPK6 inhibitor were docked against MAPK6 using AutoDock Vina. Twelve bioactive chemicals were shown to have a lower binding energy (< − 9 kcal mol −1 ) with the target protein. The binding energy of the bioactive compounds following molecular docking was found to be scattered, ranging from − 3.50 to − 10.60 kcal mol −1 , as illustrated in Fig. 3 and Table 3. The top four compounds (Assafoetidinol A (− 9.80 kcal mol −1 ), Naringin (− 9.60 kcal mol −1 ), Rutin (− 9.80 kcal mol −1 ), Tomatine (− 10.60 kcal mol −1 )) were chosen for future research based on their binding energy with the amino acid residues in the active site of MAPK6. In this study, we used a standard MAPK6 inhibitor (CID: 447077) (− 9.30 kcal mol −1 ) as control, due to previously reported inhibitory activity against Bcr-Abl-dependent cell growth with an IC 50 of 2-35 nM in different cell lines 45 .  Fig. 8 and Table 4.

Interpretation of protein-ligand interactions.
Pharmacokinetic and physicochemical properties prediction analysis. The ADME and physicochemical properties of selected bioactive compounds from Rasam were assessed through SwissADME (http:// www. swiss adme. ch/) webserver and these are presented in Table 5. From the assessed data in Of the four compounds analyzed only ASA (CID: 12041593) and standard MAPK6 inhibitor (PD-173955) was predicted to have better chances as a possible drug-relevant candidate with anticancer potential. All the four bioactive compounds are soluble in nature, except standard drug (CID: 447077). The solubility class of all the compounds are listed in the Table 5. The synthetic accessibility score was found to be > 6 except CID: 12041593 and the standard drug (CID: 447077), which indicated that other three compounds are very difficult to synthesize. The graphical representation of drug-likeness of four selected compounds are presented in the Fig. 9. The pink area within the hexagon represents the optimal range for the compounds. The recommended range for drug-like compound was insaturation (INSATU): fraction of carbons in the sp3 hybridization not less than 0.25, insolubility (INSOLU): log S not higher than 6, hydrophobicity (LIPO): between − 0.7 and + 5.0, rotatable bonds (FLEXI):  www.nature.com/scientificreports/ no more than 9 rotatable bonds, molecular weight (SIZE): between 150 and 500 g mol −1 and polar surface area (POLAR): between 20 and 130 Å 2 ). The compound CID: 12041593 (ASA) possess drug-like properties and are displayed by the red slanted hexagon within the pink shade ( Fig. 9). Other three compounds (CID: 442428, CID: 5280805 and CID: 28523) and standard drug (CID: 447077) have high polar surface area, poor solubility and high molecular weight (size) which indicated off-shoot of the vertex (polar and size). Moreover, the pharmacokinetic properties were investigated by using egg-boiled model of selected four compounds and standard MAPK6 inhibitor is depicted in Fig. 10. The egg-boiled model is helpful to predict two key pharmacokinetic properties simultaneously, i.e., the passive gastrointestinal absorption and blood brain barrier (BBB) penetration. The egg-shaped organisation plot shows that the compound present in the yolk (i.e. yellow region) represent the highly probable BBB permeation whereas albumin (i.e. white region) represent the highly probable human intestinal absorption. From the Fig. 10, the compound ASA (CID: 12041593) and standard drug (CID: 447077) found in albumin (white region) elucidated the good absorption in gastrointestinal region. Remaining three compounds were found to be outside of the boiled-egg region. From the above observed results, it can be interpreted that the compound ASA (CID: 12041593) have sufficient potential to be drug.

Analysis of toxicity.
In silico toxicity prediction of the selected four compounds has been performed using pkCSM-pharmacokinetics web-based server. The server has identified drug-induced hERG toxicity, AMES toxicity, LD 50 , hepatotoxicity, skin sensitization, Tetrahymena pyriformis (TP) toxicity, and minnow toxicity which was listed in Table 6.  www.nature.com/scientificreports/ RMSD is an important parameter to analyse the equilibration of MD trajectories and check the stability of complex systems during the simulation process. RMSD of the protein backbone atoms were plotted against time to assess its variations in structural conformation. Initially, the 7AQB-ASA complex showed variations in backbone RMSD till 30 ns ranging from 0.15 to 0.44 nm. The stable conformation was attained in the time period between 21 and 50 ns with no considerable deviations in the values (Fig. 11). 7AQB-NAR complex showed variations in backbone RMSD till 20 ns ranging from 0.17 to 0.43 nm. The stable conformation was attained in the time period between 21 and 50 ns with no considerable deviations in the values (Fig. 11). The 7AQB-RUT complex showed variations in backbone RMSD till 30 ns ranging from 0.13 to 0.35 nm. The stable conformation was attained in the time period between 31-50 ns with no considerable deviations in the values (Fig. 8). The 7AQB-TOM complex showed variations in backbone RMSD till 35 ns ranging from 0.14 to 0.43 nm. The first stable conformation was attained in the time period between 36 and 50 ns with no considerable deviations in the   , pocket 4 (12 AA with predicted score 3.43), pocket 5 (13 AA with predicted score 3.10), pocket 6 (11 AA with predicted score 3.0), pocket 7 (10 AA with predicted score 2.14), pocket 8 (12 AA with predicted score 2.04), pocket 9 (08 AA with predicted score 1.90), pocket 10 (11 AA with predicted score 1.77) and pocket 11 (08 AA with predicted score 1.75). www.nature.com/scientificreports/ values (Fig. 11). This clearly specifies that the protein underwent small structural changes in all the complexes during simulations. RMSF is an another crucial parameter while examining the stability and flexibility of complex systems during simulation 47 . RMSF was examined to analyse the changes in the behaviour of amino acid residues of target protein on binding to a ligand 48,49 . The RMSF values for Cα atoms of the protein were calculated and plotted with respect to the residues. In all the complexes, examined, the amino acid residues showed minimal fluctuations throughout the simulation. The amino acids of MAPK6 which interacted with ASA during docking showed minimal fluctuation values during MD simulation viz. CYS28, GLY29, LYS185 and LYS229, with NAR it showed low fluctuation values during MD simulation viz GLY29 and LEU192, with RUT showed minimal fluctuation values during MD simulation viz. GLY29, ARG70, LYS229 and ASN269 and with TOM it showed moderate fluctuation values during MD simulation viz. GLY29, LYS185, SER189, TYR266 and PRO301 (Fig. 12). These results revealed that binding of the ligands actuated no major effects on the flexibility of the protein.
Further, Radius of gyration (Rg) of the complex systems were also analysed. Rg is the root mean square distance of the atoms of the protein from the axis of rotation 49 . It is one among the important parameter that represents the overall change in the protein structure compactness and its dimensions during the simulation 50 . Higher Rg values characterize the protein as less compact and flexible while low values depict the high compactness and rigidity 47 . Rg values of backbone atoms of protein were plotted against time to examine the changes in structural compactness. Binding of ASA decreased the backbone Rg values till 30 ns. In the time period between 31 and 50 ns there were no considerable fluctuations and almost constant value of ~ 1.98 nm was maintained. Till end, the Rg values were found to be in the range between 1.95 and 1.99 nm. Complete analysis revealed that, in the initial stage the trajectory had shown its peak value of ~ 2.12 nm. Later this high value was never displayed again which shows the stability of protein in the complex (Fig. 13). Binding of NAR decreased the backbone Rg values till 15 ns. In the time period between 16 and 45 ns there were no considerable fluctuations and almost constant value of ~ 2.04 nm was maintained. Till end, the Rg values were found to be in the range of 2.02-2.05 nm. Complete analysis revealed that, in the initial stage, the trajectory showed its peak value of ~ 2.09 nm. Later, this high value was never displayed again which shows the stability of protein in the complex (Fig. 13). Binding of RUT decreased the backbone Rg values till 31 ns. In the time period between 32 and 50 ns there were no considerable fluctuations and almost constant value of ~ 2.03 nm was maintained. Till end, the Rg values were found to be in the range of 2.00-2.05 nm. Complete analysis revealed that, in the initial stage, the trajectory exhibited its peak value of ~ 2.10 nm. Later, this high value was never displayed again which shows the stability of protein in the  www.nature.com/scientificreports/ complex (Fig. 13). Binding of TOM decreased the backbone Rg values till 10 ns. In the time period between 11 and 50 ns there were no considerable fluctuations and almost constant value of ~ 1.96 nm was maintained. Till end, the Rg values were found to be in the range of 1.94-2.00 nm. Complete analysis revealed that, in the initial stage, the trajectory had shown its peak value of ~ 2.10 nm. Later, this high value was never displayed again which shows the stability of protein in the complex (Fig. 13). The complete interpretation revealed that both the molecules induced no major structural changes in the protein.
Moreover, analysis of Solvent Accessible Surface Area (SASA) for all the complexes was implemented. SASA is the substantial criterion to examine the extent of exposure of receptor to the surrounding solvent molecules during simulation 47,51 . In general, binding of ligand may induce the structural changes in the receptor and hence the area in contact with the solvent also may vary 49 . SASA values of protein was plotted against time to www.nature.com/scientificreports/   www.nature.com/scientificreports/ estimate the changes in surface area. For SASA complex, the trajectory showed decrease in the values till 15 ns. Except few time intervals, minute fluctuations were observed throughout the simulation period (Fig. 14). The average SASA value was found to be ~ 138 nm 2 and were in the range of 150-130 nm 2 . For NAR complex, the trajectory showed decrease in the values till 10 ns. Except few time intervals, minute fluctuations were observed throughout the simulation period (Fig. 14). The average SASA value was found to be ~ 142 nm 2 and were in the range of 149-134 nm 2 . For RUT complex, the trajectory showed decrease in the values till 10 ns. In the time interval of 11-28 ns, minute fluctuations were observed and from 29 to 34 ns a moderate fluctuation was observed (Fig. 14). The average SASA value was found to be ~ 140 nm 2 and were in the range of 154-133 nm 2 . For TOM complex, the trajectory showed decrease in the values till 10 ns. Except few time intervals, minute fluctuations were observed throughout the simulation period (Fig. 14). The average SASA value was found to be ~ 148 nm 2 and were in the range of 154-140 nm 2 . Overall, the analyses revealed that the surface area of protein in complexes were shrunken during the simulation.
To examine the binding energy of the ligands with the target protein, the MD trajectories were analyzed to interpret the extent of hydrogen bond formation during the entire course of simulation and was depicted in Fig. 15. SASA had formed good number of H-bonds with the receptor protein with a maximum of five bonds at several time frames indicating the stronger affinity towards the target. Consistency was maintained in forming almost two hydrogen bonds for the entire simulation time which signifies the stability of the complex. For the  www.nature.com/scientificreports/ NAR complex, the consistency was maintained in forming three hydrogen bonds with maximum of six bonds at certain time periods. For RUT complex, the consistency was maintained in forming four hydrogen bonds with maximum of nine bonds at certain time periods. For the TOM complex, the consistency was maintained in forming two hydrogen bonds with maximum of nine bonds at certain time periods. This clearly signifies that the top phytochemicals have the stronger affinity with the target protein.

Discussion
The purpose of this research work was to look at the cancer-preventive impact of bioactive chemicals found in the south Indian cuisine Rasam by using graph theoretical network and pharmacoinformatics analysis. Pharmacoinformatics is a collection of in silico molecular modeling tools for screening the bioactive substances based on their binding affinities, pharmacokinetics, and pharmacodynamic features 52 . By enabling researchers to narrow down the biological and synthetic research impacts, pharmacoinformatics has sped up the discovery of bioactive substances. Several substances have their positive effects predicted using pharmacoinformatics research, which were then validated by in vitro and in vivo activities. Understanding how chemicals bind, interact, and inhibit/ stimulate a certain protein might help researchers find therapeutic options for certain disease conditions. www.nature.com/scientificreports/ Initially, a graph theoretical network was developed using centrality metrics, and it was suggested for metabolic networks that included enzymatic cascades and synergistic ligand-enzyme interactions. Biological networks, which are made up of a number of vertices (or nodes) linked in a pattern by a set of edges (or connections), are designed to mimic the structure of genuine biological systems. MAPK6/ERK3 was identified as a receptor (target) for ligand (bioactive substances) binding in ROS-induced oxidative stress that leads to malignancies, according to the network analysis study. MAPK6 pathway has been shown to be impacted not just by receptor ligand interactions, but also by various other stresses. MAPK6 is shown to be a key factor for organogenesis, cancer cell growth and invasiveness. MAPK6 overexpression has been detected in numerous human cancers, including squamous cell lung carcinoma, oral squamous cell carcinoma, gastric cancer, breast cancer and melanoma 53 . DNA microarray studies have produced inconsistent information about the regulation of MAPK6 expression in cancer. Several investigations demonstrated that expression of MAPK6 mRNA is down-regulated in brain tumors, ovarian carcinoma and cutaneous melanoma, and up-regulated in leukemias, adrenocortical carcinoma, squamous cell lung carcinoma, salivary adenoid cystic carcinoma, tongue squamous cell carcinoma and cervical cancer 54 . Furthermore, since MAPK signaling cascade regulate both mitogen-and stress-activated signals, the regulation of both pathways by ROS has drawn researchers' attention 55 . The goal of the current research work was to look Figure 10. The EGG-BOILED model for the selected bioactive and standard MAPK6 inhibitor drug. The EGG-BOILED represents for intuitive evaluation of passive gastrointestinal absorption (HIA) white part and brain penetration (BBB) yellow part as well as substrates (PGP +) and non-substrates (PGP-) of the permeability glycoprotein (PGP) are represented by blue and red color circles, respectively, of the selected bioactive compound and standard MAPK6 inhibitor in the WLOGP-versus-TPSA graph. The grey region is the physicochemical space of compounds predicted to exhibit high intestinal absorption. www.nature.com/scientificreports/ into the detoxification/neutralization of ROS by employing bioactive chemicals found in Rasam spices to protect cells against cancer. A total of sixty-six bioactive compounds were chosen from twelve spices using the IMPPAT database, as well as previously published publications on their effects against different human malignancies. All the chemicals chosen were docked against MAPK6 protein kinases, with binding energy ranging from − 3. Bioactivity of a substance is largely governed by its absorption, distribution, metabolism, and excretion (ADME) characteristics, all of which are connected to its pharmacokinetic characteristics. The bioavailability of dietary phytochemicals to target cells, as well as their absorption and metabolism in the human body are certain key aspects in promoting their bioactivity and maintaining body health 57 . The small intestine absorbs some, but not all, of the components of dietary phytochemicals into the circulatory system. Some phytochemical compounds that were absorbed by the colon and altered by the gut microbiota; the microbial metabolites that were released back into the circulation showed significant activity 58 . In order for any molecule to permeate the membrane, phytochemicals/test substance must break hydrogen bonds in the aqueous environment and partition across the membrane 59 . The polar surface area (PSA) of a chemical is connected to its hydrogen-bonding potential, whereas molecular mass and lipophilicity are associated to membrane permeability 60 . As a consequence, the ADME properties must be assessed at the earlier stages of drug design and discovery process in order to pass the standard clinical studies required to be considered as prospective therapeutic candidate 61 . In this study, all the discovered phytoconstituents were confirmed in terms of usual pharmacokinetic properties using multiple  www.nature.com/scientificreports/ bioinformatics methods. Phytochemicals are naturally derived from variety of plants that are often consumed by humans and are usually considered safe to consume. While most phytochemicals are not regulated by the Food and Drug Administration (FDA) in the United States, their potential toxicity is unknown. Phytochemicals are utilized as supplements in conjunction with illness therapy all around the world, but the users do not necessarily inform this to their physicians 62 . Substance toxicity refers to the property of any compound to be poisonous and to cause harm to an organism. Toxicity testing of a substance necessitates in vitro  www.nature.com/scientificreports/ and in vivo animal experiments, which is a time-consuming, expensive, and complicated technique. Because there are no animal trials, and improved precision, accessibility, and speed, in silico toxicity assessment has become very popular in recent times and it can offer information on any synthetic or natural molecule. In this work, in silico approaches were used to estimate the toxicity levels of four chemicals. The non-carcinogenic and non-skin irritating properties of four substances were determined using in silico testing. Three compounds, Assafoetidinol A, Rutin, and Tomatine were shown to be negative in Ames testing. Toxicity tests revealed that the four phytochemicals chosen had no negative side effects (hERG). The LD 50 (median fatal dosage) indicates the immediate or acute toxicity of substances that were determined to be the most effective in the investigation. Hence, the complexes of these compounds were subjected to molecular dynamics simulations and the results were analysed with the results of apo form of MAPK6 (7AQB). The complexes were validated by interpreting the RMSD, RMSF, Rg, SASA and the lead phytochemical complexes were found to be stable during the simulations.

Conclusion
Traditionally, home-cooked meals have been shown to help avoid chronic illnesses, improve health, and save treatment costs while also boosting quality of life. This study looked at the antioxidant properties of bioactive chemicals found in the south Indian cuisine, Rasam against oxidative stress-induced human malignancies. In the human body, ROS is a metabolic by-product of cellular respiration. Oxidative stress and overexpression of MAPK6 protein are caused by an increase in ROS levels. MAPK6 overexpression causes a cascade of events in cells, including mutations and carcinogenesis. Through a thorough pharmacoinformatics-based molecular docking investigation of bioactive substances against MAPK6, the antioxidant potential of Rasam has been proven in the current work. In silico molecular docking investigations found that the four lead phytochemicals (Assafoetidinol A, Naringin, Rutin, and Tomatine) may suppress MAPK6/ERK3 expression. In addition, MD stimulation tests and in silico pharmacokinetic prediction analyses gives the safety profile of four lead compounds as well as the stability of the protein-ligand complex, although, in order to determine the Rasam's effectiveness, further in vitro and in vivo animal research work will be necessary.