The Abundant Phytocannabinoids in Rheumatoid Arthritis: Therapeutic Targets and Molecular Processes Identified Using Integrated Bioinformatics and Network Pharmacology

The endocannabinoid system consists of several phytocannabinoids, cannabinoid receptors, and enzymes that aid in numerous steps necessary to manifest any pharmacological activity. It is well known that the endocannabinoid system inhibits the pathogenesis of the inflammatory and autoimmune disease rheumatoid arthritis (RA). To the best of our knowledge, no research has been done that explains the network-pharmacology-based anti-rheumatic processes by focusing on the endocannabinoid system. Therefore, the purpose of this study is to further our understanding of the signaling pathways, associated proteins, and genes underlying RA based on the abundant natural endocannabinoids. The knowledge on how the phytocannabinoids in Cannabis sativa affect the endocannabinoid system was gathered from the literature. SwissTarget prediction and BindingDB databases were used to anticipate the targets for the phytocannabinoids. The genes related to RA were retrieved from the DisGeNET and GeneCards databases. Protein–protein interactions (high confidence > 0.7) were carried out with the aid of the string web server and displayed using Cytoscape. The Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway analysis was used to perform enrichment analyses on the endocannabinoid–RA common targets. ShinyGO 0.76 was used to predict the biological processes listed in the Gene Ontology (GO) classification system. The binding affinity between the ligand and the receptors was precisely understood using molecular docking, induced-fit docking, and a molecular dynamics simulation. The network pharmacology analyses predicted that processes like response to oxygen-containing compounds and peptodyl-amino acid modification are related to the potential mechanisms of treatment for RA. These biological actions are coordinated by cancer, neuroactive ligand–receptor interaction, lipids and atherosclerosis, the calcium signaling pathway, and the Rap1 signaling pathway. According to the results of molecular docking, in the context of RA, phytocannabinoids may bind to important target proteins such PIK3CA, AKT1, MAPK9, PRKCD, BRAF, IGF1R, and NOS3. This entire study predicted the phytocannabinoids’ systemic biological characteristics. Future experimental research is needed, however, to confirm the results so far.


Introduction
Although inflammation is a typical reaction to tissue damage, it can become out of control and cause further problems. Specific soluble pro-inflammatory mediators, including cytokines or chemokines, prostaglandins, and leukotrienes, are produced as a result of an inflammatory reaction [1]. Many acute and chronic disorders are made worse by persistent inflammation. Rheumatoid arthritis (RA) is just one of the many chronic inflammatory

Identification of Bioactive Compounds and Targets
The abundant phytocannabinoids of Cannabis sativa used in this investigation were derived according to reports by Morales et al. [14]. In addition, the canonical forms and 2D structures of the phytocannabinoids were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/, accessed on 1 December 2022). The targets of the phytocannabinoids were found using the web servers Swiss Target Prediction and Bind-ingDB (https://www.bindingdb.org/bind/index.jsp, accessed on 1 December 2022). In addition, the canonical forms and 2D structures of the phytocannabinoids were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov, accessed on 1 December 2022). The RA targets were derived from two databases: DisGeNET (www.disgenet.org, accessed on 1 December 2022) and GeneCards (www.genecards.org, accessed on 1 December 2022). The UniProt database (www.uniprot.org, accessed on 1 December 2022) was utilized to standardize the target names into official gene symbols. The keywords "Rheumatoid arthritis" and the organism "Homo sapiens" were used to seek the targets, which are verified by UniProt ID. To obtain the Venn diagram of the shared targets, Venny2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/index.html, accessed on 1 December 2022) was used.

Collection of Protein-Protein Interaction (PPI) Data
The chosen targets were added to the STRING 11.0 database (available at https: //string-db.org, accessed on 1 December 2022) to create a PPI network using "Homo sapiens" as the species and "highest confidence > 0.7" as the minimal interaction threshold. Additionally, the detached nodes were kept a secret. Cytoscape software was used to visualize these networks. To identify PPI subnetworks and clusters with a high degree of interconnection within the PPI network, we employed the Molecular Complex Detection (MCODE) method (version 1.5.1). As threshold settings, we set the Cytoscape plug-in MCODE's maximum depth to 100, node score at 0.2, and K-score to 2 [21].

GO and KEGG Enrichment Analyses
The KEGG metabolic pathway analysis was used to perform enrichment analyses on C. sativa and RA mutual targets, and ShinyGO 0.76 (available at http://bioinformatics. sdstate.edu/go76/, accessed on 1 December 2022) was used to predict Gene Ontology (GO) biological processes. The top 20 biological pathways and processes (False Discovery Rate cut-off of 0.05) were selected based on fold enrichment, and they were then displayed as bubble charts with the adjusted p-values, fold enrichment, gene counts, etc. The GO analysis predicted the various biological processes, cellular components, and molecular functions. A network was constructed using the key pathways obtained from the KEGG analysis for the discovery of the important target proteins that play a significant role in the treatment effects of C. sativa [22]. For the purpose of predicting the likely mechanism of action of the bioactive components of C. sativa in RA, all the information on the important targets and associated pathways was gathered through online software.

Molecular Docking Validation
Molecular docking studies of 13 phytocannabinoids were done using Glide XP (extraprecision) implemented in the Schrödinger suite. The energy-minimized structures of the ligands were docked with the crystal structure of PIK3CA (PDB ID: 4L23, resolution of 2.50 Å), AKT1 ((Gene: Akt1) (AlphFold)), MAPK9 ((Gene: MAPK9) (AlphFold)), PRKCD ((Gene: PRKCD) (AlphFold)), BRAF (PDB ID: 6P3D), IGF1R (PDB ID: 5FXS), and NOS3 (PDB ID: 6NH5) with validated grid parameters [23]. The RCSB protein data bank (PDB) was used to download the X-ray crystal coordinates of different proteins. From the Al-phaFold, the predicted structures of the proteins AKT1, MAPK9, and PRKCD were obtained. We used the protein preparation wizard to prepare the proteins. The water molecules that lacked three hydrogen bonds were eliminated. Then, considering the correct ionization states for both basic and acidic amino acid residues, hydrogen bonds corresponding to pH 7 were added. The force field OPLS3 (Optimized Potentials for Liquid Simulations 3) was used to minimize the energy of the protein 3D structure. The optimum protein binding site was identified using Sitemap analysis for proteins predicted by AlphaFold. Additionally, the best binding site's center was chosen to generate the grid box. The centroid of each of the co-crystallized ligands for PIK3CA, BRAF, IGF1R, and NOS3 (X6K, 0LI, OZN, and KLD, respectively) was chosen to create the grid box. Finally, docking of all low-energy conformations into the Sitemap-predicted site was done using Glide in XP mode in the absence of any constraints.

Molecular Dynamics (MD) Simulation
The Desmond module of the Schrödinger software was used to simulate the docked protein-ligand complexes in an explicit water environment. Using a system builder panel, an orthorhombic box was constructed. The simulation of the bound ligand-protein complex was performed using a simple point charge (SPC) model of water. The MD simulation was performed using a constant-temperature, constant-pressure ensemble (NPT) at 310 K and 1.013 bar as an atmospheric pressure intended to last 100 ns. The Simulation Interactions Diagram Report of the Desmond software was used to thoroughly examine the results of the molecular dynamics research [24,25].

Induced-Fit Docking (IFD)
Induced-fit docking was carried out on the best binding energy-containing bioactives with their respective receptors by employing the IFD panel of the Schrödinger suite [26]. Here, the IFD computations were carried out using the OPLS3 force field and the initial sample methodology, which generated up to 20 docked complexes utilizing automatic docking settings. The best poses produced during the standard Glide XP (extra-precision) docking of bioactive substances with their corresponding receptors were used as the input structures for IFD, and the centroid of the ligands was chosen to generate the grid box. No additional constraints were applied. A ring-conformational sampling of the ligand was done around a 2.5 kcal/mol energy window. Initial-glide docking was performed by setting both the receptor and ligand's van der Waals (vdW) scaling to 0.5. The amino acid residues were also cleaned up using the Schrödinger Prime module within 5 Å of the ligand. The Prime-refined receptor-ligand complexes within 30 kcal/mol were then redocked using the Glide XP mode. The default settings for other parameters were maintained [27].

Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) Calculation
The binding free energy was calculated from the structural information of the MM-GBSA method. These approaches employ the solvent accessibility method, molecular mechanics, and the generalized Born model to evade the free energy simulations' convolution [28]. The energy difference between the bound and unbound complex of the protein-ligand complex was calculated with the help of MM-GBSA. The binding free energies were computed by the equation as follows: Here, the MM-GBSA calculation [29,30] was done by considering the VSGB solvation model with OPLS3 Force Field with a minimized sampling method.

Physicochemical and ADME/T Studies
The QikProp module of Schrödinger suite (QikProp, 2017_2) [31] and SwissADME [32,33] webserver (available at http://www.swissadme.ch/, accessed on 3 December 2022) were used to determine the physicochemical and pharmacokinetic properties of all the phytocannabinoids, which helps to determine the physicochemical significant descriptors and pharmacokinetically important properties of the phytocannabinoids.

Identification of Bioactive Compounds and Targets
A total of 13 phytocannabinoids, found in Cannabis sativa, were retrieved from the literature [14]. A total of 140 targets of the various phytocannabinoids were identified using the BindingDB (Supplementary Table S1) and SwissTargetPrediction (Supplementary Table S2) databases. Additionally, it was discovered that there were 2723 and 4881 RA-related targets predicted from two distinct databases, namely the GeneCards and DisGeNET databases (Supplementary Tables S3 and S4). After eliminating the repeating gene symbols, 384 targets linked to phytocannabinoids and 5728 targets associated to RA were obtained from these databases. The intersection of targets linked to phytocannabinoids' bioactive properties and targets for RA yielded 244 similar targets (Supplementary Table S5) ( Figure 1A).
were used to determine the physicochemical and pharmacokinetic properties of all the phytocannabinoids, which helps to determine the physicochemical significant descriptors and pharmacokinetically important properties of the phytocannabinoids.

Identification of Bioactive Compounds and Targets
A total of 13 phytocannabinoids, found in Cannabis sativa, were retrieved from the literature [14]. A total of 140 targets of the various phytocannabinoids were identified using the BindingDB (Supplementary Table S1) and SwissTargetPrediction (Supplementary Table S2) databases. Additionally, it was discovered that there were 2723 and 4881 RArelated targets predicted from two distinct databases, namely the GeneCards and Dis-GeNET databases (Supplementary Tables S3 and S4). After eliminating the repeating gene symbols, 384 targets linked to phytocannabinoids and 5728 targets associated to RA were obtained from these databases. The intersection of targets linked to phytocannabinoids' bioactive properties and targets for RA yielded 244 similar targets (Supplementary Table  S5) ( Figure 1A).

Phytocannabinoids-Rheumatoid Arthritis Target and PPI Network Analysis
The phytocannabinoids and common targets were added to Cytoscape 3.9.1 for the construction of a network ( Figure 1B,C). This constructed network demonstrated the existence of many and intricate impacts of phytocannabinoids on RA. Based on the gradation, the utmost common targets are mentioned in the Supplementary Material (Supplementary Table S6). The results predicted that the most prominent RA-associated genes were PIK3CA, PIK3CB, AKT1, AKT2, MAPK8, MAPK9, MAPK10, PRKCB, MAPK14, IGF1R, BRAF, PRKCZ, PRKCA, and NOS3 due to interaction of phytocannabinoids. The results of module analysis of the MCODE also predicted that the involvement of the above

Phytocannabinoids-Rheumatoid Arthritis Target and PPI Network Analysis
The phytocannabinoids and common targets were added to Cytoscape 3.9.1 for the construction of a network ( Figure 1B,C). This constructed network demonstrated the existence of many and intricate impacts of phytocannabinoids on RA. Based on the gradation, the utmost common targets are mentioned in the Supplementary Material (Supplementary Table S6). The results predicted that the most prominent RA-associated genes were PIK3CA, PIK3CB, AKT1, AKT2, MAPK8, MAPK9, MAPK10, PRKCB, MAPK14, IGF1R, BRAF, PRKCZ, PRKCA, and NOS3 due to interaction of phytocannabinoids. The results of module analysis of the MCODE also predicted that the involvement of the above genes which supports the present findings (Figure 2A-C). The phytocannabinoids-RA target network constructed in Cytoscape ( Figure 3A) predicted that there were 244 common target genes between 13 phytocannabinoids and RA. Further, 20 nodes of the common genes were found to be majorly involved as per the investigation of the topological network in the Cytohubba plug-in ( Figure 3B). genes which supports the present findings (Figure 2A-C). The phytocannabinoids-RA target network constructed in Cytoscape ( Figure 3A) predicted that there were 244 common target genes between 13 phytocannabinoids and RA. Further, 20 nodes of the common genes were found to be majorly involved as per the investigation of the topological network in the Cytohubba plug-in ( Figure 3B).  genes which supports the present findings (Figure 2A-C). The phytocannabinoids-RA target network constructed in Cytoscape ( Figure 3A) predicted that there were 244 common target genes between 13 phytocannabinoids and RA. Further, 20 nodes of the common genes were found to be majorly involved as per the investigation of the topological network in the Cytohubba plug-in ( Figure 3B).

GO and KEGG Enrichment Analyses
A total of 4255 terms linked to the therapeutic effects of phytocannabinoids on RA were found as a result of GO enrichment analysis on common phytocannabinoids-RA targets. Figure 4A-C displays bubble charts for the top 20 terms in each of the three categories-biological process, cellular component, and molecular function. The results of the GO analysis showed that phytocannabinoids mainly suppressed cellular oxidation, protein phosphorylation, and changes to peptidyl amino acids. The top 20 pathways, which were ranked by the number of genes engaged in the expressed pathways ( Figure 4D-F), were identified by the KEGG pathway enrichment analysis and were found to be involved in the interactions of phytocannabinoids to treat RA (Table S7). The findings indicated that the primary pathways involved in the current study were those that were associated with cancer, neuroactive ligand-receptor interaction, lipids and atherosclerosis, the calcium signaling pathway, and the Rap1 signaling pathway. These pathways are primarily associated with angiogenesis and inflammation, which are the defining characteristics of RA. By selecting the top 10 signaling pathways based on the number of genes, this study was able to identify the main primary targets of phytocannabinoids that are responsible for the pathogenesis of RA. Additionally, PIK3CA, AKT1, MAPK9, PRKCD, BRAF, IGF1R, and NOS3 were among the common target genes that were involved in reducing the severity of RA when different phytocannabinoids were used as treatment. These genes are also known to play a role in the pathogenesis of RA as shown by Cytohubba's visualization of these genes in Cytoscape.

GO and KEGG Enrichment Analyses
A total of 4255 terms linked to the therapeutic effects of phytocannabinoids on RA were found as a result of GO enrichment analysis on common phytocannabinoids-RA targets. Figure 4A-C displays bubble charts for the top 20 terms in each of the three categories-biological process, cellular component, and molecular function. The results of the GO analysis showed that phytocannabinoids mainly suppressed cellular oxidation, protein phosphorylation, and changes to peptidyl amino acids. The top 20 pathways, which were ranked by the number of genes engaged in the expressed pathways ( Figure 4D-F), were identified by the KEGG pathway enrichment analysis and were found to be involved in the interactions of phytocannabinoids to treat RA (Table S7). The findings indicated that the primary pathways involved in the current study were those that were associated with cancer, neuroactive ligand-receptor interaction, lipids and atherosclerosis, the calcium signaling pathway, and the Rap1 signaling pathway. These pathways are primarily associated with angiogenesis and inflammation, which are the defining characteristics of RA. By selecting the top 10 signaling pathways based on the number of genes, this study was able to identify the main primary targets of phytocannabinoids that are responsible for the pathogenesis of RA. Additionally, PIK3CA, AKT1, MAPK9, PRKCD, BRAF, IGF1R, and NOS3 were among the common target genes that were involved in reducing the severity of RA when different phytocannabinoids were used as treatment. These genes are also known to play a role in the pathogenesis of RA as shown by Cytohubba's visualization of these genes in Cytoscape.

Molecular Dynamics
The molecular dynamics simulations of protein-ligand complexes are thought to be the best developed method for assessing the structure-to-function relationships of macromolecules in computer-aided drug design. To clarify the cannabinoid's potency against the six proteins associated with RA, an impartial molecular dynamics investigation was conducted in the current study. The stability and viability of the best docked ligands against all the target proteins were examined using MD simulation studies, which were run for 100 ns. The root mean square deviation (RMSD), the root mean square fluctuation (RMSF), and protein-ligand interactions were generated with the production of 100 ns MD simulations, and their relevance in the solidity of the protein-ligand complexes was thoroughly examined. Protein-ligand interactions and protein-ligand RMSDs were investigated and visualized during this period.
For the PIK3CA-cannabigerol complex, their RMSD for the first 50 nanoseconds was more than 2.5 Å, and from 50 nanoseconds on, the ligand started to stabilize with the protein at above 3.5 Å. It was discovered through protein-ligand contact analysis that roughly 50% of the amino acids were responsible for the hydrophobic interactions with cannabigerol. The remaining 50% were engaged in either the H-bond interactions or the water bridges. Hydrophobic interactions with cannabigerol were largely caused by the residues Met772, Pro778, Trp780, Ile800, Tyr836, Ile848, Val850, Met922, and Ile932, of which Met772 and Ile932 were the most accountable at 0.15 ns and 0.20 ns fractions of the interaction time, respectively. Glu849 and Val851 were responsible for the highest H-bond contacts at 0.10 ns and 0.25 ns fractions, while Met772, Ser919, and Asn920 were responsible for forming highest water bridge contacts at 0.20 ns and 0.35 ns fractions.
For the AKT1-cannabitriol complex, their RMSD for the first 60 nanoseconds was stabilized at above 4 Å, while for the next 40 nanoseconds, both the ligand and protein RMSDs were slightly fluctuating at 4.8 Å. Through protein-ligand contact analysis, it was discovered that the residues Phe161, Val164, Ala177, Lys179, and Leu181 were highly responsible for forming hydrophobic contact with cannabitriol, among which Val164 and Lys179 were the most accountable ones at 0.3 ns and 0.6 ns fractions of the interaction time, respectively. Glu234 and Asp439 were responsible for the highest water bridge forming contacts at 0.3 ns and 0.5 ns fractions. Lys179 and Asp439 were the H-bonding residues, whereas Glu234 and Asp439 were the residues forming ionic contacts.
The protein-ligand RMSD plots for the MAPK9-cannabinodiol complex showed that RMSD was greater than 4 Å for the first few nanoseconds, increased to 12 Å after a while, and then increased to 16 Å before stabilizing there. Through protein-ligand contact analysis, it was found that the residues Val40, Ala53, Leu57, Ile86, Met108, Leu110, Val158, and Leu168 were responsible for hydrophobic contacts with cannabinodiol, among which Val40, Ala53, Val158, and Leu168 were the most accountable ones at 0.40 ns, 0.30 ns, 0.35 ns, and 0.30 ns fractions of the interaction time, respectively. Ile32 and Lys55 were responsible for the highest water bridging at 0.35 ns, while, on the other hand, Ile32, Lys55, Asn114, Ser155, and Asn156 were responsible for H-bonding contacts.
The RMSD for the PRKCD-cannabitriol complex was greater than 4 Å for the first few nanoseconds before increasing to 8 Å and then stabilizing there. According to the protein-ligand contact study, Ser173, Gln284, and Lys285 in H-bond interactions at 0.8, 1.0, and 0.8 ns of the interaction period, respectively, accounted for the majority of the protein-ligand connections with cannabitriol. The remaining amino acid residues were in charge of either creating water bridges or creating hydrophobic interactions.
The RMSD for the BRAF-cannabigerol complex was greater than 1.6 Å for the first few nanoseconds before increasing to 2.4 Å, then 2.8 Å, and then stabilizing at that value. Trp531, Phe583, and Phe595 were determined to be primarily in charge of hydrophobic interactions with cannabigerol at 0.22 ns, 0.15 ns, and 0.14 ns fractions of the interaction time, respectively, according to the protein-ligand contact study. The majority of the residues either contributed to the formation of water bridges or H-bonds. Among these, Gln461 and Cys532 played a significant role in the formation of H-bonds. Ile463, Glu533, and Gly534 played a significant role in the formation of water bridges.
The RMSD for the IGF1R-∆ 9 -tetrahydrocannabinol complex was more than 3.0 Å for the first few nanoseconds, but with time, it increased to 4.8 Å, where it stabilized. Based on protein-ligand contact analysis, it was found that Met1082 was responsible for the highest H-bond contact with ∆ 9 -tetrahydrocannabinol at the 1.75 ns fraction of the interaction time. The amino acid residues Leu1005, Val1013, Ala1031, Met1142, and Ile1160 were mostly responsible for hydrophobic contacts. Leu1005 and Thr1157 were responsible for forming water bridges (Figures 6 and S2).
The RMSD for the NOS3-cannabinol complex was higher than 2.4 Å for the first few nanoseconds, then it increased to 2.8 Å, and ultimately it decreased to 2.4 Å and stabilized there. The amino acid residues Trp178 and Phe353 were predicted to have the most hydrophobic interactions with cannabinol based on the protein-ligand contact study. Most of the water bridging was done by Gly186 and Val185. It was Cys184 that made the H-bond interaction. Figure S3) of all the proteins with their best docked pose-forming phytocannabinoids (obtained via standard XP-mode docking) predicted that almost all the proteins had improved Glidescores, with the exception of PIK3CA and PRKCD (Supplementary Table S10). However, for IFD of PIK3CA with cannabigerol, the number of H-bonds increased from one (in the normal XP-mode Glide docking) to two (in the IFD-generated pose). From IFD of MAPK9 with cannabinodiol it was found that the number of H-bonds decreased from two (in the normal XP-mode Glide docking) to one (in the IFD-generated pose). For IFD of PRKCD with cannabitriol, the number of H-bonds increased from three (in the normal XP-mode Glide docking) to five (in the IFD-generated pose). For IFD of BRAF with cannabigerol, the number of H-bonds increased from one (in the normal XP-mode Glide docking) to two (in IFD-generated pose). Similarly, from IFD of IFG1R with ∆ 9 -tetrahydrocannabinol it was found that the number of H-bonds increased from one (in the normal XP-mode Glide docking) to two (in IFD-generated pose). In case of IFD of NOS3 with cannabinol, the number of π-π stacking interactions increased from one (in the normal XP-mode Glide docking) to two (in IFD-generated pose). Additionally, one H-bond interaction with Ala177 was observed from the IFD-generated pose which was absent in the normal XP-mode Glide docked pose. Furthermore, the prime MM-GBSA binding free energies of the best five docked complexes were calculated and are given in the Supplementary Material (Table S11).

Physicochemical and ADME/T Studies
Data generated by QikProp predicted that all the expected attributes for phytocannabinoids fell inside the acceptable range. The SwissADME-predicted data showed that the phytocannabinoids had significant values for the studied attributes, and based on their physicochemical features, they demonstrated drug-like qualities. All phytocannabinoids obeyed Veber's rule. The phytocannabinoids cannabigerol and cannabichromene failed to pass/satisfy the blood-brain barrier (BBB) filter ( Figure S4), but other phytocannabinoids, such as cannabinol, cannabivarin, cannadinodiol, and cannabitriol, were discovered to be P-glycoprotein (PGP) substrates. It was projected that the majority of phytocannabinoids, except for cannabielsoin and cannabitriol, would be CYP450 enzyme inhibitors. Both CYP2C9 and CYP2D6 subtypes were found to be inhibited by cannabielsoin; however only CYP2D6 subtypes were found to be inhibited by cannabitriol (Supplementary Table S12).

Discussion
The chronic inflammatory condition known as arthritis is characterized by immunemediated inflammatory sinusitis, which causes hyperplasia, vasculogenesis, pannus development, the deterioration of cartilage and bone, joint deformity, and functional impairment. Many phytocannabinoids, cannabinoid receptors, and enzymes that are part of the endocannabinoid system all contribute to the pharmacological activity of cannabis in different ways. It is well known that the endocannabinoid system inhibits the pathogenesis of the inflammatory and autoimmune disease RA. Due to the abundance of natural phytocannabinoids, the goal of this study was to build on the signaling pathways, associated proteins, and genes that cause RA using integrated bioinformatics and network pharmacology. The general workflow scheme of the present study is depicted in Figure 7. of the inflammatory and autoimmune disease RA. Due to the abundance of natural phytocannabinoids, the goal of this study was to build on the signaling pathways, associated proteins, and genes that cause RA using integrated bioinformatics and network pharmacology. The general workflow scheme of the present study is depicted in Figure 7. We gathered 2723 and 4881 predicted targets for phytocannabinoids from two distinct databases, namely GeneCards and the DisGeNET database, respectively. Thirteen phytocannabinoids and 244 common targets were obtained by the intersection of RA targets and phytocannabinoids (with the aid of venny 2.1.0). The PPI enrichment analysis results were reported with 244 nodes, 699 edges, highest confidence (0.90), and the PPI enrichment p-value of <1.0 × 10 −16 . As would be predicted, the proteins interacted with one another quite a bit. The results of the enrichment analysis showed that the proteins were, collectively, at least loosely physiologically interconnected. For the hub PPI, we chose nodes whose median edge count was greater than double that of all other nodes.
The results of the PPI network and GO analysis predicted that the key proteins involved in the numerous phytocannabinoids-RA common targets were primarily those involved in the pathways for cancer, neuroactive ligand-receptor interaction, lipids and atherosclerosis, calcium signaling, and RAP signaling. In cancer, there is an abnormality in the signaling of PI3K/Akt/mTOR due to which many types of research have been done to find the inhibitors of PI3K/Akt/mTOR signaling which may act as anticancer agents. RA is an autoimmune-mediated inflammatory disease where the abnormal proliferation of immune cells, macrophages, monocytes, dendritic cells, and synovial fibroblasts may significantly overlap with the abnormal growth of cancer cells. The results of some recent studies in arthritis using PI3K signaling inhibitors suggest that small molecule inhibitor strategies directed at PI3K signaling may be useful for immune-mediated arthritis [34]. The present study predicted that PIK3CA and AKT are the most common genes found in the top pathways of phytocannabinoids to combat arthritis. A previous study suggested that phytocannabinoids are responsible for the regulation of PI3K/AKT/GSK-3 in the brain [35]. Hence, PI3K and AKT may be the target of phytocannabinoids in RA. A previous study revealed that MAPK9 causes an increase in cell proliferation, migration, invasion, inflammatory response, facilitates apoptosis in RA patients, and promotes the progression of fibroblast-like synoviocytes [36]. Protein kinase C like PRKCD is known to have a cru- We gathered 2723 and 4881 predicted targets for phytocannabinoids from two distinct databases, namely GeneCards and the DisGeNET database, respectively. Thirteen phytocannabinoids and 244 common targets were obtained by the intersection of RA targets and phytocannabinoids (with the aid of venny 2.1.0). The PPI enrichment analysis results were reported with 244 nodes, 699 edges, highest confidence (0.90), and the PPI enrichment p-value of <1.0 × 10 −16 . As would be predicted, the proteins interacted with one another quite a bit. The results of the enrichment analysis showed that the proteins were, collectively, at least loosely physiologically interconnected. For the hub PPI, we chose nodes whose median edge count was greater than double that of all other nodes.
The results of the PPI network and GO analysis predicted that the key proteins involved in the numerous phytocannabinoids-RA common targets were primarily those involved in the pathways for cancer, neuroactive ligand-receptor interaction, lipids and atherosclerosis, calcium signaling, and RAP signaling. In cancer, there is an abnormality in the signaling of PI3K/Akt/mTOR due to which many types of research have been done to find the inhibitors of PI3K/Akt/mTOR signaling which may act as anticancer agents. RA is an autoimmune-mediated inflammatory disease where the abnormal proliferation of immune cells, macrophages, monocytes, dendritic cells, and synovial fibroblasts may significantly overlap with the abnormal growth of cancer cells. The results of some recent studies in arthritis using PI3K signaling inhibitors suggest that small molecule inhibitor strategies directed at PI3K signaling may be useful for immune-mediated arthritis [34]. The present study predicted that PIK3CA and AKT are the most common genes found in the top pathways of phytocannabinoids to combat arthritis. A previous study suggested that phytocannabinoids are responsible for the regulation of PI3K/AKT/GSK-3 in the brain [35]. Hence, PI3K and AKT may be the target of phytocannabinoids in RA. A previous study revealed that MAPK9 causes an increase in cell proliferation, migration, invasion, inflammatory response, facilitates apoptosis in RA patients, and promotes the progression of fibroblast-like synoviocytes [36]. Protein kinase C like PRKCD is known to have a crucial role in the pathogenesis of autoimmune diseases like RA [37]. The current study suggested the involvement of MAPK and the protein kinase C pathway in combating arthritis by phytocannabinoids. Previous in vitro studies in cell lines reveal that phytocannabinoids cause cytotoxicity in leukemic cells through the MAPK pathway [38]. The cannabinoid receptor 2 also has an inhibitory role in synovial fibroblasts in RA through protein kinase [39]. These results support our current findings. The BRAF gene was found to transform the rheumatoid synovial fibroblasts and lead to the destruction of articular cartilage and bone in RA [40]. Previous studies revealed that IGF-1R signaling leads to inflammation in RA due to the formation of various cytokines like IL-6 and activation of T cells [41]. A previous study revealed that cannabinoid receptors promote the chronic intermittent hypoxia-induced breast cancer metastasis via IGF-1R/AKT/GSK-3β. It indicates that phytocannabinoids may ameliorate RA through cannabinoid-receptor-induced IGF-1R and AKT. There is a common link of pathogenesis of RA and atherosclerosis. Patients with RA are at approximately 1.5-fold risk of atherosclerotic cardiovascular disease (CVD) compared with the general population [42]. Pro-inflammatory cytokines such as tumor necrosis factor alpha and interleukin-6, involved in the pathogenesis of RA, are also independently predictive of subsequent cardiovascular disease (CVD). In RA, inflammation alters HDL constituents and the concentration of LDL and HDL, thus facilitating atherosclerosis and CVD events. On the other hand, the increase of oxidative processes, frequently observed in RA, also induces atherosclerosis [43]. Cannabinoids were reported to cause a decrease in inflammation during atherosclerosis [44] which correlates with the current findings. A previous study reported that the alterations in calcium signaling may cause the functional T cell abnormalities associated with RA [45,46]. The RAP signaling pathway is also responsible for the occurrence of RA. RAP1 is deregulated in T lymphocytes, which leads to the overproduction of reactive oxygen species, which promotes the hyporesponsiveness of synovial T lymphocytes, leading to inflammation [47]. A previous study revealed that the cannabinoid receptor (CB1) causes neurite outgrowth through the Rap1 pathway [48]. This indicates that phytocannabinoids may alter the RAP1 pathway through the cannabinoid receptor (CB1). Nitric oxide is produced by nitric oxide synthetase (NOS). Important cytokines produced in RA are IL-1 and TNF-α, which stimulate NOS in inflammatory cells and increase the production of nitric oxides. There is an increase in nitric oxide levels in arthritis, as shown in preclinical and clinical studies [49]. A previous study reported that phytocannabinoids may prevent RA conditions by decreasing nitric oxide levels and cytokines [50] which supports the current findings. The results of the present studies correlate with the results of the previous similar studies that substantiate the present findings [51,52].
In light of the current findings of network pharmacology, PIK3CA, AKT1, MAPK9, PRKCD, BRAF, IGF1R, and NOS3 were chosen as the seven key target proteins closely related to the progression of RA. These proteins were further subjected to molecular docking with phytocannabinoids. Negative binding energy indicates effective binding to the target. Lower binding energy suggests a stable ligand-receptor complex better explained as <−5 kcal/mol called good binding activity and <−7 kcal/mol called strong binding activity. Except for AKT1 and PRKCD, the Glidescores of the top three compounds against the other five proteins were in the range of −8.888 to −6.751 kcal/mol. The top three active ingredients were shown to have a strong binding interaction with their respective proteins. It entails unexpected bias from the materiality, which could hasten the anonymous fallacy in in vitro or in vivo studies. However, the results of the molecular docking experiment suggest the potential therapeutic approaches and could guide trials in animals. Network pharmacology is currently used extensively for various traditional medicines, such as Indian traditional medicine, traditional Chinese medicine, Kampo, traditional Korean medicine, etc. Network pharmacology research is being conducted owing to its robust potential to investigate convoluted relationships among multiple targets. Additionally, network pharmacology accelerates the modernization of the Indian traditional system of medicine by usefully illuminating the synergistic effects of biochemical and systemic biology methodologies. Researchers may identify prospective active ingredients for upcoming in vitro and in vivo studies by using molecular docking to lay out an effective binding mode of the active ingredients with disease-related important target proteins. However, a number of the Ayurvedic formula's constituents can have antagonistic and synergistic effects on the body while being treated because they target the same protein. It is important to consider these implications and to continue to assess them. The findings of MM-GBSA (binding free energy) and IFD support the molecular dynamics findings that the ligand is stable within the protein. One of the finest molecular dynamics results was obtained with the Akt1-cannabitriol complex, despite its Glidescore of −6.035 kcal/mol. Among the seven proteins, the phytocannabinoids cannabigerol, canabichromene, and ∆ 9 -tetrahydrocannabivarin had the highest Glidescores and the most drug-like qualities, as well as acceptable ADME/T profiles.
The present study evaluated the integrated bioinformatics and network pharmacology to identify therapeutic targets and molecular processes involving the endocannabinoid system in patients RA. The work used a combined bioinformatics and network pharmacology approach to further examine signaling pathways, associated proteins, and genes that may implicate naturally occurring phytocannabinoids in the pathogenesis of RA. The current research findings suggest that the abundant phytocannabinoids can exert their curative effect through a synergistic effect of multi-components, multi-targets, and multi-pathways in the treatment of RA. However, due to the limitation of network pharmacology and molecular docking, a more complete evaluation of the role of the endocannabinoid system in the immune landscape in patients with RA and a comprehensive analysis of patients with RA and normal controls would be necessary. Hence, further evaluation of in vitro and in vivo studies needs to be carried out to validate the current findings, which is the future scope of the present study. The probable mechanisms of phytocannabinoids in rheumatoid arthritis are depicted in Figure 8. medicine by usefully illuminating the synergistic effects of biochemical and systemic biology methodologies. Researchers may identify prospective active ingredients for upcoming in vitro and in vivo studies by using molecular docking to lay out an effective binding mode of the active ingredients with disease-related important target proteins. However, a number of the Ayurvedic formula's constituents can have antagonistic and synergistic effects on the body while being treated because they target the same protein. It is important to consider these implications and to continue to assess them. The findings of MM-GBSA (binding free energy) and IFD support the molecular dynamics findings that the ligand is stable within the protein. One of the finest molecular dynamics results was obtained with the Akt1-cannabitriol complex, despite its Glidescore of −6.035 kcal/mol. Among the seven proteins, the phytocannabinoids cannabigerol, canabichromene, and Δ 9tetrahydrocannabivarin had the highest Glidescores and the most drug-like qualities, as well as acceptable ADME/T profiles. The present study evaluated the integrated bioinformatics and network pharmacology to identify therapeutic targets and molecular processes involving the endocannabinoid system in patients RA. The work used a combined bioinformatics and network pharmacology approach to further examine signaling pathways, associated proteins, and genes that may implicate naturally occurring phytocannabinoids in the pathogenesis of RA. The current research findings suggest that the abundant phytocannabinoids can exert their curative effect through a synergistic effect of multi-components, multi-targets, and multi-pathways in the treatment of RA. However, due to the limitation of network pharmacology and molecular docking, a more complete evaluation of the role of the endocannabinoid system in the immune landscape in patients with RA and a comprehensive analysis of patients with RA and normal controls would be necessary. Hence, further evaluation of in vitro and in vivo studies needs to be carried out to validate the current findings, which is the future scope of the present study. The probable mechanisms of phytocannabinoids in rheumatoid arthritis are depicted in Figure 8.

Conclusions
The network pharmacology analysis results predicted that the processes such as response to oxygen-containing compounds and peptodyl-amino acid modification are related to the potential treatment mechanisms for rheumatoid arthritis. These biological actions may be coordinated by cancer, neuroactive ligand-receptor interaction, lipids and atherosclerosis, the calcium signaling pathway, and the Rap1 signaling pathway. According to the results of molecular docking, in the context of rheumatoid arthritis, phytocannabinoids may bind to important target proteins such as PIK3CA, AKT1, MAPK9, PRKCD, BRAF, IGF1R, and NOS3. Due to a lack of experimental validation, this study has shortcomings, although it can inspire in vitro or in vivo studies. However, due to the limitations of network pharmacology and molecular docking, a more complete evaluation of the role of the endocannabinoid system in the immune landscape in patients with RA and a comprehensive analysis of patients with RA and normal controls would be necessary. This entire study predicted the phytocannabinoids' systemic biological characteristics. Future experimental research is still needed to confirm the results so far.
Author Contributions: Y.N.D. and K.K.R. contributed to the generation of the concept and supervision of the research work. A.N. contributed significantly in the network analysis and drafting the research paper. A.D. assisted in the molecular docking and dynamics analysis. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The datasets supposing the current study are available in public databases from IMPPAT, PUBMED, PubChem, Swiss Target Prediction, BindingDB, DisGeNET, GeneCards, OMIM, UniProt, STRING, ShinyGO, RCSB PDB, and SwissADME. All the analysed data during the current study are available in the manuscript and Supplementary Material.