Integrated network pharmacology and molecular docking approaches to reveal the synergistic mechanism of multiple components in Venenum Bufonis for ameliorating heart failure

In vivo constituents of VB

A total of 42 in vivo constituents of VB have been reported, including 17 prototypes and 25 metabolites (He et al., 2012; Liang et al., 2008; Miyashiro, Nishio & Shimada, 2008; Ning et al., 2010; Tao et al., 2017; Xia et al., 2010; Xin et al., 2016; Zhu et al., 2013). The molecular 2D files of these constituents were downloaded from the ChemSpider database (http://www.chemspider.com/) and were saved in mol format (step 1 in Fig. S1). The details were listed in Table S1.

The prediction of VB-related targets

The predicted proteins targeted by the in vivo constituents of VB were screened from MedChem Studio (MedChem Studio, 3.0; Simulations Plus, Inc, Lancaster, CA, USA, 2012). This software could efficiently capture the FDA-approved drugs that have similar chemical structures to the components in TCM (Yu et al., 2016). We picked out VB compound-drug pairs with high confidence scores (≥0.6) and considered the target proteins of the known drugs as the VB-related targets. Other parameters were set as the default values (step 1 in Fig. S1).

Known therapeutic targets for HF

In this study, two databases were employed in acquiring pathological targets for HF, by use of “heart failure” as the query. One was the DrugBank database (http://www.drugbank.ca/, version 5.1.1), which could provide detailed information on the drug targets and their links with human diseases (Griesenauer, Schillebeeckx & Kinch, 2019). Only drug–target interactions for FDA-approved anti-HF drugs and potential human protein targets associated with HF were selected for further analysis. The second platform was the Online Mendelian Inheritance in Man (OMIM) database (http://www.omim.org/, updated on May 4, 2018), a constantly updated database of human genetic diseases and genes (step 1 in Fig. S1) (Hamosh et al., 2005).

Protein–protein interaction (PPI) data

The PPI information related to the putative targets of VB constituents and known therapeutic targets for HF was harvested by use of STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database (http://string-db.org/). This database could provide a global perspective of proteins and their functional interactions and associations (Jensen et al., 2009). Results were limited to “Homo sapiens” and protein interactions with a confidence score greater than or equal to 0.4 would be selected. Other parameters were set as the default values (step 2 in Fig. S1).

Network construction and analysis

In order to illustrate the interaction among the ingredients, targets and diseases, a “ingredient–target–disease” network was built by introducing the information of candidate compounds of VB, putative targets of VB and HF-associated targets into Cytoscape software (version 3.6.0, Boston, MA, USA). This software is an efficient open source bioinformatics tool for visualizing and analysing the complex biological networks (Shannon et al., 2003). Three important topological characteristics (“degree”, “betweenness”, and “closeness”), which have been described in our previous publication (Yu et al., 2018), were calculated to assess the central attribute of each hub node in the network by use of the Cytoscape plugin “Network Analyzer”. And “degree” > median degree centrality, “betweenness” > median betweenness centrality and “closeness” > median closeness centrality were adopted as the screening criteria to acquire the critical targets (Luo et al., 2020). Other parameters were set as the default values (step 2 in Fig. S1).

Pathway enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis and Gene Ontology (GO) enrichment analysis were undertaken to explore the potential functions of the pivotal target proteins involved in the VB-mediated treatment of Hf (Nguyen et al., 2019a; Nguyen et al., 2019b), by use of the Database for Annotation, Visualization and Integrated Discovery (DAVID) system (http://david.abcc.ncifcrf.gov/home.jsp/, v6.8) (Dennis Jr et al., 2003). Relevant pathways with the false discovery rate (FDR)-corrected P-value < 0.05 were considered statistically significant. Other parameters were set as the default values (step 3 in Fig. S1).

Molecular docking simulation

Molecular docking studies were conducted to further verify the reliability of the potential targets using CDOCKER module implemented in Discovery Studio 2016 (DS 2016). CDOCKER is a semi-flexible molecular docking analysis method based on the CHARMm force field, which can produce high-precision docking results, and it provides information on the interaction binding energy and ligand–receptor docking mode. The three-dimensional (3D) structures of the candidate compounds were generated using Chem3D Pro 12.0 and the crystallographic structures of the proteins encoded by the candidate target genes were obtained from the PDB database (http://www.rcsb.org/pdb/home/home.do), which was then decorated by removing the ligands, adding hydrogen, removing water, optimizing and patching amino acids. The binding site was defined by the ligand atoms, and the radius range was automatically generated. After each compound was docked, the 10 best conformations were obtained (Yang, 2020). Finally, CDOCKER interaction energies (CIEs) were used to assess the binding affinities between the core targets and the corresponding compounds. And the conformation corresponding to the lowest CIE was selected as the most probable binding conformation. All parameters used in calculation were default except for explained (step 4 in Fig. S1).

Putative targets for VB

As shown in Table S2, a total of 260 potential targets of the in vivo components of VB were obtained by MedChem Studio. The results showed that the candidate compounds could act on multiple targets, and one target could also be linked to multiple components.

Known therapeutic targets for HF

We collected 109 and 199 HF-related targets from DrugBank database and OMIM database, respectively. We then checked the data and eliminate redundant entries, leaving a final dataset of 292 targets associated with HF (Table S3). Among them, 22 targets were both the VB- and HF-related targets, including NR1I2, MT-CO1, CA2, NR3C1, SHBG, ATP1A1, CYP17A1, PGR, NR3C2, MME, AR, ACE, PRKCA, CYP11B2, RXRA, PPARG, PIK3CG, PPARD, PDE4B, KCNMA1, MED1 and HIF1A.

Network and pathway analysis

To facilitate scientific interpretation of the complex relationships between VB and HF, a “chemical target-disease associated gene” network, comprising the VB-related targets and the HF-related targets, was constructed based on the PPI data from STRING database. As listed in Table S4, this network was composed of 461 nodes and 5009 edges. After computing the values of the topological features of all hubs, 158 major hubs were identified because they satisfied the criteria (degree cutoff = 19, betweenness cutoff = 0.001659, and closeness cutoff = 0.381426). Among them, 93 hubs were VB-related targets, 65 hubs were HF-related targets. These major hubs may play a critical role in the entire interaction network. The specific information about the major hubs was shown in Table S5.

Potential mechanisms of VB in treating HF

To elucidate the biological process (BP), molecular function (MF), cellular components (CC) of these major hubs involved in, GO enrichment analysis was conducted on the major hubs. As shown in Figs. 2A–2C, the top 10 significant GO entries (P < 0.05) were selected on the basis of P value.

To deeply determine the function and systemic association of the main hubs, KEGG pathway enrichment analysis were conducted. The results indicated that these hubs were enriched in 101 significant pathways (P < 0.05). As shown in Fig. 2D, the top 11 pathways were considered as the main biological processes involved in the treatment. To further identify the functional mechanisms of VB on HF, the candidate compounds-main VB-related targets-main pathways network diagram was generated and elucidated in Fig. 3. The hubs can be mainly divided into the following three functional modules (circles): signal transduction (including adrenergic signalling in cardio-myocytes, calcium signalling pathway, cAMP signalling pathway, dilated cardiomyopathy, long-term potentiation and cGMP-PKG signalling pathway), cardiovascular system (including cardiac muscle contraction, vascular smooth muscle contraction and hypertrophic cardiomyopathy), and neural regulation (including dopaminergic synapse and circadian entrainment). Interestingly, adrenergic signalling in cardio-myocytes was highly enriched in KEGG pathway analysis, which played a critical role in the regulation of cardiac muscle contraction (the top-ranked GO: Biological Process terms), suggesting that VB may impart therapeutic effects on HF majorly through adrenergic signalling in cardio-myocytes.

As shown in Table S6, the VB putative targets associated with adrenergic signalling in cardio-myocytes include guanine nucleotide-binding protein G(s) subunit alpha (GNAS), adenylate cyclase 2 (ADCY2), ADCY5, protein phosphatase 1 catalytic subunit gamma (PPP1CC), protein phosphatase 2 catalytic subunit alpha (PPP2CA), protein phosphatase 2 catalytic subunit beta (PPP2CB), calcium voltage-gated channel subunit alpha1 C (CACNA1C), CACNA1D, ATPase Na⁺/K⁺ transporting subunit alpha 1 (ATP1A1), mitogen-activated protein kinase 1 (MAPK1), protein kinase C alpha (PKCα, encoded by PRKCA). Figure 4 depicts a graphical overview of adrenergic signalling in cardio-myocytes influenced by main putative targets of VB components.

Molecular docking

Table 1 displayed the CIEs (top 5 for each target) of VB hit constituents against the active sites of the screened targets in adrenergic signalling, including ATP1A1, GNAS, ADCY2, ADCY5, PPP1CC, PPP2CA, PPP2CB, CACNA1C, CACNA1D, MAPK1 and PRKCA. The results indicated that the VB-related components had been docked successfully with ATP1A1, GNAS, MAPK1 and PRKCA, which may be the key targets involved in VB for the treatment of HF. Herein, we selected four representative pairs of binding interactions to illustrate how the four targets bound to their corresponding components (Fig. 5). The interplay between ATP1A1 and bufotalin was depicted in Figs. 5A and 5B. The hydroxyl groups on bufotalin could form three hydrogen bonds with SER209, ARG191 and VAL712. Another key residue which involved in interaction was MET157. The binding mode of GNAS and cinobufaginol was depicted in Figs. 5C and 5D. The hydroxyl groups could bind with LEU171, MET255, ASN254 and LYS300 by forming hydrogen bonds. The carbonyl on the lactone ring could also form hydrogen bond with LYS305. In addition, the lactone ring could bind with GLU164 and LYS305 via pi-anion and pi-cation interactions. Other interactions including alkyl and pi-alkyl were connected with ALA303, TYR163 and LEU296. The action mode of MAPK1 and cinobufaginol was depicted in Figs. 5E and 5F. Cinobufaginol could form five hydrogen bonds with LYS112, LYS52 and TYR34. In addition, the lactone ring could bind with ASP109 via pi-anion interaction. Other interactions including alkyl and pi-alkyl were connected with VAL37. The interplay between PRKCA and 19-oxo-bufalin was depicted in Figs. 5G and 5H. The hydroxyl groups of 19-oxo-bufalin could form two hydrogen bonds with PRO202 and LYS230. Another key residue which involved in interaction was LEU200.