Investigation of the mechanism of action of chemical constituents of celery seed against gout disease using network pharmacology, molecular docking, and molecular dynamics simulations

Abstract Celery (Apium graveolens L.) has long been considered as a potential herbal medicine for the prevention and treatment of gout. However, the relationship between the chemical constituents and pharmacological activities of this medicinal plant has not been fully investigated yet. Therefore, this study aims to apply network pharmacology, molecular docking and molecular dynamics to explore the relationship between the chemical constituents of celery seed and its biological effects in the treatment of gout. Network pharmacology was built and analyzed based on the data collected from GeneCards, OMIM databases and SwissTargetPrediction web server using Cytoscape 3.9.0 software. The GO and KEGG pathway analysis of the potential targets of celery seed related to gout disease was performed using the ShinyGO v0.75 app. Molecular docking and molecular dynamics were carried out using Autodock vina and NAMD 2.14 software, respectively. The network analysis identified 16 active compounds and thirteen key targets of celery seed in the treatment of gout. The GO analysis and the KEGG pathway enrichment analysis suggested that the mechanism of action of the chemical constituents of celery seed might be involved in several pathways, notably the PI3K-Akt signaling pathway, Ras signaling pathway, and HIF-1 signaling pathway, respectively. Molecular docking and molecular dynamics revealed that apiumetin might be an important chemical that plays a key role in the pharmacological effect of celery seed. These results might be useful to select the Q-markers to control the quality of the products from celery seeds. Communicated by Ramaswamy H. Sarma


Introduction
Gout (ICD-11 code: FA25, according to the International Classification of Diseases), which is associated with hyperuricemia, is a disease resulting from the deposition of monosodium urate crystals in tissues (Chi et al., 2020;World Health Organization, 2019).According to recent reports, the prevalence of gout in the world ranged from below 1 to 6.8%, while the incidence ranged from 0.58 to 2.89/1000 people per year and tends to increase globally (Dehlin et al., 2020).
The general objectives of gout treatment are to lower uric acid and to reduce inflammation (Cronstein & Terkeltaub, 2006).Although the current drugs used for the treatment of gout, such as allopurinol, colchicine, corticosteroids and NSAIDs are highly effective, the undesirable effects of these drugs when used in the long term are still a big problem to overcome (Cronstein & Terkeltaub, 2006).On the other hand, a large number of herbal medicines were proven to be safe and achieved certain effectiveness in treating gout (Chi et al., 2020).Celery (Apium graveolens L.) has long been considered as a potential herbal medicine for the prevention and treatment of gout because of its anti-hyperuricemic (Soliman et al., 2020), xanthine oxidase (XO) inhibitory (Dolati et al., 2018), anti-inflammatory and analgesic effects (Li et al., 2019).However, the relationship between the chemical constituents and pharmacological activities of this medicinal plant has not been fully investigated yet.
In recent years, network pharmacology has become a new and effective approach to drug research and development with the help of big data and bioinformatics.Network pharmacology describes the complex interrelationships between disease, human biological systems, and drugs through networks.From there, it is possible to understand the mechanism of action as well as to identify the molecular targets related to a specific disease of a drug.Because of its effectiveness and advantages, network pharmacology has quickly been widely applied in the study of active ingredients as well as to explore the mechanism of action of herbal medicines (Shao & Zhang, 2013).Therefore, our study was carried out to discover the potential chemical constituents and molecular targets of celery seed in the treatment of gout using network pharmacology combined with molecular docking and molecular dynamics.

Data preparation
Gout-related targets were collected from the GeneCards database (https://www.genecards.org)(Stelzer et al., 2011) and OMIM database (https://www.omim.org/)(Hamosh et al., 2000) with "gout" as the keyword while the information about the chemical constituents of celery seed was obtained from published literature.Then the structure of each compound was revealed from the PubChem database (https:// pubchem.ncbi.nlm.nih.gov/)(Kim et al., 2016).After that, the ADMETlab web server (https://admetmesh.scbdd.com/service/screening/index)(Xiong et al., 2021) was used to calculate the oral bioavailability of obtained constituents and expressed as a value of F(30%).Compounds with F(30%) � 0.3 would be chosen for further study.
In the next step, the targets of these compounds which might or might not relate to gout were collected using the SwissTargetPrediction web server (http://www.swisstargetprediction.ch) (Gfeller et al., 2013) in which the selected species was Homo sapiens.After the nomenclatures of targets were standardized from the Uniprot database (https://www.uniprot.org)(Consortium, 2019;The UniProt Consortium, 2023) to facilitate analysis and comparison, the potential targets for gout treatment of celery seed were identified as the intersection of the "gout-related targets" dataset and the "targets of compounds presented in celery seed" dataset.

Network building and analysis
The compound-target-pathway network was built and analyzed using Cytoscape 3.9.0software (Shannon et al., 2003).Then the active compounds of celery seed were identified based on the value of the degree (DC), betweenness centrality (BC) and closeness centrality (CC).A protein-protein interaction (PPI) network was also established using the STRING database (https://string-db.org/) in which the organism was set as Homo sapiens and the interaction score must be equal to or greater than 0.900.After that, the DC, BC and CC of each node (protein) were calculated.The larger the degree centrality, betweenness centrality and closeness centrality of a node, the more important that protein was in the protein-protein interaction network (Shao & Zhang, 2013).The cut-off values were defined as the average of DC, BC and CC values, respectively.

Gene ontology (GO) enrichment analysis and the Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis
The GO and KEGG pathway analysis of the potential targets of celery seed related to gout disease was performed using the ShinyGO v0.75 app with the species set as Homo sapiens.
The GO analysis predicted the processes most likely to be relevant to the potential targets, belonging to one of three groups: biological processes, cellular components, and molecular function while KEGG pathway analysis suggested the pathways most likely to be involved in the potential effects of celery in the treatment of gout with a p-value cutoff (FDR) of 0.05.2RAM, gene: RELA) were downloaded from Protein Data Bank (PDB, www.rcsb.org).After that, these proteins were added hydrogen, followed by being assigned the partial charges.Before carrying out the molecular docking procedure, the protocol was validated by redocking the co-crystal ligands and then compared with the experiment.The method was considered accurate and reliable if the RMSD between the co-crystal and re-docked ligands were less than 2 Å (Huang & Zou, 2006).

Molecular docking
In the next steps, the structures of ligands were downloaded from the PubChem database and then converted to pdbqt file.Then the small molecules were docked into the active site of the targets using Autodock vina software (Trott & Olson, 2010).The Discovery Studio Visualizer software was used to analyze the docking results.

Molecular dynamics
In this study, molecular dynamics simulations were performed using NAMD (Nanoscale Molecular Dynamics) version 2.14 (Phillips et al., 2020).First, the proteins were prepared using VMD (Visual Molecular Dynamics) version 1.9.7 while the ligand topology file was obtained from the CHARMM-GUI server (Input Generator module) (Lee et al., 2016).Then, a water cubic box including 9143 TIP3 water molecules (in the case of HSP90a-apiumetin complex) or 9226 TIP3 water molecules (in the case of MAPK14-apiumetin complex) and the corresponding complexes were obtained, respectively.After being neutralized with Na þ and Cl À 0.15M, each system was minimized in 1000 steps followed by being simulated for 40 ns (20 million steps, 2fs/step) in NPT ensemble (Pressure: 1 bar and Temperature: 300 K).The results were recorded after every 0.5 ps.Finally, the number of hydrogen bonds,  the RMSD of the protein backbone and the RMSF of Ca were calculated from the MD simulations using VMD 1.9.7.

Data preparation
Using GeneCards, OMIM database and published literature (Gupta et al., 2019;Hussain et al., 2013;Lin et al., 2007;Qiao et al., 2020), 1551 genes associated with gout disease (1541 from GeneCards, 57 from OMIM with 47 duplicated genes) and 155 chemical constituents of celery seed were collected.Of these, only 68 compounds met the requirements for oral bioavailability, corresponding to 779 targets.The Venn diagram (Figure 1) showed that the intersection of the two datasets was 185.In other words, the components of celery seeds could interact with 185 potential targets in the treatment of gout.

Compound-target-pathway network
The compound-target-pathway network illustrating the effect of celery seed on gout disease consisted of 238 nodes, corresponding to 54 compounds, 185 genes and 65 pathways.The total number of edges in the network is 1364, corresponding to 1364 potential target-compound and pathwaytarget interactions.

Protein-protein interaction network
From 185 potential targets related to gout disease of celery seed, the protein-protein interaction network was obtained and was shown in Figure 2a.
After that, 13 intersectional targets between these nineteen core genes and thirty potential targets identified from the compound-target-pathway network were considered the most promising targets in the gout treatment of celery seed.The DC, BC and CC values of these genes in the PPI network were presented in Table 2.

GO analysis and KEGG pathway analysis
The results of GO analysis were presented in Figure 3, including ten biological processes, ten cellular components, and ten molecular functions most likely to be relevant to 185 potential targets related to the gout disease of celery seed.
In particular, biological processes included: response to oxygen-containing compound, regulation of response to external stimulus, regulation of cell population proliferation, regulation of cell communication, regulation of signaling, regulation of response to stress, cell motility, localization of cell and regulation of locomotion.
The cellular components predicted to participate in the pharmacological effects of celery in the treatment of gout included: receptor complex, membrane raft, membrane microdomain, cell surface, side of membrane, anchoring junction, integral component of plasma membrane, focal adhesion and cell-substrate junction.
The molecular functions revealed by GO analysis were: transmembrane receptor protein tyrosine kinase activity, protein tyrosine kinase activity, protein kinase activity, phosphotransferase activity, kinase activity, transferase activity, adenyl ribonucleotide binding, adenyl nucleotide binding and ATP binding (Figure 3).
The KEGG pathway analysis identified 65 pathways involved in the pharmacological effects of celery in the treatment of gout, of which the 20 first ones were shown in Figure 4.
The obtained results suggested that the pharmacological activity of celery seed in the treatment of gout was based on two main mechanisms: affecting the inflammatory response and the excretion of uric acid in the kidney.Among them, the most prominent was the effect on the inflammatory process.These effects met two goals in the treatment of gout: lowering uric acid and anti-inflammation.The results of the study were also consistent with the reported uric acid-lowering and anti-inflammatory activities of celery seed (Dolati et al., 2018;Li et al., 2019;Soliman et al., 2020).It could also be concluded that the effect of celery on gout was not caused by one substance but was the synergistic effect of many components on different targets.

Molecular docking
In this study, molecular docking was performed to explore the mechanism of binding between 16 potential compounds and core targets.First, the docking protocol was validated by redocking the co-crystal ligand into the active site of the proteins, namely PI3Ka and ERK2.The results were illustrated in Figure 5, which suggested that the docking procedure was accurate and reliable with the RMSD of co-crystal and re-docked ligands of 0.63 and 0.71 Å, in the cases of ERK2 and PI3Ka, respectively.After being validated, molecular docking between 16 potential constituents (from C1 to C16) of celery seed and 12 core proteins was carried out.The results were shown in Figure 6.
As could be seen from Figure 6, AKT1, HSP90a, PPARa, nitric oxide synthase and MAPK14 were the proteins targeted the most while the inhibition of integrin b1 might not be the major mechanism of action of celery seed in the prevention and treatment of gout.Among 16 potential constituents presented in celery seed, apiumetin (C1) and diflufenican (C5) were the compounds that showed the highest activities, followed by cicaprost (C2), ormetoprim (C10) and 3 0 -methoxy-E,E-dienoestrol (C11).According to Long et al. (2022) and Khan and Lee (2022), the binding affinity < À 4.25 kcal/mol suggested that ligands could bind to proteins while the binding energy < À 5 kcal/mol indicated good binding affinity and binding energy < À 7.0 kcal/mol showed strong binding activity.The results from Figure 6 revealed that most of the promising compounds could interact well with the core target with the docking scores ranging from À 4.3 to À 10.1 kcal/mol.The binding modes of apiumetin with HSP90a (DG affinity ¼ À 9.2 kcal/mol) and MAPK14 (DG affinity ¼ À 9.6 kcal/mol), diflufenican with nitric oxide synthase (DG affinity ¼ À 10.1 kcal/mol) and HSP90a (DG affinity ¼ À 9.6 kcal/mol) were illustrated in Figure 7.The results suggested that Met98, Leu107, Val150, Phe138 and Trp162 might be the key residues involving the pharmacological effects of HSP90a.In the case of MAPK14 and nitric oxide synthase, the important residues were identified as Val30, Tyr35, Val38, Ala51, Lys53, Thr106, His107, Leu108, Phe169 and Gln263, Arg266, Pro350, Gly371, Trp372, Arg381, Arg388, respectively.

Molecular dynamics
To explore the stability of the protein-ligand complexes, molecular dynamics were performed.Although Ming Qiao et al. reported that diflufenican (C5) was identified in celery seeds, we supposed that this compound was just a contaminant since it had been widely used as an herbicide.Therefore, in this case, we would only focus on the HSP90aapiumetin and MAPK14-apiumetin complexes.The number of hydrogen bonds, RMSF of Ca and RMSD of these two complexes during 40 ns of simulation were illustrated in Figure 8 and Figure 9, respectively.As could be seen from Figure 8, HSP90a-apiumetin complex was stable with the RMSD values < 2.5 Å during 40 ns of simulations and the average RMSD of HSP90a-apiumetin complex of 1.80 Å.The insignificant difference in the RMSDs of the initial and final states also suggested that apiumetin (C1) could be a good inhibitor of HSP90a.In the case of MAPK14-apiumetin complex, a significant fluctuation was observed with the average RMSD of 2.28 Å.After 30 ns of simulation, the system reached equilibrium.Therefore, the simulation results after 30 ns could be used to analyze the MAPK14-apiumetin complex.
The RMSF of Ca of HSP90a in complex with apiumetin indicated that the complex was stable when there was no large fluctuation observed.On the other hand, all the residue with RMSF > 2.5 Å, including Lys74 and Glu223 were not the key residues involving the pharmacological effects of this protein which suggested that the binding modes of apiumetin and HSP90a revealed by molecular docking were reliable.Similar results were observed when four residues with large fluctuations (RMSF > 3 Å), including Val183, Ala184, Pro351 and Pro352 were not the important residues of MAPK14 that were identified from molecular docking.

Discussion
Since the concept was first introduced by Andrew L. Hopkins, network pharmacology has become an effective approach to identifying the active ingredients and exploring the mechanism of action of herbal medicines which often have complex chemical compositions and unclear modes of action (Shao & Zhang, 2013).This method could provide an overview of the main constituents as well as the complex interaction between the herbal drugs and the targets, thereby helping to design further in vitro and in vivo studies (Chandran et al., 2017).This is one of the advantages of network pharmacology over other virtual screening methods such as ligand-based screening.However, this approach also has some drawbacks that can be mentioned as the data to build the network is collected from databases or published papers.Therefore, the quality of the research depends a lot on the accuracy of the source of information.On the other hand, the results obtained may be incomplete, especially when the number of studies on medicinal plants/diseases is limited (Chandran et al., 2017).
Our study also revealed thirteen prominent targets of constituents of celery seed for treatment of gout, including PIK3R1, PIK3CA, MAPK1, MAPK3, HSP90AA1, AKT1, EGFR, RELA, MAPK14, ITGB1, MAPK8, PPARA and NOS2, respectively.All of them were demonstrated to be involved in inflammatory responses, especially inflammation triggered by urate crystals.In addition, a study conducted in 2003 by Feher et al. (2003) reported that no case developed gout when using PPARa activator, which had also been proved to significantly reduce serum uric acid concentration.The molecular docking and molecular dynamics simulations results suggested that among these core proteins, AKT1, HSP90a, PPARa, nitric oxide synthase and MAPK14 might be the primary targets of celery seed for gout.
The KEGG pathway analysis suggested that the pharmacological activity of celery seed in the treatment of gout involved a variety of biological processes in the body, notably the PI3K-Akt signaling pathway, the Ras signaling pathway, and the HIF-1 signaling pathway.Abnormalities of the PI3K/Akt signaling pathway were implicated in various pathological conditions, including inflammatory processes (Fruman et al., 2017).In addition, inhibition of the PI3K/Akt pathway might reduce uric acid-induced upregulation of ABCG2, thereby increasing renal excretion of uric acid (Chen et al., 2018).Another notable metabolic pathway was the Ras signaling pathway, in which Ras inactivation could decrease gout inflammation and affect IL-6 secretion in vivo (Oh et al., 2021).The HIF-1 signaling pathway also played a key role in inflammatory expression by directly or indirectly participating in the regulation of inflammatory factors such as IL-1b and NLRP3 (Zhang et al., 2013).Other pathways, such as the MAPK pathway, had also been demonstrated to play a key role in AP-1 and NF-jB activation which produced proinflammatory mediators during the acute gout attack (Dalbeth & Haskard, 2005).However, similar to other virtual screening methods, the results given by network pharmacology, molecular docking and MD simulations are only preliminary.Therefore, it is necessary to carry out experimental validation to prove that the obtained results were accurate and reliable.

Figure 1 .
Figure 1.Potential targets related to gout of celery seed.

Figure 2 .
Figure 2. The protein-protein interaction network of a) 185 potential protein targets and b) 19 core targets.

Figure 3 .
Figure 3. GO analysis: biological processes (a), cellular components (b) and molecular function (c) of gout treated by celery seed.

Figure 5 .
Figure 5.The superimposition of co-crystal and re-docked (yellow) ligand of a. ERK2 and b.PI3Ka complex.

Figure 4 .
Figure 4. KEGG pathway analysis of gout treated by celery seed.

Figure 6 .
Figure6.The binding energy of potent compounds of celery seed and core targets of gout disease.The lower the binding energy was, the stronger the interaction between the protein and ligand was.

Figure 8 .
Figure 8.The number of hydrogen bonds (a), RMSF of Ca (b) and RMSD (c) of HSP90a-C1 complexes during 40 ns of simulation.

Figure 9 .
Figure 9.The number of hydrogen bonds (a), RMSF of Ca (b) and RMSD (c) of MAPK14-C1 complexes during 40 ns of simulation.

Table 1 .
Chemical structure of 16 potent compounds in celery seed.

Table 2 .
The BC, CC and DC values of 13 intersectional targets.