Computational pharmacology profiling of borapetoside C against melanoma

Abstract Melanoma,also known as a ‘black tumor’, begins in the melanocytes when cells (that produce pigment) grows out of control. Immunological dysregulation, which raises the risk for multiple illnesses, including melanoma, may be influenced by stress tiggered through viral infection, long term effects of ultraviolet radiation, environmental pollutants etc. Borapetoside C is one of the phytoconstituents from Tinospora crispa, and its biological source has been reported for its antistress property. Network pharmacology and KEGG pathway analysis of borapetoside C-regulated proteins were conducted to identify the hub genes involved in melanoma development. Further, a molecular docking was performed between borapetoside C and targets involved in melanoma. Further, the top 3 complexes were selected based on the binding energy to conduct molecular dynamics simulations to evaluate the stability of ligand-protein complex followed by principal component analysis and dynamic cross-correlation matrix. In addition, borapetoside C was also screened for its pharmacokinetics and toxicity profile. Network Pharmacology studies and KEGG pathway analysis revealed 8 targets involved in melanoma. Molecular docking between borapetoside C and targets involved in melanoma identified 3 complexes with minimum binding i.e. borapetoside C- MAP2K1, MMP9, and EGFR. Further, molecular dynamics simulations showed a stable complex of borapetoside C with MMP9 and EGFR. The present study suggested that borapetoside C may target MMP9 and EGFR to possess an anti-melanoma property. This finding can be useful in developing a novel therapeutic agent against melanoma from a natural source. Communicated by Ramaswamy H. Sarma


Introduction
Melanocytes produce the ultraviolet (UV)-absorbing pigment melanin and are located at the basal level of the epidermis.
Melanocyte-stimulating hormone is secreted by keratinocytes in response to UV light exposure, and this hormone subsequently binds with the melanocortin 1 receptor to produce melanin (Williams et al., 2011).Melanocytes produce both eumelanin and pheomelanin.A decreased risk of developing skin cancer is associated with a higher concentration of the UV-protective pigment eumelanin, which is present in darker skin tones.The formation of pheomelanin not only results in carcinogens (Seiberg, 2001;Morgan et al., 2013;Mitra et al., 2012) but also provides less protection from UV light.There is evidence that pheomelanin increases UV-induced reactive oxidative species, which in turn increases deoxyribonucleic acid (DNA) damage (Seiberg, 2001;Morgan et al., 2013;Gajula & Gaddameedhi, 2015;Robles-Espinoza et al., 2016).Cancer data from the Centers for Disease Control and Prevention show that 22.1% of every 100,000 Americans develop melanoma (malignant tumor originating from melanocytes).Despite only making up 4% of skin cancer incidences, it is responsible for 75% of skin cancer-related fatalities.Therefore, there is a need for new chemotherapeutic agents against melanoma with minimal side effects.
Tinospora species are frequently employed in multiple forms of traditional medicine for a variety of reasons (Gray-Schopfer et al., 2007).Previously, the use of T. cordifolia lotion has been recommended against Sarcoptes scabiei infection to combat scabies (Castillo et al., 2013).In addition, T. cordifolia exerts antiosteoporosis, anti-diabetic, hypolipidemic, anticancer, anti-HIV, antitoxic, immunomodulating, wound healing, and antioxidant effects (Sharma et al., 2019).The adaptogenic capacity of T. cordifolia was reported to boost physical performance while at the same time inhibiting excessive activation of the sympathetic nervous system (Salve et al., 2015).
It has been reported that stress hormones like norepinephrine promote the function of cytokines like interleukin 6 and 8, which are proangiogenic and assist tumor growth.It is believed that stress, in conjunction with genetic and environmental variables, may contribute to developing melanoma and its progression (Sinnya & De'Ambrosis, 2013).In addition, previously, T. crispa has been reported to inhibit MMP-13 and also check the migration of squamous cell carcinoma cell lines (Phienwej et al., 2015) followed by antiproliferative activities on different human cancer cell lines (Zulkhairi et al., 2008).Further, borapetoside C (a furanoid diterpene glycoside) is one of the active phytoconstituent from T. crispa (Ruan et al., 2012;Hossen et al., 2016).Based on these study backgrounds, we attempted to screen borapetoside C (Figure 1) as a potential anti-melanoma compound.

Enrichment analysis
Borapetoside C-modulated melanoma proteins were enriched in a search tool for the retrieval of interacting genes/proteins (STRING) database (Snel et al., 2000; https://string-db.org/) in a full network at 0.4 confidence and 5% FDR stringency to trace gene ontology terms and probably modulated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways at 5% false discovery rate stringency for Homo sapiens.The protein-protein interactions were evaluated based on known interactions (curated databases and experimentally determined), predicted interactions (gene neighbourhood, gene fusions, and gene co-occurrence), and miscellaneous (text mining, co-expression, and protein homology) by assuming whole genome statistical background.

Disease-protein network construction and analysis
Before constructing the network, all the KEGG disease records were evaluated for 'cancer' by using 'melanoma ', 'sarcoma', 'cancer', 'glioma', 'leukemia', and 'carcinoma' keywords.The disease-proteins-borapetoside C network was constructed using Cytoscape (Shannon et al., 2003; https:// cytoscape.org).The protein-protein network was evaluated by treating it directed.Also, the disease (KEGG pathways) -protein-borapetoside C network was evaluated by treating it undirected by setting the node size as 'low values to small sizes' and node color from 'low values to bright colors' (Dwivedi et al., 2021).

Ligand molecule preparation
Borapetoside C (PubChem CID: 101697033) in .sdfformat was retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/)database, which was further converted into .pdbformat using Discovery Studio Visualizer ver.2021, energy was minimized and later used as a ligand.

Prediction of the active site in receptors
The active site of macromolecules was identified using the Indian Institute of Technology, Delhi's supercomputing facility for bioinformatics and computational biology server (http://www.scfbio-iitd.res.in/dock/ActiveSite.jsp) to search all of the possible binding cavities.The cavity possessing the highest volume was selected to dock the ligand as it would be the most suitable site for molecular docking.The cavity points were noted for performing docking with Auto Dock Tools 1.5.6 (Morris et al., 2009; https://autodock.scripps.edu/).

Molecular docking of borapetoside C with proteins involved in melanoma
Borapetoside C was docked against BRAF (x: À 7.  (x: 34.362, y: 44.387, and z: 24.453) using Auto Dock Tools at 60 dimensions (x, y, and z) and 0.5 Å space.After docking, the ligand posewith the minimum binding energy was selcted to visualize the ligand-protein interaction in Biovia Discovery Studio 2021 (Bhattacharya et al., 2022).

Molecular dynamics simulation
To examine docked complexes' structural and intermolecular interaction stabilities, an all-atom explicit molecular dynamics (MD) simulation for 100 ns was carried out.The GROMACS (https://www.gromacs.org)2021.3 software package with Amber ff99SB-ildn force field was employed during the simulation as demonstrated previously (Van Der Spoel et al., 2005;Dwivedi et al., 2022).The topological parameters of the ligand and the whole complex were generated using the xleap module of AmberTools (https://ambermd.org/AmberTools.php), and the partial charges of the small molecules were obtained by doing quantum calculations using an antechamber with a 'bcc' charge model.The prepared system was solvated in a rectangular box with 10.0 Å boundary conditions from the protein's borders in all directions using the 3-site water (TIP3P) model.The charges on the prepared system were neutralized by introducing the required number of counter ions.To find the near-global state least energy conformations, the steepest descent followed by the conjugate gradient energy minimization method was applied.Canonical (NVT) and isobaric (NPT) ensembles were used to equilibrate the system for 1 ns.A modified Berendsen thermostat method was used in NVT equilibration to maintain the volume and temperature constant (300 K).Parrinello-Rahman barostat was used to maintain the pressure at 1 bar constant during NPT equilibration.In addition, the Particle Mesh Ewald approximation was used with a cut-off value of 1 nm to calculate the long-range electrostatic interactions, van der Waals, and coulomb interactions.Similarly, bond length was constrained using the LINear Constraint Solver algorithm.Finally, each complex was simulated for a 100 ns production run with coordinates at every 2 fs.The trajectories produced were examined using the built-in gromacs utilities.The MD plots were constructed using xmgrace.

Evaluation of binding affinity using molecular mechanics poisson-boltzmann surface area
In MD simulations and thermodynamic calculations, the relative binding energy of a ligand-protein complex was employed to evaluate the binding affinities.In the present study, the relative binding energy and its contribution to individual residues were calculated via the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method using the 'g mmpbsa' tool (Kumari et al., 2014).The parameters from past studies (Bhandare et al., 2019;Bhandare & Ramaswamy, 2018) were considered while calculating the binding energy.The binding energy was determined throughout the stable trajectory observed between 50 and 100 ns using 50 representative snapshots.

Principal component analysis
Using MD trajectories, principal component analysis investigates the molecular motion i.e. translational and rotational mobility of the molecule using the 'least square fit' to the reference structure (Amadei et al., 1993;Amadei et al., 1996;Moharana et al., 2022).The direction of the molecule's motion is reflected in a group of eigenvectors that are produced by diagonalizing a covariance matrix which was created by a linear transform of cartesian coordinate space.The energy contribution of each eigenvector to the motion shown by the eigenvalue is associated with the eigenvector.
The 'time-dependent movements' that the components carry out in a certain vibrational mode are demonstrated by projecting the trajectory onto a particular eigenvector.The atomic vibrational components' contribution towards this form of coordinated motion is demonstrated by the projection's time average.Using the built-in gromacs utilities 'g_covar', the eigenvectors and eigenvalues of the trajectory were produced by computing and diagonalizing the covariance matrix.Additionally, the eigenvectors were examined and illustrated using the 'g_anaeig' program.

The dynamic cross-correlation matrix
To determine if motion between atom pairs is correlated (positive or negative), the dynamic cross-correlation matrix was used to measure the magnitude of all pairwise cross-correlation coefficients (Khanal et al., 2021;Khanal et al., 2022).
In this section, using MD trajectory we examined each dynamic cross-correlation matrix component.C ij ¼ 1 denotes that the fluctuations of i and j have the same period and phase (positive correlation), C ij ¼ 0 denotes that there is no correlation, and C ij ¼ 1 denotes that the fluctuations have a negative correlation.

Pharmacokinetic, and toxicity profile and targets of borapetoside C
Borapetoside C was predicted to be soluble in water with low gastrointestinal absorption (according to boiled egg theory) and 2.62 LogP.In addition, it was predicted as a p-glycoprotein substrate, a noninhibitor of CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4 and non-blood-brain barrier permeable.Further, borapetoside C side effects were predicted as nephrotoxic with a pharmacological activity (Pa) value of 0.377 (Pa value 0 to 1 represents low to high toxicity).

Gene ontology analysis
The protein-protein interaction (PPI) of 65 proteins (nodes) had 306 different interactions (edges) with a 9.42 average node degree, 0.585 average local clustering coefficient, and 97 expected number of edges at enrichment p < 1.0e-16 (Figure 3a).In addition, CASP3, EGFR, MMP9, MAPK3, and MMP2 were identified as major hub genes (Figure 3b).The detailed scores of each protein for average short path length, clustering coefficient, closeness centrality, eccentricity, stress, degree, neighbourhood connectivity, number of directed edges, topological coefficient, Edge count, in degree, and outdegree are summarized in Figure 4.

Active site of macromolecule and molecular docking
Borapetoside C was docked against macromolecule at the following identified active sites i.e.BRAF (cavity point¼ À 7. Borapetoside C was predicted to have the minimum binding energy of À 12.15 kcal/mol with MAP2K1 with an inhibition constant of 1.25 nM.Borapetoside C was predicted to have 2 H-bonds in the docked complex with amino acids Leu54 and His119 of MAP2K1.This was followed by a complex between borapetoside C and MMP9 with a binding energy of À 12.03 kcal/mol by interacting with Asp185 and Gly186 and a complex between borapetoside C and EGFR having a binding energy of À 11.77 kcal/mol interacted with Met98.The least affinity was predicted between borapetoside C and MMP2 with a binding energy of À 8.86 kcal/mol and interacted with Asp26 (Table 1).Further, the borapetoside C interaction with BRAF, CASP3, EGFR, MAP2K1, MAPK1, MAPK3, MMP2, and MMP9 is presented in Figure 9.

Stability of EGFR-borapetoside C complex
The EGFR-borapetoside C complex trajectory demonstrated stable dynamics throughout the 100 ns MD simulation (Figure 10).The average root means square deviation (RMSD) value for the backbone and complex was observed to be �2.5 Å and �3.2 Å, respectively.The complex reached stability after an equilibration period of � 10 ns i.e. from 10 to 100 ns, and the complex RMSD achieved the optimum geometric conformation.The N and C terminal amino acid residues (Gly1 to Glu16 and Arg282 to Met292) forming a loop showed maximum fluctuation (� 8 Å and � 6.2 Å) respectively, and the residues involved in ligand interaction (Thr766,Cys751,Thr830,Asp831,Leu820,Cys773,Val702,Leu768,Leu694,and Val704) showed minimal residual fluctuations up to 2 Å.The radius of gyration (Rg) explains the structural folding and compactness of the molecule.The Rg value revealed stable folding during the simulation by exhibiting a steady decrease in Rg value from 20.5 to 20 Å and showed stable complex formation.Similarly, the solventaccessible surface area (SASA) was analyzed to distinguish the protein compactness behavior.The initial and final average surface area occupied by EGFR-borapetoside C docked complex was 160 nm 2 and 155 nm 2 respectively.The complex formed a compact globular shape as revealed by a steady decrease in both Rg and SASA values.This complex formed 7 stable H-bonds, of which 4 were consistent throughout the simulations.Further, the relative binding affinity between borapetoside C and EGFR was investigated using the MMPBSA approach by the g_mmpbsa tool.The estimated relative binding energy was À 66.401 kJ/mol.Further, the estimated van der Waal, electrostatic, polar solvation, and SASA energy was À 135.646,À 149.381, 234.865, and À 66.401 kJ/mol.The residues contributing most to the binding energy were identified by calculating residue decomposition energy.The residues Ser696, Gly697, Phe699, Gly700, Thr701, Val702, and Ile735 favored the stable complex formation.However, Lys721, Glu738, and Asp831 residues did not favor the interactions.Among these residues, Val702 showed significant contributions to the binding energy as it had the least contribution energy (À 6.447 kJ/mol).However, other residues also possessed the contribution energy between À 2.581 to À 3.65 kJ/mol and participated in stable complex formation.

Stability of MAP2K1-borapetoside C complex
In the MAP2K1-borapetoside C complex, the RMSD of backbone atoms was found to be stable (� 9.0 Å) throughout the 100 ns simulation (Figure 11).However, the RMSD of the complex was observed to be stable till 30 ns, and a large fluctuation was observed throughout 100 ns (� 16 Å).The N terminal residues showed maximum fluctuation, and the residues involved in ligand interaction (Gly61, Leu63, Ala132, Lys57, His119, and Leu54) showed minimal residual fluctuations up to 2.5 Å.However, the ligand was bound adjacent to the N-terminal residues (the N-terminal residues forming the longest loop possessed higher fluctuation) which allowed the ligand to get escape from the binding pocket.The Rg and SASA explain the structural folding and compactness of the molecule.The Rg value of protein revealed stable folding during the simulation by exhibiting a steady decrease in Rg value from 25 to 23 Å and revealed stable complex formation.
Similarly, the solvent-accessible surface area (SASA) was analyzed, and the initial and final average surface area occupied by MAP2K1-borapetoside C docked complex was 220 nm 2 and 210 nm 2, respectively.From 0 to 10 ns, a sharp decrease in the Rg and SASA indicates an unstable binding of borapetoside C to MAP2K1.This complex formed 8 H-bonds, of which none of the interactions was consistent throughout the simulations.Further, the relative binding affinity between borapetoside C and MAP2K1 was investigated using the MMPBSA approach by the g_mmpbsa tool.The estimated relative binding energy was 1926.877.401kJ/mol.Further, the estimated van der Waal, electrostatic, polar solvation, and SASA energy were À 0.002, 0.122, 1926.671, and 0.086 kJ/mol.The positive relative binding energy of the complex is due to the unstable binding mode of borapetoside C with MAP2K1.The residues contributing most to the binding energy were identified by calculating residue decomposition energy.All the residues exhibited significant positive energy contribution for opposing the borapetoside C binding to the MAP2K1.None of the residues participated in ligand binding and stable complex formation.

Stability of MMP9-borapetoside C complex
In the MMP9-borapetoside C complex, a steady increase in the RMSD of backbone and complex was observed till 50 ns (from � 3 to 7.5 Å) and was found to be stable throughout 100 ns with slight fluctuation at 70 ns (Figure 12).The steady increase in the RMSD could be due to the longest loop region (Arg95 to His119) adjacent to the binding pocket showing relatively higher fluctuation compared to ligand binding residues, which allowed ligands to get buried into the binding pocket for stable complex formation.This was confirmed by the steady decrease in the Rg and SASA values.Interestingly, both backbone and complex atoms maintain a similar trend throughout the 100 ns MD simulation, indicating a formation of stable complex formation.The residues involved in ligand interaction (Tyr48, Tyr52, Leu39, Asp1185, Gly186, Thr37, and Asn38) showed minimal residual fluctuations up to 3 Å.The residues forming loop 'Ser300 to Arg312' showed large fluctuation (� 4.0 to 7.5 Å).The Rg and SASA were also found to be stable throughout the 100 ns simulation.The overall Rg value range was between 26 to 25 Å.The initial Rg value was found to be � 26 Å and gradually decreased to 25 Å during simulation.Similarly, the SASA range was within � 230 to 210 nm 2 with an average SASA of � 220 nm 2 .The initial SASA value was found to be � 230 nm 2 and gradually decreased to 210 nm 2 during simulation.A stable complex was formed during the MD simulation and established a compact globular shape as revealed by a decrease in the Rg and SASA value.This complex formed 7 stable H-bonds, of which 5 were consistent throughout the simulations.Further, the relative binding affinity between borapetoside C and MMP9 was investigated using the MMPBSA approach by the g_mmpbsa tool.The estimated relative binding energy was À 124.553 kJ/mol.Further, the estimated van der Waal, electrostatic, polar solvation, and SASA energy was À 257.702, À 109.904, 267.117, and À 24.064 kJ/mol.The residues contributing most to the binding energy were identified by calculating residue decomposition energy.The residues Leu44, Glu46, Tyr48, Ala93, Met94, Arg95, and Met422 favored the stable complex formation.However, Thr96, Lys92, and Asp103 residues did not favor the interactions.Among these residues, Tyr48 showed significant contributions to the binding energy as it had the least contribution energy (À 6.994 kJ/mol).However, other residues also possessed the contribution energy between À 2.23 to À 5.39 kJ/mol and participated in stable complex formation.

Principal component analysis of complexes
The maximum collective motion of the complex was investigated by the first 2 principal components (PC1 and PC2) and captured by the first 50 eigenvectors.It was observed that the borapetoside C-MAP2K1 complex showed higher conformational flexibility (PC1: À 20 nm to 10 nm and PC2: À 10 nm to 15 nm) and a larger diversity of conformations (eigenvalue: 35 nm 2 ) during the simulations (indicated in red in Figure 13B1 and B2).Similarly, borapetoside C-EGFR complexes showed lower conformational flexibility and (PC1: À 8 nm to 3 nm and PC2: À 6 nm to 4 nm) smaller diversity of conformations (eigenvalue: 4 nm 2 ) (indicated in black in Figure 13A1 and A2).Likewise, borapetoside C-MMP9 showed moderate conformational flexibility (PC1: À 13 nm to 15 nm and PC2: À 7 nm to 7 nm) and moderate diversity of conformations (eigenvalue: 37 nm 2 ) during the simulations (indicated in green in Figure 13C1 and C2).

Dynamic cross-correlation matrix
Structural information on the coordinated motion of ligandbinding domains was gained from the observation of dynamic cross-correlation in complexes and the concerted residual motion in all the simulated complexes is presented in Figure 14.In complex borapetoside C-EGFR, the binding site residues show anti-correlation with the N-terminal domain of the EGFR.The amplitude of anticorrelation is maximum compared to borapetoside C with MAP2K1 and MMP9 complexes.Borapetoside C -EGFR complex showed a strong correlation between the residues 700 to 830.The cooperative motion expressed by the binding pocket residues ranging from 700 to 830 with the N-and C-terminal region revealed the significance of the active site residues in stabilizing the borapetoside C -EGFR complex whereas, in the borapetoside C-MAP2K1 complex, binding pocket correlation was lost.However, in borapetoside C-MAP2K1 complex, binding pocket residues ranging from 130 to 270 revealed a strong correlation.Thus, we propose binding of borapetoside C with EGFR and MMP9 would favor the conformational transition and promote the stable complex formation with enhanced non-bonded interactions compared to the borapetoside C-MAP2K1 complex.

Discussion
Melanoma is developed in the melanocytes when these pigment-producing cells start proliferating uncontrollably (Rotte & Bhandaru, 2016).It is one of the most dreaded types of skin cancer, with the highest number of cases recorded in Australia in 2020 (https://www.wcrf.org/cancer-trends/skincancer-statistics/).The current pharmacotherapy for melanoma involves the use of surgery, radiation, and chemotherapy, including drugs like cisplatin, temozolomide, and dacarbazine (Velho, 2012).Tinospora species (mainly T. crispa) is the major source of borapetoside C, and also reports suggest polysaccharide fraction T. cordifolia to be effective in lowering the metastatic potential of B16F-10 melanoma cells (Sharma et al., 2019;Leyon & Kuttan, 2004).In addition, human cancer cell lines HepG2, HL-60 (human promyelocytic leukemia cells), and Hep3B cell lines were suppressed by the methanol extract of T. crispa stem and displayed a dosedependent activity (Ibahim et al., 2011;Ahmad et al., 2016).Since, borapetosideC is one of the important active bioactive of T. crispa, the present study aimed to identify the effect of borapetoside C against melanoma through multiple computationalapporaches like gene ontology enrichment analysis, molecular docking, molecular dynamic simulations, MMPBSA analysis etc.
In the early stage of drug discovery, system biology tools play an evident role in proposing the possible mechanism of test agents against disease (Dwivedi et al., 2022;Patil et al., 2022a).Similarly, in the present study, we retrieved targets modulated by borapetoside C which were further crossmatched with the melanoma-linked targets to perform enrichment analysis.Further, KEGG pathway analysis was performed to identify hub genes and evident pathways modulated by borapetoside C against melanoma.In addition, molecular docking and simulation studies were performed  on the identified hub genes to retrieve ligand-protein interaction and stability.Enrichment analysis identified borapetoside C to modulate 65 genes (61.90% of total borapetoside C-modulated) in melanoma pathogenesis.Further, the gene ontology analysis predicted CASP3, EGFR, MMP9, and MAPK3 as major hub genes, and membrane raft (GO: 0045121) was identified as the lead cellular component; membrane rafts are small (10-200 nm), heterogeneous, highly dynamic, sterol-and sphingolipid-enriched membrane domains that compartmentalize cellular processes (Pike, 2006).Similarly, catalytic activity, acting on a protein (GO:0140096) was predicted as the lead molecular function, and protein metabolic process (GO:0019538) as the lead biological process indicates the effect of borapetoside C via catalytic activity on proteins at the membrane rafts of the cell.Further, KEGG pathway analysis predicted pathways in cancer (hsa05200) to be the major KEGG pathway.Further, the melanoma pathway was identified to be modulated by eight genes (MAPK1, CDK4, MAPK3, IGF1R, EGFR, BRAF, MAP2K1, and FGFR1).
Molecular docking revealed borapetoside C-dual specificity mitogen-activated protein kinase 1 complex to possess the highest binding affinity.High somatic mutation rates in melanoma are known to exhibit the distinctive features of UVinduced DNA repair (Sample & He, 2018).Recently, it has been discovered that Spitz neoplasms are significantly impacted by structural rearrangements in MAPK genes.According to reports, non-canonical BRAF mutations in melanomas were associated with gain-of-function mutations in MAP2K1 and MAP2K2 (MEK1 and MEK2, respectively), leading to constitutive ERK phosphorylation and increased resistance to MEK inhibitors.Recurrent somatic MAP2K1 and MAP2K2 mutations were found in a larger sample of melanoma patients and occurred at an average frequency of 8% (Sunshine et al., 2020).
The top 3 hub genes were identified through molecular docking were further subjected to molecular dynamics simulations of 100 ns.In the present study, molecular docking revealed the interaction of borapetoside C with the active site residue of EGFR, MAP2K1, and MMP9.During 100 ns MD simulation, borapetoside C was identified to interact stably with the EGFR and MMP9 but not with MAP2K1.In the borapetoside C-EGFR complex, initially, borapetoside C was predicted to interact with Thr766, Cys751, Thr830, Asp831, Leu820, Cys773, Val702, Leu768, Leu694, and Val704 residues.However, after MD simulation, the residual contribution energy revealed borapetoside C to interact stably with   Ser696, Gly697, Phe699, Gly700, Thr701, Val702, and Ile735.Among these, Val702 was found to be a common interactive residue, and interestingly other residues were also found to be within the active site region.Secondly, in the borapetoside C-MAP2K1 complex, borapetoside C initially formed interaction with Gly61, Leu63, Ala132, Lys57, His119, and Leu54, and during simulation, the trajectory analysis revealed an unstable binding mode of borapetoside C to the active site .Hence, it may not be a suitable lead candidate against MAP2K1.Similarly, in the borapetoside C-MMP9 complex, borapetoside C initially formed interaction with Tyr48, Tyr52, Leu39, Asp1185, Gly186, Thr37, and Asn38, and during MD simulation, borapetoside C showed stable contacts with Leu44, Glu46, Tyr48, Ala93, Met94, Arg95, and Met422.Among these, Tyr48 was found to be a common interactive residue and scored the least contribution energy of À 6.994 kJ/mol.The MD analysis concludes that borapetoside C is a potent lead hit toward the active residues of EGFR and MMP9.Additionally, the principal component analysis revealed that the maximum collective motion of the borapetoside C with EGFR and MMP9 complexes and showed lower conformational flexibility with a smaller diversity of conformations compared to MAP2K1.Thus, we propose that borapetoside C with EGFR and MMP9 could be more effective.
Further, an examination of the EGFR gene expression in melanoma using fluorescence in situ hybridization (FISH) revealed that EGFR gene amplification is associated with a worse prognosis and that EGFR protein expression is more frequently seen in patients who have a positive sentinel lymph node.Thicker tumors were linked to the presence of EGFR polysomy (Sunshine et al., 2020;Vidal et al., 2020).Additionally, studies have shown that MEK inhibitors like trametinib and cobimetinib can help with chemotherapy resistance.Thus, to improve the efficacy of current treatment regimens and combat the formation of acquired resistance, new molecular targets and treatment approaches are required.Such a target might be EGFR, which appears to be a key player in the process.Thus, using its inhibitors in anti-melanoma therapy may be advantageous for some melanoma patients (Lim et al., 2017;Chen & Davies, 2014;Patil et al., 2022b).A type of peptidase known as matrix metalloproteinases (MMPs) can alter the extracellular matrix by encouraging tumor-invasive activities.MMP9 is the peptidase that contributes the most to the growth of cancer, including melanoma.Abnormal overexpression of MMP9 can be brought on by dysfunctions in the PI3K/Akt and MAPK signaling pathways.Non-coding RNAs and other proteins, such as osteopontin and tissue inhibitors of metalloproteinases, that activate or inhibit MMP9 are also linked to the abnormal synthesis of the MMP9 (Niland et al., 2022).

Conclusions
The current study utilized a series of computational pharmacology tools to evaluate borapetoside C as an anti-melanoma agent.Network pharmacology and KEGG pathway analysis identified 8 hub genes to be modulated with borapetoside C treatment.Further molecular docking of these genes with borapetoside C suggested MAP2K1to possess the highest binding affinity followed by MMP9 and EGFR.The stability of all the 3 complexes was further evaluated using molecular dynamics simulations.Further, borapetoside C-MMP9 and borapetoside C-EGFR complexes were found to be stable, suggesting MMP9 and EGFR as probable targets for borapetoside C to combat against melanoma.In addition, the present study is solely based on the computational approach and hence, these findings need to be further validated with suitable wet-lab experiments.
Edit, Review, Revise, Finalize.PSRD, NRC, RKC, SD, AC: Review, Draft.All authors have read and approved the presented version of the research.

Figure 3 .
Figure 3. a: Protein-protein interaction (visualized in a: STRING and b: Cytoscape after analysis) of the borapetoside C-triggered protein.Node color; colored nodes: query proteins and first shell of interactors, white nodes: second shell of interactors, node content; empty nodes: proteins of unknown 3D structure, filled nodes: some 3D structure is known or predicted, known interactions; from curated databases, experimentally determined, predicted interactions; gene neighbourhood, gene fusions, gene co-occurrence & others; text mining, co-expression, protein homology.

Figure 5 .
Figure 5. Gene ontology analysis of borapetoside C regulated genes for cellular components, molecular function, and biological process.The gene ontology terms were evaluated at 5% false discovery rate stringency for Homo sapiens.

Figure 6 .
Figure 6.Interaction representation of the disease-gene association.The interaction was evaluated at a 5% false discovery rate stringency for Homo sapiens against a whole genome statistical background run.

Figure 9 .
Figure 9. Interaction of borapetoside C with (a) BRAF, (b) CASP3, (c) EGFR, (d) MAP2K1, (e) MAPK1, (f) MAPK3, (g) MMP2, and (h) MMP9.The green bond presents the H-bond interaction between the ligand and the respective protein whereas the other bond represents the hydrophobic interactions.The pink and green shade around the ligand represents the H-bond donor and acceptor respectively.A 2D sketch of borapetoside C interaction with the above-mentioned targets is provided in the supplementary figure

Figure 11 .
Figure 11.MD simulation of MAP2K1 -borapetoside C complex for 100 ns.(a) the RMSD plot of backbone and complex, (b) RMSF plot of C-a, (c) protein Rg plot, (d) protein SASA, (e) number of H-bond interactions formed between MAP2K1 and Borapetoside C, and (f) contribution energy plot highlighting the importance of the binding pocket residues in stable complex formation.

Figure 12 .
Figure 12.MD simulation of MMP9 -borapetoside C complex for 100 ns.(a) the RMSD plot of backbone and complex, (b) RMSF plot of C-a, (c) protein Rg plot, (d) protein SASA, (e) number of H-bond interactions formed between MMP9 and Borapetoside C, and (f) contribution energy plot highlighting the importance of the binding pocket residues in stable complex formation.

Figure 13 .
Figure 13.Principal component analysis of complexes of borapetoside C with (A) EGFR, (B) MPA2K1, and (C) MMP9.A1, B1, and C1 represent projections of MD trajectories on two eigenvectors corresponding to the first 2 principal components, respectively and A2, B2, and C2 represent the first 50 eigenvectors plotted versus eigenvalue for A), B), and C) respectively.

Table 1 .
Binding affinity, inhibitory constant, and interaction of borapetoside C with targets.