Pharmacological mechanisms of carvacrol against hepatocellular carcinoma by network pharmacology and molecular docking

Abstract

Background

Preclinical studies have demonstrated that carvacrol possesses various biological and pharmacological properties, including anti-hepatocellular carcinoma (HCC) effects. However, the molecular basis of its therapeutic action on HCC remains unclear.

Objective

The aim of this study was to investigate and further validate the multi-target therapeutic mechanism of carvacrol against HCC.

Methods

The chemical structure of carvacrol was obtained from the PubChem database, and its potential targets were identified using SwissTargetPrediction, HERB, and BATMAN-TCM. HCC-specific genes were screened from the TCGA-LIHC cohort. The therapeutic targets of carvacrol against HCC were determined through the intersection of these datasets. Subsequently, a multivariate Cox regression prognostic model was established. Molecular docking was performed to analyze the interactions between carvacrol and its therapeutic targets. Additionally, molecular dynamics simulations were conducted to validate the molecular docking results using Discovery Studio 2019 software.

Results

A total of 223 carvacrol targets and 882 HCC-specific genes were identified. Fifteen therapeutic targets of carvacrol against HCC were obtained, including CA2, AR, ALB, AURKA, ALPL, EPHX2, BCHE, IL1RN, AGRN, CRP, DMGDH, APOA1, SOX9, HPX, and CHKA. The prognostic model accurately and independently predicted survival outcomes. AGRN and AURKA were significantly associated with HCC overall survival. Molecular docking and molecular dynamics simulations demonstrated that carvacrol exhibited strong potential for stable binding to the therapeutic targets AGRN and AURKA.

Conclusion

Our findings elucidate the multi-target mechanism of action of carvacrol against HCC, providing a foundation for future research on its application in HCC management.

Keywords

hepatocellular carcinoma carvacrol network pharmacology molecular docking AGRN AURKA

1 Introduction

Liver cancer remains the third leading cause of cancer death globally, with 905,677 new cases as well as 830,180 deaths reported worldwide in 2020.¹ HCC is the most common histological subtype of liver cancer.² HCC cases are typically not detected until intermediate or advanced stages.³ Hence, decision making about systemic treatment remains crucial. Surgery (resection, liver transplantation, etc.) and local-regional treatment (LRT; ablation, etc.) are suitable for the early stages.⁴ additional LRTs (transarterial chemoembolization or radioembolization, etc.) and systemic treatment are required for cases with advanced stage or metastases.⁵ Sorafenib, a tyrosine kinase inhibitor, has been approved by the FDA (Food and Drug Administration) as the only systemic option against inoperable HCC. However, novel systemic therapies, particularly immunotherapy, have recently shown promise against HCC.⁶ Phase III clinical trials have not consistently demonstrated survival benefit of immune checkpoint inhibitors in unresectable HCC under the first- and second-line clinical settings.⁷ Moreover, the objective response rate remains only 15–20% when using these drugs as monotherapy for HCC.⁸ Several Phase III trials (ORIENT-32, COSMIC-312, HIMALAYA, etc.) have combined anti-programmed death 1 or anti-programmed death ligand 1 drugs in combination with anti-cytotoxic T-lymphocyte associated protein 4 or anti-angiogenic drugs, and showed more favorable clinical outcomes in comparison to sorafenib monotherapy against unresectable HCC.^9–11 Therefore, novel therapeutic options are essential for treating HCC.

Agents with a single target are often ineffective in alleviating HCC, while traditional Chinese medicine demonstrates unique advantages in HCC therapy due to its multi-target activity and pleiotropic effects. Carvacrol is a monoterpene phenol (5-isopropyl-2-methylphenol) extracted from various aromatic plants, including thyme and oregano,¹² it's currently used in low concentrations as a food flavoring agent, preservative, and fragrance component in cosmetics.¹³ Preclinical studies have shown that carvacrol possesses diverse biological and pharmacological properties, including anti-HCC activity. For instance, carvacrol reduces the viability of HCC cells by downregulating SLC6A3.¹⁴ Additionally, this compound exerts anti-proliferative and pro-apoptotic effects on HCC cells.¹⁵ Carvacrol reduces lipid peroxidation, mitigates hepatic cell injury, and protects the antioxidant system in diethylnitrosamine-induced HCC.¹⁶ It also demonstrates potential anticancer activity by inhibiting cellular proliferation and preventing metastasis in diethylnitrosamine-induced HCC.¹⁷ Carvacrol protects against paracetamol-induced toxicity in HCC cells.¹⁸ However, to date, no systematic study has been conducted to identify the therapeutic targets of carvacrol against HCC.

Pharmacotranscriptomics is a powerful approach for assessing the therapeutic efficacy of drugs and identifying new drug targets.¹⁹ By utilizing network pharmacology and molecular docking approaches, this study elucidated the mechanism of action of carvacrol against HCC, thereby laying a foundation for further research and the potential clinical application of carvacrol in HCC therapy.

2 Materials and methods

2.1 Prediction of carvacrol targets

The PubChem database (https://pubchem.ncbi.nlm.nih.gov) is a widely used chemical information resource that supports drug discovery and chemical biology research.²⁰ The 2D chemical structure of carvacrol was obtained from the PubChem database. SwissTargetPrediction (http://www.swisstargetprediction.ch) is a web tool that predicts the most probable protein targets of small molecule compounds based on similarity principles via reverse screening.²¹ The HERB (http://herb.ac.cn/), a high-throughput, experiment- and reference-guided web server for traditional Chinese medicine, provides 12,933 targets and 28,212 diseases to 7263 herbs and 49,258 compounds.²² The targets of carvacrol were identified using the SwissTargetPrediction and HERB databases. BATMAN-TCM (http://bionet.ncpsb.org/batman-tcm) is an online bioinformatics analysis tool for evaluating the molecular mechanisms of traditional Chinese medicine,²³ It includes predicting the targets of traditional Chinese medicine compounds; functional enrichment analysis of targets, establishment of compound-target-pathway/disease relevant networks as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched by targets, and comparison of different traditional Chinese medicine.²³ Utilizing the BATMAN-TCM tool, carvacrol targets were predicted in accordance with the a score cutoff of 5 and an adjusted p-value cutoff of 0.05. The Universal Protein Resource (UniProt; http://sparql.uniprot.org/) is an integrative resource for protein sequence and annotation data, covering over 120 million proteins.²⁴ The gene name of targets was matched and corrected with the UniProt. The targets predicted by SwissTargetPrediction, HERB, and BATMAN-TCM were merged and deduplicated.

2.2 Data acquisition of TCGA-LIHC cohort

RNA sequencing (RNA-seq) profiles of HCC patients were obtained from The Cancer Genome Atlas (TCGA) dataset (https://portal.gdc.cancer.gov/), including 374 HCC and 50 normal liver tissue specimens. Additionally, clinicopathological data from 377 HCC cases were downloaded from the TCGA-LIHC cohort.

2.3 Analysis of HCC-specific genes

Transcriptional levels of genes were compared in 374 HCC versus 50 normal liver tissue specimens using the limma package.²⁵ Based on the screening criteria of an adjusted p-value < 0.05 and |log2 fold change (FC)| ≥ 1.0, genes with significant differential expression were identified and considered HCC-specific genes. A volcano plot and heatmap of the HCC-specific genes were generated.

2.4 Screening shared targets

Carvacrol targets and HCC-specific genes were imported into the Venny online tool (version 2.1; https://bioinfogp.cnb.csic.es/tools/venny/index.html). The overlapping targets were identified as therapeutic targets of carvacrol against HCC. Using the ggpubr package, the relative expression values of these targets were calculated. Based on the median expression value of each target, HCC cases were divided into two groups. Kaplan-Meier survival curves were plotted, and p-values were calculated using the log-rank test.

2.5 Analysis of protein-protein interactions

The STRING resource (https://string-db.org/) was used to evaluate and integrate protein-protein interactions, including both direct (physical) and indirect (functional) interactions.²⁶ Therapeutic targets of carvacrol against HCC were imported into the STRING online tool. The criteria were set to hide unconnected targets and to apply the highest confidence threshold of 0.900. The resulting protein-protein interaction network was obtained. Hub genes were identified using the cytoHubba plugin in Cytoscape software (version 3.7.2).²⁷

2.6 Functional enrichment analysis

ClusterProfiler package offers the process of biological-term classification as well as enrichment analysis of of specific gene sets.²⁸ Using the ClusterProfiler package, Gene Ontology (GO) and KEGG pathway enrichment analyses of the therapeutic targets of carvacrol against HCC were performed.

2.7 Prognostic model construction

A prognostic model of the therapeutic targets of carvacrol against HCC was constructed using multivariate Cox regression analysis, a prognostic model of therapeutic targets of carvacrol against HCC was constructed. The risk score for each HCC case was calculated using the following formula:

risk score = \sum_{n = 1}^{n} c o e f f i c i e n t (n) * e x p r e s s i o n v a l u e (n),

where coefficient (n) represented the regression coefficient, while expression value (n) indicated the relative expression value of each therapeutic target of carvacrol against HCC in the prognostic model. The 377 HCC cases with clinicopathological information were evenly divided into training and validation sets. With the median risk score, HCC cases in training, verification and entire sets were categorized into high- and low-risk groups, respectively. The survival status (alive or deceased) was visualized for both groups. Kaplan-Meier survival curves were plotted for each group, and p-values were calculated using the log-rank test with the “survminer” and “survival” packages. The predictive specificity and sensitivity of the prognostic model were assessed using Receiver Operating Characteristic (ROC) curves, and the area under the curve (AUC) was calculated using the “survival ROC” package. Univariate Cox regression analysis was performed to investigate the associations of the prognostic model and clinicopathological factors with overall survival. Risk factors were defined as those with a hazard ratio >1 and p-value < 0.05, while protective factors were defined as those with a hazard ratio <1 and p-value < 0.05. Additionally, independent prognostic factors were identified using multivariate Cox regression analysis.

2.8 Gene Set Enrichment Analysis (GSEA)

GSEA, a computational approach, can determine whether a specific gene set show a difference between two biological phenotypes.²⁹ Based on the enrichment score, GO and KEGG pathways associated with each therapeutic target of carvacrol against HCC were evaluated. The reference gene sets “c5.go.v7.5.1.symbols” and “c2.cp.kegg.v7.5.1.symbols” were downloaded from the Molecular Signatures Database (MSigDB).³⁰

2.9 Molecular docking

The structural formula of carvacrol was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The corresponding 3D structures were generated using Chem3D software and exported in mol2 format. The Protein Data Bank (PDB; http://www.rcsb.org/) is the global archive for biological macromolecular structure data.³¹ The domains of each therapeutic target of carvacrol against HCC were obtained from the PDB and saved in pdb format. PyMOL software was used to remove water and phosphate groups from the protein. AutoDockTools was used to convert the mol2 format of carvacrol and the pdb format of the domains to pdbqt format, followed by the identification of the active pocket. The Vina script was then run to calculate molecular binding energy and visualize the molecular docking. Additionally, Discovery Studio 2019 was used to identify docking sites and compute the LibDockScore for flexible binding. The molecular docking results were imported into PyMOL for visualization of molecular docking conformations. If the binding energy was less than 0, spontaneous binding between the ligand and receptor could occur. A Vina binding energy ≤ −5.0 kcal·mol⁻¹ and LibDockScore >50 indicated stable docking. The 3D and 2D molecular docking results of ligand-receptor complexes were visualized to assess predictive reliability.

2.10 Molecular dynamics simulation

Molecular dynamics simulation was performed to verify the binding of component-protein complexes with the best molecular docking results, using the Simulation and Standard Dynamics Cascade modules of Discovery Studio 2019 software. During the simulation, the molecular parameters of the ligands were based on the Charmm36 force field, while the receptors used Charmm positioning. The protein-ligand complexes were solvated using the Solvation module. The molecular dynamics simulation was then executed. Using the Analyze Trajectory module, structural properties of the molecular dynamics trajectories, such as geometric properties, the number of non-bonding interactions per simulation frame, root mean square deviation (RMSD), and root mean square fluctuation (RMSF), were analyzed. Additionally, non-bonded interactions between peptides and proteins were detected.

3 Results

3.1 The chemical structure of carvacrol and its potential targets

We first retrieved the 2D chemical structure of carvacrol from the PubChem database, as shown in Figure 1A. To identify the targets of carvacrol, we used the SwissTargetPrediction, HERB, and BATMAN-TCM tools. As a result, 106, 13, and 122 carvacrol targets were predicted by SwissTargetPrediction, HERB, and BATMAN-TCM, respectively. After merging and removing duplicates, 223 carvacrol targets were identified. Figure 1B shows the distribution of carvacrol targets: 51% were predicted by BATMAN-TCM, 44% by SwissTargetPrediction, and 5% by HERB. We also analyzed the target classes for the top 25 carvacrol targets predicted by the SwissTargetPrediction web server. These included kinase (20%), enzyme (16%), family A G protein-coupled receptor (16%), oxidoreductase (8%), nuclear receptor (8%), lyase (8%), eraser (8%), voltage-gated ion channel (4%), cytochrome P450 (4%), and secreted protein (4%), as shown in Figure 1C.

Figure 1.

Identification of therapeutic targets of carvacrol against HCC. (A) The 2D chemical structure of carvacrol from the SwissTargetPrediction web server. (B) The distribution of carvacrol targets across the SwissTargetPrediction, HERB as well as BATMAN-TCM tools. (C) The target classes for the first 25 carvacrol targets predicted by the SwissTargetPrediction web server. (D) Volcano diagram for the genes with differential expression between 374 HCC and 50 normal liver tissue specimens in the TCGA-LIHC cohort. Red represents up-regulation; green represents down-regulation; and black represents no difference in gene expression. (E) Heatmap illustrating the genes with significant differential expression on the basis of adjusted p-value. Red represents up-regulation; and blue represents down-regulation. (F) Venn diagram for shared genes between carvacrol targets and HCC-specific genes.

3.2 Selection of HCC-specific genes

Based on the criteria of an adjusted p-value < 0.05 and |log2 FC| ≥ 1.0, 882 differentially expressed genes between 374 HCC and 50 normal liver tissue specimens from the TCGA-LIHC cohort were identified as HCC-specific genes, including 487 upregulated and 395 downregulated genes (Figure 1D). The heatmap illustrates the expression patterns of all HCC-specific genes across the 374 HCC and 50 normal liver tissue samples (Figure 1E).

3.3 Identification of therapeutic targets of carvacrol against HCC

After intersecting carvacrol targets with HCC-specific genes, fifteen therapeutic targets of carvacrol against HCC were identified: CA2, AR, ALB, AURKA, ALPL, EPHX2, BCHE, IL1RN, AGRN, CRP, DMGDH, APOA1, SOX9, HPX, and CHKA (Figure 1F).

3.4 Interactions between protein products of therapeutic targets of carvacrol against HCC

Using the STRING online tool, interactions between the protein products of the therapeutic targets of carvacrol against HCC were assessed, as shown in Figure 2A. Hub genes (CA2, AR, ALB, AURKA, ALPL, EPHX2, BCHE, IL1RN, AGRN, CRP, DMGDH, APOA1, SOX9, HPX, and CHKA) were selected based on the cytoHubba plugin of Cytoscape software (Figure 2B). Table 1 lists the degree of the therapeutic targets of carvacrol against HCC.

Figure 2.

Interactions between protein products of therapeutic targets of carvacrol against HCC and their biological functions and pathways. (A) Protein-protein interaction network of protein products of therapeutic targets of carvacrol against HCC. (B) Hub genes among protein products of therapeutic targets of carvacrol against HCC. (C, D) Biological processes, (E, F) Cellular components, (G, H) molecular functions, and (I, J) KEGG pathways enriched by protein products of therapeutic targets of carvacrol against HCC.

Table 1.

The degree of therapeutic targets of carvacrol against HCC.

Gene Symbol	logFC	Adjusted p-value	Degree
CA2	−1.20645	1.16E-11	6
AR	−1.0512	4.82E-09	6
ALB	−1.72895	1.86E-25	13
AURKA	2.678679	4.13E-29	3
ALPL	−1.44342	1.05E-11	4
EPHX2	−1.42558	1.6E-23	2
BCHE	−1.75093	9.54E-16	7
IL1RN	−2.03246	2.19E-36	4
AGRN	1.052369	1.14E-08	5
CRP	−1.06507	0.001827	8
DMGDH	−1.51954	3.14E-20	2
APOA1	−1.24839	2.75E-09	9
SOX9	1.367765	4.41E-07	6
HPX	−1.66305	8.01E-21	4
CHKA	1.230007	8.24E-13	5

3.5 Biological functions and pathways of therapeutic targets of carvacrol against HCC

Using the clusterProfiler package, the biological functions and pathways of the therapeutic targets of carvacrol against HCC were explored. Biological processes such as ammonium ion metabolic process, cell-cell adhesion involved in gastrulation, epithelial cell proliferation in prostate gland development, male sex determination, negative regulation of heterotypic cell-cell adhesion, positive regulation of reproductive processes, prostate gland growth, regulation of cell-cell adhesion in gastrulation, response to estrogen, and response to steroid hormones were significantly associated with the therapeutic targets of carvacrol against HCC (Figure 2C, D). Additionally, these targets were involved in various cellular components, including blood microparticle, chylomicron, collagen-containing extracellular matrix, endocytic vesicle lumen, endoplasmic reticulum lumen, low-density lipoprotein particle, meiotic spindle, nuclear envelope lumen, pronucleus, and spindle pole centrosome, as shown in Figure 2E, F. Additionally, they exhibited molecular functions such as amyloid-beta binding, beta-catenin binding, carboxylic ester hydrolase activity, drug binding, dystroglycan binding, lipoprotein particle binding, lipoprotein particle receptor binding, protein-lipid complex binding, steroid binding, and toxic substance binding, etc. (Figure 2G, H). KEGG pathways such as African trypanosomiasis, chemical carcinogenesis-receptor activation, collecting duct acid secretion, folate biosynthesis, glycine, serine and threonine metabolism, nitrogen metabolism, oocyte meiosis, proximal tubule bicarbonate reclamation, thiamine metabolism, and vitamin digestion and absorption were significantly enriched (Figure 2I, J).

3.6 Construction of a prognostic model on the basis of therapeutic targets of carvacrol against HCC

A multivariate Cox regression model was developed for HCC based on the therapeutic targets of carvacrol against HCC. Using the regression coefficient and expression value of each therapeutic target, the risk score for 377 HCC cases in the TCGA-LIHC cohort was calculated. The cohort was divided into training and validation sets in a 1:1 ratio. Using the median risk score as the cutoff, cases in all three subsets were classified into high- and low-risk groups. Survival status in the training, validation, and entire sets was assessed. In all three sets, a higher number of deceased cases were observed in the high-risk group compared to the low-risk group (Figure 3A-C). The survival differences between the subsets were then compared. Low-risk cases showed a significant advantage in overall survival compared to high-risk cases across all three sets (Figure 3D-F). To assess the predictive specificity and sensitivity of the prognostic model, ROC curves were plotted. The AUC values were greater than 0.7 across all three sets (Figure 3G-1). This confirmed that the prognostic model demonstrated excellent specificity and sensitivity in predicting overall survival. Additionally, univariate Cox regression analysis was performed to evaluate the associations between the prognostic model and clinicopathological factors with HCC prognosis. In Figure 3J, the prognostic model, pathological stage, T stage, and M stage were associated with a more unfavorable prognosis. Further multivariate Cox regression analysis was performed to identify independent prognostic factors. As shown in Figure 3K, the prognostic model was identified as an independent risk factor for HCC prognosis.

Figure 3.

Construction of a prognostic model on the basis of therapeutic targets of carvacrol against HCC utilizing multivariate cox regression analysis. (A-C) The distribution of survival status (green, alive; red, dead) in high- and low-risk subsets in training, verification, and entire sets from the TCGA-LIHC cohort. (D-F) Kaplan-Meier curves of overall survival between high- and low-risk subsets in training, verification, and entire sets from the TCGA-LIHC cohort. P values were computed with log-rank test. (G-I) ROC curves in training, verification, and entire sets from the TCGA-LIHC cohort. (J) Univariate cox regression analysis of the associations of the prognostic model and clinicopathological factors with HCC prognosis in the TCGA-LIHC cohort. Hazard ratio, confidence interval as well as p-value were computed. (K) Multivariate cox regression analysis for selecting the independent prognostic factors in the TCGA-LIHC cohort.

3.7 Clinical implications of therapeutic targets of carvacrol against HCC: AGRN and AURKA

The clinical implications of each therapeutic target of carvacrol against HCC were further investigated. As a result, both AGRN and AURKA were significantly associated with poor overall survival in HCC cases, while the other therapeutic targets did not show significant correlations with HCC prognosis (Figure 4A, B). Next, we examined their interactions with pathological stage. No significant differences in AGRN expression were observed across different pathological stages (Figure 4C). However, AURKA expression was higher in stage III cases compared to stage I (Figure 4D). This indicated that AURKA might contribute to disease progression.

Figure 4.

Clinical implications and relevant molecular mechanisms of therapeutic targets of carvacrol against HCC: AGRN and AURKA. (A) Kaplan-Meier curves of overall survival between high and low AGRN expression in the TCGA-LIHC cohort. P-value was computed with log-rank test. (B) Kaplan-Meier curves of overall survival between high and low AURKA expression in the TCGA-LIHC cohort. (C, D) Distribution of AGRN and AURKA expression values among distinct pathological stages in the TCGA-LIHC cohort. (E, F) GO enrichment results of AGRN and AURKA. (G, H) KEGG pathways linked to AGRN and AURKA.

3.8 Molecular mechanisms underlying therapeutic targets of carvacrol against HCC: AGRN and AURKA

Further analysis was conducted to investigate the molecular mechanisms underlying the therapeutic targets of carvacrol against HCC, specifically AGRN and AURKA. In GO enrichment analysis, AGRN was positively associated with processes such as the negative regulation of cellular amide metabolic process, intermediate filament, intermediate filament cytoskeleton, mRNA 3'UTR binding, and RNA binding involved in post-transcriptional gene silencing (Figure 4E). Meanwhile, AURKA showed positive correlations with processes such as RNA-mediated gene silencing by inhibition of translation, sensory perception of smell, snRNA metabolic process, snRNA processing, and olfactory receptor activity (Figure 4F). In KEGG pathway enrichment analysis, AGRN was negatively linked to pathways such as aminoacyl-tRNA biosynthesis, fatty acid metabolism, glycine, serine and threonine metabolism, neuroactive ligand-receptor interaction, and primary bile acid biosynthesis (Figure 4G). Moreover, AURKA was negatively associated with fatty acid metabolism, glycine, serine and threonine metabolism, neuroactive ligand-receptor interaction, olfactory transduction, and primary bile acid biosynthesis (Figure 4H).

3.9 Molecular docking of carvacrol and therapeutic targets of carvacrol against HCC

To verify the interaction between carvacrol and its therapeutic targets against HCC, molecular docking was performed. Using Chem3D software, the 3D structure of carvacrol was generated in mol*2 format. Additionally, we downloaded the 3D structures of the therapeutic targets of carvacrol against HCC from the PDB database in pdb format. Using AutoDockTools software, we converted the formats of carvacrol and its therapeutic targets into pdbqt format to identify the active pocket. After running the Vina script, the binding energies of the ligands and receptors were calculated. As shown in Table 2, the binding energies of most docking complexes formed between the therapeutic targets of carvacrol against HCC and carvacrol were lower than −5.0 kcal·mol⁻¹, indicating stable docking. Additionally, we used Discovery Studio 2019 software to dock carvacrol with its targets and calculate the LibDockScore. In Table 2, carvacrol can perform semi-flexible docking with the receptor-ligand of the targets. The LibDockScores of the docking models formed by the targets and carvacrol were all greater than 50. Among all the targets, AGRN formed the most stable complex with carvacrol, followed by AURKA. Due to the formation of different hydrogen bonds and hydrophobic interactions at various amino acid sites, different docking sites were observed between carvacrol and the respective targets. Finally, we performed 3D (Figure 5A-O) and 2D (Figure 6A-O) molecular docking between carvacrol and therapeutic targets of carvacrol against HCC.

Figure 5.

3D model of molecular docking between carvacrol and therapeutic targets of carvacrol against HCC. (A) CA2; (B) AR; (C) ALB; (D) AURKA; (E) ALPL; (F) EPHX2; (G) BCHE; (H) IL1RN; (I) AGRN; (J) CRP; (K) DMGDH; (L) APOA1; (M) SOX9; (N) HPX; and (O) CHKA.

Figure 6.

2D model of molecular docking between carvacrol and therapeutic targets of carvacrol against HCC. (A) CA2; (B) AR; (C) ALB; (D) AURKA; (E) ALPL; (F) EPHX2; (G) BCHE; (H) IL1RN; (I) AGRN; (J) CRP; (K) DMGDH; (L) APOA1; (M) SOX9; (N) HPX; and (O) CHKA.

Table 2.

Molecular docking of carvacrol and therapeutic targets of carvacrol against HCC.

Domain	Vina (kcal·mol⁻¹)	RMSD	DS (LibDockScore)	Hydrogen bond interaction	Hydrophobic interaction
CA2(3F6U)	−5.5	1.299	82.5727	-	SER:214,VAL:213,TYR:228,
					CYS:191,ALA:190,TRP:215
AR(2F2 K)	−7.4	1.624	84.9137	-	LYS:262,PRO:215,PRO:211,
					LEU:212,SER:210,ILE:260
ALB(4LB2)	−5.9	2.053	63.1487	-	PRO:180.LEU:179,ALA:176
AURKA(5DOS)	−6.8	1.412	70.5223	-	TYR:295,LYS:258,LEU:296,
					TRP:277,CYS:290
ALPL(1ZEF)	−5.2	17.864	67.3863	-	ALA:294,LEU:253,ARG:250,
					PRO:286,MET:254
EPHX2(3ANS)	−6.4	1.740	66.7548	LEU:408,TRP:525	PHE:497,LEU:417
BCHE(4AQD)	−6.4	1.650	54.8526	ILE:69	LEU:274,PRO:491
IL1RN(1IRA)	−3.3	7.646	66.067	ASP:128	LEU:12,ARG:22,ILE:89,PRO:25,
					PRO:23,PHE:127,TYR:124
AGRN(3V64)	−5.0	0.727	61.9421	HIS:1861	GLU:1866,SER:1868
CRP(2ZMF)	−7.3	1.368	72.3969	ASP:305,CYS:287	VAL:356,TYRR:362,PHE:304,
					ALA:288,ARG:286,VAL:385
DMGDH(5L46)	−7.1	0.853	58.9898	TRP:727,SER:599	PHE:537,ALA:726,TYR:676
APOA1(4UYG)	−6.4	1.064	61.8957	-	HIS:387,HIS:352,MET:453,
					PRO:454
SOX9(4EUW)	−4.8	1.767	75.1767	-	DA:8,ARG:162,PRO:108,
					ALA:111,VAL:155,ALA:158
HPX(1SU3)	8.7	0.500	78.9208	-	HIS:218,VAL:215,LEU:181,
					ALA:182,CYS:92
CHKA(4BR3)	−7.5	1.867	74.5174	ASN:345,	TYR:440,TRP:420,TRP:423,
					TYR:333

3.10 Molecular dynamics simulation analysis of carvacrol and therapeutic targets of carvacrol against HCC (AGRN, and AURKA)

The optimal conformations produced by the docking of carvacrol with its therapeutic targets against HCC (AGRN and AURKA) were verified through molecular dynamics simulation analysis. To simulate the physiological environment, 4928 water molecules, 17 sodium ions, and 13 chloride ions were added to the AGRN-carvacrol ligand-receptor complex, while 5700 water molecules, 15 sodium ions, and 20 chloride ions were added to the AURKA-carvacrol complex. To assess the structural stability of the protein-ligand complexes during molecular dynamics simulation, RMSD values were calculated for the complexes over a 100 ns simulation period. In Figure 7A, Both complexes remained stable after the 100 ns molecular dynamics simulation. The RMSD values of the AGRN-carvacrol complex fluctuated primarily between 1.35592 and 1.7014, with a mean RMSD value of 1.46382, while the RMSD values of the AURKA-carvacrol complex fluctuated between 1.22534 and 1.68954, with a mean RMSD value of 1.49804. The RMSD fluctuations of both complexes were within a reasonable range, indicating that the structures reached equilibrium after the simulation, confirming that the AGRN-carvacrol and AURKA-carvacrol complexes remained stable throughout the molecular dynamics simulation process. The mean RMSD values of both complexes were less than 2.0, with AGRN-carvacrol showing greater stability in binding compared to the AURKA-carvacrol complex. To further analyze the flexibility of individual amino acids in the complexes during the simulation, RMSF values of all amino acids were calculated. In Figure 7B, C, In the AGRN-carvacrol complex, large fluctuations were observed around the following amino acids: Ile1777, Glu1778, Ser1779, Glu1780, Ala1816, Glu1817, Arg1818, Gln1900, Lys1901, Leu1902, Pro1903, Thr1947, and Pro1948. Similarly, the AURKA-carvacrol complex showed large fluctuations around Ser283, Ser284, Cys290, Pro349, Asp350, Phe351, Val352, Ser388, Lys389, and Pro390. Despite these fluctuations, the overall structures tended to be stable. The RMSF values of amino acid residues in the AGRN-carvacrol complex were smaller than those in the AURKA-carvacrol complex, indicating that the reduced fluctuation in the AGRN-carvacrol complex contributed to its greater stability. The heatmaps of hydrogen bonds during molecular docking are shown in Figure 7D, E. Additionally, hydrogen bond interactions were observed in all conformations (especially the red bars), indicating that these hydrogen bonds were persistent and stable. The simulation results of the AGRN-carvacrol and AURKA-carvacrol protein-ligand complexes are illustrated in Figure 8A-F and Figure 9A-F.

Figure 7.

Molecular dynamics simulation analysis of carvacrol and therapeutic targets of carvacrol against HCC (AGRN, and AURKA). (A) RMSD value change during molecular dynamics simulation process (green: AGRN-carvacrol protein ligand complex; red: AURKA-carvacrol protein ligand complex). (B) Change of RMSF values of AGRN-carvacrol protein ligand complex. (C) Change of RMSF values of AURKA-carvacrol protein ligand complex. (D) Heat map of hydrogen bond of AGRN-carvacrol protein ligand complex. (E) Heat map of hydrogen bond of AURKA-carvacrol protein ligand complex.

Figure 8.

Molecular dynamics simulation results of AGRN-carvacrol protein ligand complex. (A) 0 ns; (B) 20 ns; (C) 40 ns; (D) 60 ns; (E) 80 ns; (F) 100 ns.

Figure 9.

Molecular dynamics simulation results of AURKA-carvacrol protein ligand complex. (A) 0 ns; (B) 20 ns; (C) 40 ns; (D) 60 ns; (E) 80 ns; (F) 100 ns.

4 Discussion

Research in traditional Chinese medicine has increasingly focused on high-throughput transcriptomic screening to explore the molecular effects of compounds.³² Natural compounds can provide HCC patients with improved survival outcomes, while offering lower systemic toxicity and fewer side effects.³³ Limited preclinical studies have suggested the anti-HCC properties of carvacrol. Therefore, research into the underlying mechanisms of carvacrol's therapeutic effects may provide deeper insights into HCC and offer a research model for its application in HCC therapy. Using network pharmacology and molecular docking approaches, this study identified fifteen therapeutic targets of carvacrol against HCC, including CA2, AR, ALB, AURKA, ALPL, EPHX2, BCHE, IL1RN, AGRN, CRP, DMGDH, APOA1, SOX9, HPX, and CHKA, highlighting the multi-target therapeutic potential of carvacrol in HCC treatment.

A prognostic model based on the fifteen therapeutic targets of carvacrol against HCC was developed, which accurately and independently predicted survival outcomes. Previous research has highlighted the functional roles of these targets in HCC. For example, CA2 (carbonic anhydrase 2) plays a critical role in regulating ion transport and pH balance. CA2 has been shown to inhibit epithelial-mesenchymal transition and metastasis in HCC and is associated with favorable overall survival outcomes.³⁴ Additionally, CA2 has the potential to predict HCC recurrence and is correlated with clinicopathological factors, including α-fetoprotein levels, microvascular invasion, and TNM stage.³⁵ Quantitative secretome analysis also reveals that serum CA2 independently predicts recurrence and overall survival in HCC patients.³⁶ Given the pronounced sexual dimorphism of HCC, sex hormone receptors have been implicated in the pathophysiological processes of the disease.³⁷ AR (androgen receptor) has dual and opposing roles in VETC-dependent and invasion-dependent metastasis of HCC.³⁸ ALB (albumin) is associated with aggressive metastasis in HCC cases, and it plays a role in reducing HCC invasion and metastasis.³⁹ Downregulation of EPHX2 (epoxide hydrolase 2) is associated with HCC progression and independently predicts overall survival in patients.⁴⁰ BCHE (Butyrylcholinesterase) can accurately predict a sustained complete response following transarterial chemoembolization and serves as a prognostic marker for HCC.⁴¹ IL1RN (interleukin 1 receptor antagonist) is expressed at low levels in HCC and is associated with tumor immunity.⁴² Consistent with previous research,⁴³ AGRN (Agrin) is associated with poorer overall survival in HCC. AURKA (aurora kinase A), a serine/threonine kinase, plays a crucial role in chromosomal separation and mitotic spindle stability by interacting with the centrosome during mitosis. The current study suggests that AURKA upregulation is linked to poor overall survival and disease progression. Depletion of AURKA, in combination with HSF1 suppression, reduces proliferation and induces apoptosis in HCC by activating endoplasmic reticulum stress.⁴⁴ In inflammatory hepatocellular adenoma, serum CRP (C-reactive protein) is often diffusely expressed in tumoral hepatocytes, with a clear demarcation from the surrounding non-tumor liver tissue.⁴⁵ DMGDH (Dimethylglycine dehydrogenase) attenuates HCC metastases.⁴⁶ Downregulation of APOA1 (apolipoprotein A1) is associated with DNA methylation and poor overall survival in HCC.⁴⁷ SOX9 (SRY-box transcription factor 9) maintains cancer stem cell features as well as accelerates metastases in HCC.^48,49 Serum HPX (hemopexin) is regarded as an underlying marker of HCC.⁵⁰ CHKA (choline kinase alpha) modulates the interaction between EGFR and the mechanistic target of rapamycin complex 2 in HCC, thereby accelerating resistance and tumor progression.⁵¹

Our molecular docking studies demonstrated the interaction between carvacrol and its therapeutic targets against HCC. Molecular dynamics simulation analysis confirmed the stable binding between carvacrol and two of these targets (AGRN and AURKA), indicating direct pharmacodynamic associations of carvacrol with AGRN and AURKA for the treatment of HCC.

5 Conclusion

Collectively, by integrating network pharmacology and molecular docking approaches, this study examined and further validated the multi-target therapeutic mechanism of carvacrol against HCC, particularly focusing on AGRN and AURKA. These findings may provide a foundation for further research and the clinical application of carvacrol in HCC therapy.

Footnotes

Abbreviations

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability

The data that support the findings of this study are available in PubChem database (https://pubchem.ncbi.nlm.nih.gov), SwissTargetPrediction (http://www.swisstargetprediction.ch), HERB (http://herb.ac.cn/), BATMAN-TCM (http://bionet.ncpsb.org/batman-tcm) and The Universal Protein Resource ().

References

Sung

Ferlay

Siegel

, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2021; 71: 209–249.

Liu

Niu

Qiu

. A five-gene signature based on stromal/immune scores in the tumor microenvironment and its clinical implications for liver cancer. DNA Cell Biol 2020; 39: 1621–1638.

Tan

DJH

Lin

, et al. Clinical characteristics, surveillance, treatment allocation, and outcomes of non-alcoholic fatty liver disease-related hepatocellular carcinoma: a systematic review and meta-analysis. Lancet Oncol 2022; 23: 521–530.

Marrero

Kulik

Sirlin

, et al. Diagnosis, staging, and management of hepatocellular carcinoma: 2018 practice guidance by the American association for the study of liver diseases. Hepatology 2018; 68: 723–750.

Altayar

O'Shea

, et al. AGA Clinical practice guideline on systemic therapy for hepatocellular carcinoma. Gastroenterology 2022; 162: 920–934.

Sonbol

Riaz

Naqvi

SAA

, et al. Systemic therapy and sequencing options in advanced hepatocellular carcinoma: a systematic review and network meta-analysis. JAMA Oncol 2020; 6: e204930.

Finn

Ryoo

Merle

, et al. Pembrolizumab as second-line therapy in patients with advanced hepatocellular carcinoma in KEYNOTE-240: a randomized, double-blind, phase III trial. J Clin Oncol 2020; 38: 193–202.

Ruf

Heinrich

Greten

. Immunobiology and immunotherapy of HCC: spotlight on innate and innate-like immune cells. Cell Mol Immunol 2021; 18: 112–127.

Kudo

. Combination cancer immunotherapy with molecular targeted agents/anti-CTLA-4 antibody for hepatocellular carcinoma. Liver Cancer 2019; 8: 1–11.

10.

Finn

Qin

Ikeda

, et al. Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. N Engl J Med 2020; 382: 1894–1905.

11.

Greten

Lai

, et al. Targeted and immune-based therapies for hepatocellular carcinoma. Gastroenterology 2019; 156: 510–524.

12.

Marchese

Arciola

Coppo

, et al. The natural plant compound carvacrol as an antimicrobial and anti-biofilm agent: mechanisms, synergies and bio-inspired anti-infective materials. Biofouling 2018; 34: 630–656.

13.

Suntres

Coccimiglio

Alipour

. The bioactivity and toxicological actions of carvacrol. Crit Rev Food Sci Nutr 2015; 55: 304–318.

14.

Yin

Chen

, et al. Screening and validation of a carvacrol-targeting viability-regulating protein, SLC6A3, in liver hepatocellular carcinoma. Dis Markers 2022; 2022: 3736104.

15.

Yin

Yan

, et al. Anti-proliferative and pro-apoptotic effect of carvacrol on human hepatocellular carcinoma cell line HepG-2. Cytotechnology 2012; 64: 43–51.

16.

Jayakumar

Madankumar

Asokkumar

, et al. Potential preventive effect of carvacrol against diethylnitrosamine-induced hepatocellular carcinoma in rats. Mol Cell Biochem 2012; 360: 51–60.

17.

Subramaniyan

Krishnan

Balan

, et al. Carvacrol modulates instability of xenobiotic metabolizing enzymes and downregulates the expressions of PCNA, MMP-2, and MMP-9 during diethylnitrosamine-induced hepatocarcinogenesis in rats. Mol Cell Biochem 2014; 395: 65–76.

18.

Palabiyik

Karakus

Halici

, et al. The protective effects of carvacrol and thymol against paracetamol-induced toxicity on human hepatocellular carcinoma cell lines (HepG2). Hum Exp Toxicol 2016; 35: 1252–1263.

19.

Chen

Niu

Qiao

, et al. Characterization of interplay between autophagy and ferroptosis and their synergistical roles on manipulating immunological tumor microenvironment in squamous cell carcinomas. Front Immunol 2021; 12: 739039.

20.

Kim

Chen

Cheng

, et al. Pubchem in 2021: new data content and improved web interfaces. Nucleic Acids Res 2021; 49: D1388–d1395.

21.

Daina

Michielin

Zoete

. Swisstargetprediction: updated data and new features for efficient prediction of protein targets of small molecules. Nucleic Acids Res 2019; 47: W357–w364.

22.

Fang

Dong

Liu

, et al. HERB: a high-throughput experiment- and reference-guided database of traditional Chinese medicine. Nucleic Acids Res 2021; 49: D1197–d1206.

23.

Liu

Guo

Wang

, et al. BATMAN-TCM: a bioinformatics analysis tool for molecular mechANism of traditional Chinese medicine. Sci Rep 2016; 6: 21146.

24.

Uniprot: the universal protein knowledgebase in 2021. Nucleic Acids Res 2021; 49: D480–d489.

25.

Ritchie

Phipson

, et al. . Limma powers differential expression analyses for RNA-Sequencing and microarray studies. Nucleic Acids Res 2015; 43: e47.

26.

Szklarczyk

Gable

Nastou

, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res 2021; 49: D605–d612.

27.

Chin

Chen

, et al. Cytohubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol 2014; 8: S11.

28.

Wang

Han

, et al. Clusterprofiler: an R package for comparing biological themes among gene clusters. Omics 2012; 16: 284–287.

29.

Subramanian

Tamayo

Mootha

, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005; 102: 15545–15550.

30.

Liberzon

Birger

Thorvaldsdóttir

, et al. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst 2015; 1: 417–425.

31.

Berman

Westbrook

Feng

, et al. The protein data bank. Nucleic Acids Res 2000; 28: 235–242.

32.

Zhu

Kong

Niu

, et al. Mapping intellectual structure and research performance for the nanoparticles in pancreatic cancer field. Int J Nanomedicine 2020; 15: 5503–5516.

33.

Anwanwan

Singh

, et al. Challenges in liver cancer and possible treatment approaches. Biochim Biophys Acta Rev Cancer 2020; 1873: 188314.

34.

Zhang

Wang

Chen

, et al. Carbonic anhydrase 2 inhibits epithelial-mesenchymal transition and metastasis in hepatocellular carcinoma. Carcinogenesis 2018; 39: 562–570.

35.

Zhang

Zhuo

Zhou

, et al. Association between the expression of carbonic anhydrase II and clinicopathological features of hepatocellular carcinoma. Oncol Lett 2019; 17: 5721–5728.

36.

Xing

Yuan

Liu

, et al. Quantitative secretome analysis reveals clinical values of carbonic anhydrase II in hepatocellular carcinoma. Genomics Proteomics Bioinformatics 2021; 19: 94–107.

37.

Dauki

Blachly

Kautto

, et al. Transcriptionally active androgen receptor splice variants promote hepatocellular carcinoma progression. Cancer Res 2020; 80: 561–575.

38.

Zhou

Liu

Pan

, et al. Dual and opposing roles of the androgen receptor in VETC-dependent and invasion-dependent metastasis of hepatocellular carcinoma. J Hepatol 2021; 75: 900–911.

39.

Yang

Zhang

. Molecular mechanism of albumin in suppressing invasion and metastasis of hepatocellular carcinoma. Liver Int 2022; 42: 696–709.

40.

Zhan

Bai

Liao

, et al. Identification and validation of EPHX2 as a prognostic biomarker in hepatocellular carcinoma. Mol Med Rep 2021; 24: 650.

41.

Kim

Min

, et al. Targeted proteomics predicts a sustained complete-response after transarterial chemoembolization and clinical outcomes in patients with hepatocellular carcinoma: a prospective cohort study. J Proteome Res 2017; 16: 1239–1248.

42.

Yang

Sang

. Molecular subtypes based on immune-related genes predict the prognosis for hepatocellular carcinoma patients. Int Immunopharmacol 2021; 90: 107164.

43.

Zou

Xiao

Huang

, et al. Survey of the translation shifts in hepatocellular carcinoma with ribosome profiling. Theranostics 2019; 9: 4141–4155.

44.

Shen

Yin

Zhou

, et al. Combined inhibition of AURKA and HSF1 suppresses proliferation and promotes apoptosis in hepatocellular carcinoma by activating endoplasmic reticulum stress. Cell Oncol (Dordr) 2021; 44: 1035–1049.

45.

Bioulac-Sage

Cubel

Balabaud

, et al. Revisiting the pathology of resected benign hepatocellular nodules using new immunohistochemical markers. Semin Liver Dis 2011; 31: 91–103.

46.

Sun

Zhou

, et al. LncRNA RP11-422N16.3 inhibits cell proliferation and EMT, and induces apoptosis in hepatocellular carcinoma cells by sponging miR-23b-3p. Onco Targets Ther 2019; 12: 10943–10961.

47.

Kim

Lee

, et al. High-Speed screening of lipoprotein components using online miniaturized asymmetrical flow field-flow fractionation and electrospray ionization tandem mass spectrometry: application to hepatocellular carcinoma plasma samples. Anal Chem 2021; 93: 4867–4875.

48.

Tang

, et al. CD73 Sustained cancer-stem-cell traits by promoting SOX9 expression and stability in hepatocellular carcinoma. J Hematol Oncol 2020; 13: 11.

49.

Chen

Qian

Sun

, et al. Sox9/INHBB axis-mediated crosstalk between the hepatoma and hepatic stellate cells promotes the metastasis of hepatocellular carcinoma. Cancer Lett 2021; 499: 243–254.

50.

Morota

Nakagawa

Sekiya

, et al. A comparative evaluation of Golgi protein-73, fucosylated hemopexin, α-fetoprotein, and PIVKA-II in the serum of patients with chronic hepatitis, cirrhosis, and hepatocellular carcinoma. Clin Chem Lab Med 2011; 49: 711–718.

51.

Lin

, et al. Choline kinase α mediates interactions between the epidermal growth factor receptor and mechanistic target of rapamycin Complex 2 in hepatocellular carcinoma cells to promote drug resistance and Xenograft tumor progression. Gastroenterology 2017; 152: 1187–1202.