Abstract The present study employed a comprehensive approach of network pharmacology, molecular dynamic simulation and in-vitro assays to investigate the underlying mechanism of the anti-osteoarthritic potential of Vanda tessellata extract (VTE). Thirteen active compounds of VTE were retrieved from the literature and the IMPPAT database. All of these passed the drug likeness and oral bioavailability parameters. A total of 535 VTE targets and 2577 osteoarthritis related targets were obtained. The compound-target-disease network analysis revealed vanillin, daucosterol, gigantol and syringaldehyde as the core key components. Protein-protein interaction analysis revealed BCL2, FGF2, ICAM 1, MAPK1, MMP1, MMP2, MMP9, COX2, STAT3 and ESR1 as the hub genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed AGE-RAGE signalling pathway, HIF-1 signalling pathway and ESR signalling pathway as the major signalling pathway of VTE involved in treating osteoarthritis. Molecular docking analysis showed daucosterol and gigantol to have good binding affinity with BCL2, ESR1 and MMP9, and the results were further confirmed through molecular dynamics simulation analysis. The mechanism predicted by network pharmacology was validated in vitro on IL-1β-induced SW982 synovial cells. VTE did not show any cytotoxicity and inhibited the migration of SW982 cells. VTE inhibited the expression level of IL-6, IL-8, TNF-α, PGE-2, MMP-2 and MMP-9 in a dose-dependent manner. VTE inhibited nuclear translocation of NF- κβ and suppressed phosphorylation of p38, extracellular signal-regulated kinase (ERK), and c-Jun NH2-terminal kinase (JNK) of the mitogen-activated protein kinase (MAPK) signalling pathway. The results showed that VTE exerted an anti-osteoarthritic effect by a multi-target, multi-component and multi-signalling pathway approach. Keywords: Vanda tessellata extract, Molecular docking, MD simulation, Network pharmacology, Osteoarthritis Graphical abstract [35]Image 1 [36]Open in a new tab 1. Introduction Osteoarthritis (OA) is one of the most prevalent degenerative age-related joint diseases, affecting nearly 300 million individuals worldwide [[37]1]. It affects almost every joint of the body, including the knees, feet, hands, and hips and is characterised by joint discomfort, swelling, and stiffness, which can hamper a person's ability to carry out their routine work [[38]2,[39]3]. The clinical characteristics of OA include the breakdown of joint cartilage, thickening of subchondral bone, osteophytes development, inflammation of the synovial lining, internal joint inflammation, and changes in the structure of the joint capsule, ligaments, and surrounding muscles [[40][4], [41][5], [42][6]]. OA risk factors include lifestyle choices, dietary habits, aging, genetic predisposition, obesity, and mechanical injuries [[43]7,[44]8]. Currently, there is no cure for OA, and the drugs available in the market primarily focuses on symptoms management rather than addressing the underlying cause. Drug therapies like non-steroidal anti-inflammatory drugs (NSAID) or selective COX-2 inhibitors are frequently employed to alleviate joint pain. However, their prolonged use often leads to gastrointestinal complications [[45]9,[46]10]. Intra-articular corticosteroids offer a strong anti-inflammatory effect, which helps to reduce joint pain and swelling but exhibits extensive side effects [[47]11]. TNF-α, IL-1 and IL-6 antagonists show promising therapeutic effects, yet they make the body more prone to inflammation and induces tumour-related diseases [[48]12,[49]13]. Hence, there is an urgent need to develop a drug with significant curative properties devoid of any side effects. For centuries, plants have been utilised for their healing properties, offering a vast range of beneficial compounds that have traditionally been used without causing any side effects [[50]14]. Plant-derived products have been reported to alleviate osteoarthritis symptoms [[51]15]. The medicinal values of many orchid species have been recognized since ancient times [[52]16]. Many plants are also reported for their efficiency in treating arthritis in ayurvedic system of medicine [[53]17]. Bioactive compounds derived from orchids possess high medicinal potential and can act as a source of drugs [[54]18]. Vanda tessellata (Roxb.) Hook. ex G.Don (Syn.Vanda rouxburghii) is an epiphytic subshrub of the family Orchidaceae and is distributed across different parts of India, Nepal, Bangladesh, Burma, Vietnam, China and Sri Lanka [[55]19]. In the Unani system, the plant juice is used for treating toothaches, bronchitis, boils and as a tonic for the brain and liver [[56]20]. Several pharmacological properties of Vanda tessellata such as anti-inflammatory, anti-convulsants, anti-diarrhoea, antioxidant, anti-microbial, wound healing, neuroprotective and hepatoprotective activity, have also been previously reported [[57][21], [58][22], [59][23], [60][24], [61][25], [62][26], [63][27], [64][28]]. Different parts of the plant are used to treat pain, inflammation, arthritis, sciatica, liver disease, bronchitis, hiccough, and fever [[65]29]. The root of Vanda tessellata is reported to exhibit activity against bacterial infection, tuberculosis, abdominal diseases, tremors, and otitis [[66]30]. Root & leaf extracts of Vanda tessellata possess significant analgesic and anti-inflammatory activity [[67]28]. Phytochemical screening studies have indicated the presence of different chemical classes such as polyphenols, alkaloids, terpenoids, flavonoids, tannins, fatty acids, saponins, glycosides and sitosterol in the plant [[68]31]. The conventional drug discovery approach has some limitations due to the complex nature of diseases that involve interaction among various molecules and pathways. Hence, targeting a single molecule might not be enough to treat a disease effectively. Network pharmacology is an emerging area of study that uses complex interaction networks to explore complex biological system and their interaction with the drug [[69]32,[70]33]. Although a previous study on an arthritic rat model has suggested the potency of Vanda tessellata against arthritis [[71]34], the underlying molecular mechanisms is still not clear. Additionally, there is a lack of a comprehensive understanding of the pharmacodynamic properties of the active constituents of Vanda tessellata and the primary molecular targets responsible for its anti-osteoarthritic effects. Therefore, the present study aimed to investigate the potentiality of Vanda tessellata extract (VTE) for treating OA through a comprehensive approach using network pharmacology, molecular docking, molecular dynamic simulation and in-vitro assays. 2. Materials and methods 2.1. Network pharmacology 2.1.1. Identification and screening of drug-like candidates of Vanda tessellata The active chemical constituent of VTE was retrieved the literature [[72]35] and the Indian Medicinal Plants, Phytochemicals and Therapeutics) IMPPAT ([73]https://cb.imsc.res.in/imppat/home) database [[74]36]. The canonical SMILES (Simplified Molecular Input Line Entry System) of each active compound were searched in the PubChem database. The drug likeness parameters were predicted using the Lipinski rule (molecular weight ≤500 g, hydrogen bond acceptor ≤10, hydrogen bond donor ≤5, calculated logP >5 or MlogP >4.15) and the Abbott bioavailability score [[75]37]. The in silico ADME screening and drug-likeness evaluation was carried out using the web tool SwissADME ([76]http://www.swissadme.ch/). The components that did not pass the screening were removed from the list. The filtered compounds were submitted to the Swiss Target Prediction server ([77]http://www.swisstargetprediction.ch/) to identify the possible targets. 2.1.2. Prediction of osteoarthritis-related target genes The keyword “Osteoarthritis” was used to retrieve OA-related target genes. The Gene expression profile datasets [78]GSE55457, [79]GSE12021 and [80]GSE55235 were downloaded from Gene Expression Omnibus database ([81]https://www.ncbi.nlm.nih.gov/geo/). The [82]GSE55457 dataset has 10 normal and 10 OA individuals, the [83]GSE12021 dataset has 9 normal and 10 OA individuals, and the [84]GSE55235 dataset has 10 normal and 10 OA individuals. The Differentially Expressed Genes (DEGs) of all three respective datasets were obtained for the normal and OA patients using iDEP 0.96 online tool. Genes meeting cut-off values p < 0.05 and |log FC| ≥ 0.5 were screened as DEGs between normal and OA individuals. The R package "limma" in iDEP.96 was employed for filtering up and down-regulated DEGs. The iDEP.96 online tool was used for all three datasets to visualise the volcano plots and heatmaps containing the 100 most significant genes each. Four disease databases, namely GeneCards ([85]https://www.genecards.org/), DisGeNET ([86]https://www.disgenet.org/), OMIM ([87]https://omim.org/), and TTD ([88]https://db.idrblab.net/ttd/) were used to retrieve OA-related targets. After eliminating duplicate targets, all the OA-related target genes obtained from the three GEO datasets and disease databases were compiled into a single data file for further analysis. 2.1.3. Constructing a network of compound-disease-target network The Cytoscape v3.9.1 software is used to and visualise the relationship between each gene with active constituents and associated disease. The nodes indicate the herb, candidate active constituents, disease name and common targets, whereas the edges indicate the interconnection between the nodes. The degree of a node indicates the number of interactions it has with the nodes. 2.1.4. Protein-protein interaction of common targets The Venn diagram was constructed to obtain intersecting genes between targets of VTE and OA-associated genes, which might be potential targets for treating OA using the web tool Venny 2.0.2 ([89]https://bioinfogp.cnb.csic.es/tools/venny/). Then, these common genes were uploaded to STRING v11.5 to construct a protein-protein interaction (PPI) network. The species option was set as “Homo sapiens”, and a high threshold score of 0.7 was selected. Disconnected nodes were removed from the network. The network was exported to Cytoscape v3.9.1 for visualisation. CytoHubba, a plugin of Cytoscape, was used to filter topologically important genes based on their ranking methods, i.e degree, maximal clique centrality (MCC) and maximum neighbourhood component (MNC) score. The Top 10 genes on the basis of score were selected as hub genes for subsequent studies. 2.1.5. GO and KEGG pathway enrichment analysis The intersecting common targets were used to analyse Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genome (KEGG) analysis to understand the biological pathways associated with OA. KEGG enrichment analysis was performed using ShinyGo v0.77 ([90]http://bioinformatics.sdstate.edu/go/) web server to identify the metabolic pathway associated with target genes. GO analysis was conducted using the DAVID database ([91]https://david.ncifcrf.gov/) as it gives information about the biological pathways (BP), cellular components (CC) and molecular functions (MF). The visualisation of the enriched GO categories was done using an SR plot in the form of a bar graph. The significant pathways and functions were determined based on the threshold limit of p-value <0.05, FDR (False Discovery Rate) < 0.01 and by selecting the species as “Human”. 2.1.6. Molecular docking The molecular docking was performed between ten hub genes obtained from the PPI analysis and the top four active compounds from the compound-target-disease network. RCSB PDB database was used to download the crystal structure of each of the proteins in the PDB format. The PBD IDs for the proteins having no mutation and low resolution including BCL2 (PBD ID: 7LH7), FGF2 (PBD ID: 6L4O), ICAM 1 (PBD ID: 2OZ4), MAPK1 (PBD ID: 4ZZN), MMP1 (PBD ID: 3SHI), MMP2 (PBD ID: 7XJ0), MMP9 (PBD ID: 4WZV), COX2 (PBD ID: 5F1A), STAT3 (PBD ID: 6NJS) and ESR1 (PBD ID: 1 × 7R) were selected. The active site of the receptor protein with a higher pocket score was selected using the online server Prankweb ([92]https://prankweb.cz). The receptor molecule was processed by removing water molecules, original ligand molecules and addition of non-polar hydrogen with the help of BIOVIA Discovery Studio Visualizer. The Swiss PDB view was used to fill the missing amino acid residues. The final preparation involves the addition of polar hydrogen, gasteiger charges, and energy minimization of each protein, which was done in UCSF Chimera. The two dimensional structure (2D) structure of the core compound was obtained from the PubChem database in SDF format. Then, the refined core proteins and core compounds were uploaded to the PyRx tool for docking. Autodock tool was used to minimise the energy of the ligand molecules and conversion them to pdbqt files. Finally, the docking visualisation was performed using Autodock Vina. The suitable X, Y and Z coordinates were set for each complex in the grid box as per the pocket-1 value given in the Prankweb for active site binding of the ligand. After docking, the binding affinity score was extracted to log*.txt files and the best-docked interactions were again visualized in Discover studio BIOVIA. The validation of docking was confirmed based on low RMSD (<2 Å) of the redocked ligand from the orientation of the co-crystallized ligand and the reproduction of observed interactions from the pdb structure. The best docking position was the one with the minimum root mean square deviation (RMSD) predicted by X-ray crystal configuration, and the binding energy between ligand and receptor protein was evaluated to indicate the binding strength. 2.1.7. Molecular dynamics simulation The molecular dynamic simulation was performed for top-docked complexes using the GROMACS package available in the SiBioLead LLC ([93]https://sibiolead.com). The simulation pre-processing step was done to generate the topology of protein and ligand using AMBERTOOLS and ACPYPE tool kit. The complex was solvated in a triclinic-type periodic boundary box using a SPC water model. Before MD simulation, 0.15 mM of Na^+/Cl^− ions were added to make the system electrically neutral. Next, the energy minimization process was done with a 10000 steps steepest descent algorithm. The simulation system was equilibrated using the constant volume/constant pressure method with an equilibration time of 100 ps. The equilibrium temperature was maintained at 300 K and pressure at 1 bar throughout the simulation. Finally, the MD simulation runs of the complexes was performed with a leap-frog integrator for 100 ns each, and the interpretation of the simulation result was made under GROMACS built-in tools. 2.2. In-vitro analysis 2.2.1. Plant material and extract preparation Fresh aerial roots of Vanda tessellata were collected from Machhaghara, Gajapati district in Odisha (19° 11′ 42.252″ N, 84° 19′ 42.276″ E, 544. 8 m above the sea level) and identified by taxonomist Prof. P.C. Panda. A voucher specimen (2536/CBT Dt. September 25, 2023) was deposited to the herbarium of SOA University. The roots were shade-dried for about 15 days until the moisture was removed and ground in an mechanical homogenizer. The extract was prepared from the dried powder (30 g) with the help of a soxhlet apparatus using methanol (200 ml) at temp. of 25 °C for 48 h. The extract was filtered using Whatmann filter paper and concentrated using a rotary vacuum evaporator (Hei-VAP Core HL, Heidolph, Schwabach, Germany). The extract yield was found to be 13.15 % on dry weight basis. 2.2.2. Cell culture The human synovial sarcoma cell line (SW-982) was purchased from NCCS, Pune, India. Cells were cultured in Dulbecco's Minimum Essential Medium (DMEM) with 10 % FBS, 1 % antibiotic-antimycotic solution and 1 % L-glutamine (200 mM). Cells were maintained in the CO[2] incubator at 37^oC in an atmospheric condition of 5 % CO[2] and 18–20 % O[2]. 2.2.3. Cell viability and migration analysis of SW982 cells The MTT analysis was carried was to check the cytotoxic effect of the drug on the SW982 cell. The SW982 cells were seeded in the 96 well plates at 1 × 10^4 cells/well and treated with VTE of various doses (6.25–200 μg/ml) for 24 h. The MTT solution was added to each well, followed by the replacement of DMSO (100 μl) after 3 h to solubilise the formazan crystals. The percentage of cell proliferation was checked by measuring the absorbance with a microplate reader at 570 nm. Cells were seeded in 6 well plates at a density of 2 × 10^5 cells/well until it attained 80 % confluence as a monolayer. After 24 h of incubation, the monolayer was gently scratched using a sterile pipette tip across the centre of the well. This was done without changing the media to create a wound. After scratching, the well was washed twice with medium to remove the detached cells. Then, the cells were treated with 10 ng/ml of IL-1β followed by different concentrations of VTE (50 and 100 μg/ml) and incubated at 37^oC for 48 h. The cell images were captured at different intervals (0, 24 and 48 h). The microscope (Olympus CKX41-A32PH, Tokyo, Japan) was used for capturing the image of different views of the monolayer. The gap distance can be calculated using Image J software. The rate of wound closure was measured as follows: Wound closure (%) = (Initial wound area at 0 h-Final wound area at t h)/Initial wound area at 0 h*100. 2.2.4. Enzyme-linked immunosorbent assay Quantitative ELISA was carried out to measure he expression levels of the pro-inflammatory cytokines. For this analysis, SW982 cells were cultured in 12-well plated and incubated for 48 h for cell attachment, followed by incubation for 2 h with IL-1β (10 ng/ml) to induce inflammation. Subsequently, the cells were treated with different doses of VTE (50 and 100 μg/ml) and incubated for 24 h. After the treatment, the supernatant was collected, and ELISA was performed to quantify the secretion level of PGE2, IL-6, IL-8 and TNF-α following the manufacturer's kit (RayBiotech, Norcross, GA). 2.2.5. RT-qPCR analysis The SW982 cells were plated in a 6-well plate at a density of 2 × 10^5 cells/well and incubated at 37 °C for 24 h. Then, cells were stimulated with IL-1β (10 ng/ml) and incubated with different concentrations of VTE (6.25 & 25 μg/ml) for 24 h. The RNA was isolated from the cells using a Qiagen RNase kit as per the manufacturer's protocol. The complementary DNA (cDNA) was carried out using an IScript cDNA synthesis kit (Bio-Rad, CA). RT-qPCR reaction was performed on the Quant Studio3 system using SYBR Green Mastermix. The primer sequences used for the expression study are listed in [94]Table 1. The amplification conditions for 40 cycles were as follows: initial denaturation: 95^oC for 5 min, denaturation at 95^oC for 10 s, annealing at 60 for 20 s and extension at 72^oC for 20 s. The relative expression fold change was calculated relative to β-actin using △△Ct method. Table 1. Oligonucleotide sequence of primers used for RT-qPCR analysis. Gene names primer Sequence MMP-2 Forward ACCTGGATGCCGTCGTGGAC Reverse TGTGGCAGCACCAGGGCAGC MMP-9 Forward CAGTACCGAGAGAAAGCCTATT Reverse CAGGATGTCATAGGTCACGTAG β- actin Forward GGAGATTACTGCCCTGGCTCCTA Reverse GACTCATCGTACTCCTGCTTGCTG [95]Open in a new tab 2.2.6. Nuclear translocation assay of NF-κβ SW982 cells (2.5 x 10^5 cells/well) were cultured in a 96-well plate in CO[2] for 24 h. Then, cells were stimulated with IL-1β (10 ng/ml) and kept for 2 h, followed by treatment of VTE (50 & 100 μg/ml). Further, 0.5 mL of BD Cytofix/Cytoperm solution was added, followed by washing with 0.5 % bovine serum albumin (BSA) in 1X phosphate-buffered saline (PBS) and 0.1 % sodium azide. Then, the cells were immunostained with 10 μL of PE Mouse anti–NF– κβ p65 antibody for 30 min and counter-stain with 100 μL of DAPI solution (1 μg/ml) for 10 min. The cells were imaged using the ZEISS LSM 880 imaging system (Carl Zeiss, Oberkochen, Germany) and NF-κβ p65 expression intensity was measured using software Image J. 2.2.7. Detection of MAPK and STAT3 phosphorylation levels by ELISA SW982 cells were grown to confluence in 24-well plates. The cells were stimulated with IL-1β (10 ng/ml) followed by treatment in the absence or presence of VTE (50 & 100 μg/ml) and incubated for 24 h. Subsequently, total and phosphorylated MAPKs (ERK 1/2, p38, and JNK) and STAT3 were detected from cell lysate using a cell-based ELISA kit (RayBiotech, Norcross, GA, USA) following the manufacturer's instructions. 3. Results and discussion 3.1. Network pharmacology 3.1.1. Screening of drug-like candidates of Vanda tessellata and prediction of OA-related target genes Thirteen active components of Vanda tessellata were identified from the literature [[96]35] and the IMPPAT database. All these components successfully passed the ADME screening based on Lipinski's rule of five. This rule use to assess the drug-likeness of a compound, helps in predicting its absorption, distribution, metabolism, and excretion (ADME) properties. Larger molecules tend to have difficulty in crossing the membrane which decreases their bioavailability. High log P value indicates high lipophilicity with poor water solubility resulting in an inefficient proper absorption and distribution. Hydrogen donor or acceptor increases the water solubility but reduces its ability to pass lipid membrane thereby reducing the absorption of the compound [[97]37] ([98]Table 2). Compounds that comply with these rules are more likely to have good oral bioavailability, which is essential for oral medications. By filtering out compounds that are less likely to be orally active, compounds with a higher probability of success can be selected, thereby reducing attrition rates in drug development [[99]38]. Table 2. Drug likeness screening of Vanda tessellata constituent. Sl no. Compounds MW (g/mol) ML OGP HBA<10 HBD<5 Rotational bond<10 TPSA __________________________________________________________________ WLOGP Lipinski violation Abbott Bioavailability (in Å^2) 1 β-Sitosterol 414.71 6.73 1 1 6 20.23 8.02 1 0.55 2 Daucosterol 576.85 3.96 6 4 9 99.38 5.85 1 0.55 3 Heptacosane 380.73 8.86 0 0 24 0 10.78 1 0.55 4 1-Octacosanol 410.76 7.07 1 1 26 20.23 10.14 1 0.55 5 Tessalatin 270.28 1.87 4 2 1 58.92 2.61 0 0.55 6 Gigantol 274.31 2.26 4 2 5 58.92 2.9 0 0.55 7 Gallic acid 170.12 −0.16 5 4 1 97.99 0.5 0 0.56 8 Gamma-Sitosterol 414.71 6.73 1 1 6 20.23 8.02 1 0.55 9 Stigmasterol 412.69 6.62 1 1 5 20.23 7.8 1 0.55 10 Dihydroconiferyl dihydro-p-coumarate 330.37 2.61 5 2 9 75.99 3.22 0 0.55 11 Methyl Linoleate 294.47 4.7 2 0 15 26.3 5.97 1 0.55 12 Syringaldehyde 182.17 0.24 4 1 3 55.76 1.22 0 0.55 13 Vanillin 152.15 0.51 3 1 2 46.53 1.21 0 0.55 [100]Open in a new tab A total of 1443 potential targets corresponding to these active compounds were identified with the canonical SMILES in the Swiss Target Prediction database ([101]Table S1). After removing duplicates a total of 535 targets related to VTE was obtained. The gene expression data of 29 individuals with normal joint conditions and 30 patients with osteoarthritis (OA) were obtained from the GEO dataset ([102]GSE55457, [103]GSE12021, and [104]GSE55235). Differentially expressed genes (DEGs) between individuals with normal joint conditions and those with osteoarthritis (OA) were identified using the iDEP online tool. A log2 (Fold Change) threshold of 0.5 was used as a threshold value for analysing meaningful changes in gene expression. This threshold is suitable when looking for subtle alterations in expression. Several researchers have made use of log FC of 0.5 as a cut off criteria for selecting target genes with significantly differential expression [[105]39,[106]40]. Following a comparative analysis using volcano plots, 119 genes from [107]GSE55457, 134 genes from [108]GSE12021, and 458 genes from [109]GSE55235 were identified to be upregulated in OA patients compared to normal individuals. Additionally, 433 genes from [110]GSE55457, 486 genes from [111]GSE12021, and 455 genes from [112]GSE55235 were found to be downregulated. ([113]Fig. 1A–C and [114]Table S1). The heatmap displaying the top 100 DEGs of all three datasets are demonstrated where red represents normal patients and turquoise blue represents OA patients ([115]Fig. 1D–F). Fig. 1. [116]Fig. 1 [117]Open in a new tab Gene expression levels between normal and Osteoarthritis patients were taken from the Gene Expression Omnibus (GEO) database. Heat map of top 100 up and down-regulated genes from the GEO Dataset (A) [118]GSE55457, (B) [119]GSE12021, and (C) [120]GSE55235. The red and blue colour gradient indicates the up and down regulatory genes. The legend with the color key represents the log fold change of DEGs. Volcano plot distribution for GEO datasets (D) [121]GSE55457 (E) [122]GSE12021 (F) [123]GSE55235. Red indicates highly expressive genes, whereas blue indicates less expressive genes of OA. (For interpretation of the references to colour in this figure legend, the reader is referred to