3.1. Network pharmacology prediction
3.1.1. Potential targets identification of OCA for treating VPA-induced hepatotoxicity
As shown in Figure 1A, 462 targets related to VPA-induced hepatotoxicity were obtained from CTD database. From CTD, TargetNET, SwissTarget Prediction and Pharmmapper databases, a total of 288 targets related to OCA were screened. There were 81 overlapping targets, through which OCA might ameliorate VPA-induced hepatotoxicity.
In order to identify the potential key targets, the PPI information of intersecting targets was constructed by STRING database and Cytoscape 3.9.1. As shown in Figure 1B, the PPI network composed 80 nodes with 460 connections. After triplicate screenings (Degree > 7.5, betweenness centrality > 0.005553, and closeness centrality > 0.422508), 30 nodes were identified as potential key targets of OCA in the treatment of VPA-induced hepatotoxicity (Supplementary Table 1).
3.1.2. KEGG pathway and GO enrichment analysis
In order to explore the possible mechanisms of OCA in the treatment of VPA-induced hepatotoxicity, the 81 intersection targets were evaluated by a KEGG pathway analysis. These targets were highly enriched in 126 pathways (P < 0.05), and many of which were related to drug-induced liver injury. As shown in Figure 2A, the top 10 pathways related to VPA-induced hepatotoxicity were non-alcoholic fatty liver (NAFLD) disease, IL-17 signaling pathway, bile secretion, FoxO signaling pathway, PPAR signaling pathway, PI3K-Akt signaling pathway, TNF signaling pathway, AMPK signaling pathway, MAPK signaling pathway and apoptosis.
GO enrichment analysis was carried out to explore the potential roles of the target genes. In this study, 81 intersection targets were enriched in 267 GO items, including 197 biological processes (BP), 49 cellular components (CC), and 21 molecular functions (MF). Based on the lowest P value, the top 5 BP items were negative regulation of apoptotic process, positive regulation of transcription from RNA polymerase II promoter, response to ethanol, positive regulation of transcription, DNA-templated, response to xenobiotic stimulus, and positive regulation of pri-miRNA transcription from RNA polymerase II promoter. The top 5 CC items were nucleoplasm, chromatin, macromolecular complex, cytoplasm and nucleus. The top 5 MF items were RNA polymerase II transcription factor activity, ligand-activated sequence-specific DNA binding, enzyme binding, sequence-specific DNA binding, zinc ion binding, and identical protein binding. The top 10 items with lowest P value in every category were shown as a histogram (Figure 2B).
3.1.3. The main target-pathway network of OCA treatment for VPA-induced hepatotoxicity
According to the targets and top 10 KEGG pathways, a target-pathway network was constructed by Cytoscape 3.9.1. As shown in Figure 3, the network contained 92 nodes and 206 edges, fully revealing the multi-target mechanisms of OCA in treating VPA-induced hepatotoxicity.
3.2. GEO data analysis
Dataset GSE138810 provided the information on the effect of OCA treatment on hepatic gene expression in mice treated with VPA. By comparing the hepatic gene expression profiles between OCA and VPA+OCA group, 80 DEGs were identified with 47 up-regulated and 33 down-regulated (Figure 4). As shown in Figure 5A, these DEGs were highly enriched in bile secretion, metabolic pathways, arachidonic acid metabolism, PPAR signaling pathway, retinol metabolism, thyroid cancer, endometrial cancer, gastric cancer and melanoma. Among them, bile secretion and PPAR signaling pathway were also identified in the analysis of network pharmacology, which would be further explored.
As shown in Figure 5B, GO enrichment analysis revealed that many biological processes were involved in the effect of OCA on hepatic gene expression profiles of mice under chronic treatment with VPA, such as lipid metabolic process, transmembrane transport, cellular response to lithium ion, long-chain fatty acid metabolic process, fumarate transport and so on. These DEGs primarily localized in basolateral plasma membrane, apical plasma membrane, membrane, cell surface and extracellular space, and functioned as transmembrane transporter activity, heme binding, monooxygenase activity, oxidoreductase activity and oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular.
3.3. Molecular docking
To validate whether bile secretion and PPAR signaling pathway were involved in the treatment of VPA-induced hepatotoxicity with OCA, the potential key targets (e.g., RXRA/RXRα, CYP7A1, PPARA/PPARα, PPARG/PPARγ, CYP3A4, NR1H4/FXR and NR0B2/SHP) involved in two signaling pathways were selected for molecular docking with OCA. As shown in Table 1, free binding energies of docking results were in the range of +4.81 to -8.05 kcal/mol. OCA exhibited strong binding activity with CYP3A4, FXR and CYP7A1, as well as good binding activity with PPARγ, RXRα, and PPARα; however, it showed poor binding activity with SHP. Therefore, except for SHP, the other six proteins were the key targets in the OCA treatment of VPA-induced hepatotoxicity. The interactions and binding modes between six proteins (CYP3A4, FXR, CYP7A1, PPARγ, RXRα, and PPARα) and OCA were shown in Figure 6A-F. According to the chemical structure of OCA, three hydroxyl groups enabled these proteins to form conventional hydrogen bond with OCA; one carboxyl group enabled these proteins which contained lysine or arginine in the binding pockets to form attractive charge and / or salt bridge with OCA. In addition to these interactions, van der Waals forces and hydrophobic interactions also contributed to the binding affinities.
Table 1 Details of targets and OCA for molecular docking
Target
|
KEGG Pathway
|
PDB ID
|
Binding Energy (kcal/mol)
|
CYP3A4
|
Bile secretion
|
4D6Z
|
-8.05
|
FXR
|
Bile secretion
|
1OSH
|
-7.63
|
CYP7A1
|
Bile secretion, PPAR signaling pathway
|
3V8D
|
-7.01
|
PPARγ
|
PPAR signaling pathway
|
5Y2T
|
-6.47
|
RXRα
|
Bile secretion, PPAR signaling pathway
|
6LB4
|
-6.39
|
PPARα
|
PPAR signaling pathway
|
6KB1
|
-6.15
|
SHP
|
Bile secretion
|
6W9M
|
+4.81
|
3.4. Molecular dynamics simulations
Based on the docking results, the complexes of OCA with CYP3A4 (4D6Z), CYP7A1 (3V8D), PPARγ (5Y2T), and PPARα (6KB1) were subjected to molecular dynamics simulation analysis to assess their dynamics.
RMSD measures the extent of deviation in the position of atoms relative to their initial positions, with a lower deviation suggesting greater conformational stability. As shown in Figure 7A, the RMSD curves of all complexes initially fluctuated during the simulations, and gradually converged to a stable state below 0.5 nm. However, in comparison with other complexes, the RMSD curve of PPARγ-OCA exhibited dramatic fluctuation in the 100 ns simulation trajectory. The Rg parameter quantifies the compactness of the complexes, with a lower value indicating greater levels of compaction. The Rg curves of OCA with CYP3A4, CYP7A1, PPARγ and PPARα in Figure 7B exhibited slight fluctuations while remaining stable throughout the simulations, with averages of 2.296 nm, 2.318 nm, 2.541 nm and 1.771 nm, respectively. RMSF analysis is performed to assess the flexibility of amino acid residues within a molecular system during a molecular dynamics simulation. The RMSF curves indicated that most of the protein amino acid residues in complexes of OCA with CYP3A4, CYP7A1 and PPARα exhibited slight fluctuations, measuring less than 0.5 nm, while the PPARγ-OCA complex displayed more pronounced fluctuation ranging from 0.6 to 1.2 nm (Figure 7C). The SASA parameter quantifies the surface area of a molecule exposed to the solvent molecules during simulations. The SASA plot (Supplementary Figure 1) revealed that slight fluctuations in the curves of OCA with CYP3A4, CYP7A1 and PPARα, while the curve of OCA with PPARγ fluctuated dramatically before eventually stabilizing, which was consistent with the results of RMSD, Rg and RMSF. Additionally, the hydrogen bonds information of complexes was also acquired during the simulations. As shown in Figure 7D, the number of hydrogen bonds formed between OCA and CYP3A4, CYP7A1, PPARγ and PPARα were 5, 6, 1 and 0 respectively at 100 ns.
By integrating and analyzing the aforementioned results of molecular dynamics simulations, it was found that there were no significant conformational changes in the proteins during the simulations, indicating stable dynamic interactions between OCA and CYP3A4, CYP7A1, PPARγ and PPARα.