Oral bioavailability (OB) and drug-likeness (DL) are the two main parameters of the absorption, distribution, metabolism, and excretion (ADME) model (Barton et al., 2006) that affect drug absorption throughout the gastrointestinal tract. Hence, by using “Renshen,” “Huanglian,” and “Chuanxiong,” as keywords, we identified the active ingredients of the YZQX formula, screening for compounds with oral bioavailability (OB) ≥30% and drug-likeness (DL) ≥0.18 in the Traditional Chinese Medicine Systems Pharmacology (https://tcmspw.com/tcmsp.php) (Ru et al., 2014).

After screening the active ingredients of the YZQX formula, the next significant process was to identify their corresponding targets. The putative targets of the active ingredients in the YZQX formula were identified and collected using the DrugBank database (https://www.drugbank.ca/) (Law et al., 2014). Further, all of these obtained targets were limited to the species “Homo sapiens” in the Uniprot database (https://www.uniprot.org/) (2017). In addition, we used the Cytoscape 3.7.2 software to construct and present a Herb-ingredient-target network of the YZQX formula. In this network, nodes represent formulas, herbs, ingredients, and targets, while the edges represent their interaction relationships.

The microarray data of GSE63061 and GSE98793 were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) using their microarray platforms, GPL10558 (Illumina HumanHT-12 V4.0 expression bead chip) and GPL570 (Affymetrix Human Gene Expression Array), respectively. Among them, the GSE63061 dataset contained 139 AD samples and 134 normal samples, and the GSE98793 dataset contained 128 MDD samples and 64 normal samples. Thereafter, the probes were converted into gene symbols according to the annotation information of the corresponding platform. The DEGs between disease (AD or MDD) groups and normal groups were identified and screened as implemented in the R software (Version 3.6.0) limma package with a threshold value of p < 0.05. DEGs with p < 0.05, and |log2 fold change (FC)| > 0.1 were then presented visually using volcano plots. Venn diagrams were used to further obtain the overlapping genes from the AD and MDD microarray datasets.

To explore all the possible interactions among genes, the PPI network of YZQX-related targets and AD-MDD overlapped targets were established using the Cytoscape plugin BisoGenet (Martin et al., 2010), which integrates PPI data from the following six databases: the Biological General Repository for Interaction Datasets (BioGRID), the Biomolecular Interaction Network Database (BIND), the Database of Interacting Proteins (DIP), the Molecular INTeraction Database (MINT), the Human Protein Reference Database (HPRD), and the IntAct Molecular Interaction Database (IntAct).

Building on the previous process, two PPI networks were intersected into a merged network for further analysis using Cytoscape 3.7.2 software. Next, the CytoNCA (Tang et al., 2015) plugin of Cytoscape was used to calculate and evaluate the topological significance of each node in the network based on these parameters of degree centrality (DC), betweenness centrality (BC), closeness centrality (CC), eigenvector centrality (EC), local average connectivity-based method (LAC), and network centrality (NC). A previous study (Tang et al., 2015) specified in detail the definitions and computing methods of these parameters; in brief, a higher value represents greater importance of the node in the merged network. We calculated the degree of centrality in a merged network. It is worth mentioning that this gene is considered “a big hub” when the degree centrality of one node is more than twice the median degree centrality in a merged network. Subsequently, a sub-network was extracted based on two times the median degree centrality and used to perform the next screening. The median values of the other five parameters served as threshold values to further identify the hub genes and core network. Finally, after two screenings, the core network of YZQX treating comorbidity with AD and MDD was obtained.

In addition, we also further analyzed the attributes of network screened by CytoNCA using MCODE plugins in Cytoscape software.

ClusterProfiler (Yu et al., 2012), a new ontology-based package of R software version 3.6.0, which assigned the GO and Kyoto Encyclopaedia of Genes and Genomes (KEGG) database, was used for enrichment analysis. Of note, the parameter was set to a p-value cutoff <0.05, which was deemed to be significant. The GO term contained three categories: biological process (BP), cellular component (CC), and molecular function (MF). The results of GO, sorted from small to large according to the p-value as a parameter, and the top 10 terms of each category were selected using a bubble plot for visualization. For the top five terms, a sub-network was then established in Cytoscape software version 3.7.2 to visualize the relationships between genes and GO terms. Furthermore, based on the genes of the core network of YZQX in AD and MDD treatment, KEGG pathway analysis was performed, and a bubble plot was plotted.

To explore the possible docking modes and binding affinities between macromolecules and small molecules (ligands) in-depth, we performed molecular docking based on molecular modeling techniques using AutoDock software (Yang et al., 2016). We subsequently identified the direct-acting targets in the core network and their corresponding bioactive ingredients for molecular docking. The X-ray crystal structures of key targets were obtained from the RCSB Protein Data Bank (PDB, https://www.rcsb.org/) and saved in PDB format. The MOL2 formats of the bioactive ingredients were downloaded from the ZINC database (http://zinc.docking.org/). Before molecular docking, we used AutoDock software (version 4.2.6) to remove the ligand, hydrogenate, calculate the charge, add protein type for the macromolecules, and saved it in PDBQT format. Subsequently, we ran the AutoDock to perform molecular docking between macromolecules and ligands based on the default parameters. Blinding energy of less than “−5” represented a better binding interaction between macromolecules and small molecules (Gaillard, 2018). Finally, the PyMOL software (version 2.4.1) was used to visualize the docking results.

A total of 43 ingredients in YZQX were screened from the TCMSP database with the criteria of OB ≥ 30% and DL ≥ 0.18. These were defined as candidate bioactive ingredients, of which seven were from Chuanxiong, 14 from Huanglian, and 22 from Renshen. The detailed characteristics, including molecular ID, molecular name, structure, OB, and DL, of the selected bioactive ingredients are presented in Table 1. We retrieved the DrugBank database using candidate bioactive ingredients and identified 504 ingredient-related targets. From each of the Chinese herbs in YZQX, 39, 251, and 214 potential targets were obtained, respectively (Supplementary Table S1). After removing the reduplicated targets, the total number of putative targets associated with bioactive ingredients was adjusted to 193. In addition, these targets were annotated as gene symbols in the UniProt database.

Subsequently, using Cytoscape software, we established a network indicating the relationship among Chinese herbs, ingredients, and putative targets, as shown in Figure 2. This ingredient-target network contained 231 nodes and 541 edges. Quercetin was the most prolific, with 141 potential targets, followed by kaempferol (56), and beta-sitosterol (28).

We downloaded the AD expression datasets and selected DEGs by calculating the difference in gene expression between 134 control samples and 139 AD samples. Ultimately, 4,873 DEGs, including 2,104 downregulated and 2,769 upregulated genes, were identified from the GSE63061 dataset using the limma package in R software. Additionally, based on the cutoff criterion of p < 0.05, a total of 1,707 DEGs, containing 889 downregulated and 818 upregulated genes, were extracted and analyzed from the microarray data GSE98793. The DEGs with p < 0.05, and |log2 FC| > 0.1 in AD or MDD were visualized in a volcano plot (Figures 3A,B). And the detailed information of AD- and MDD-related DEGs was presented in Supplementary Table S2. Moreover, the overlapping genes of AD and MDD were 378, as shown in Figure 3C.

In the present study, we used the BisoGenet plugin based on the six PPI databases to generate PPI networks of YZQX putative targets and comorbidity with AD- and MDD-related targets (Figure 4). As presented in Figure 4A, the PPI network of AD with depression-related targets consisted of 162,057 edges and 6,435 nodes. The PPI network of putative targets associated with YZQX contained 154,133 edges and 6,322 nodes (Figure 4B).

To explore the pharmacological mechanisms of YZQX in treating AD and MDD comorbidities, we utilized the merge network function supplied by Cytoscape software to generate a new merged network, with 119,660 edges and 4,122 nodes, consisting of overlapping genes from the two networks mentioned above (Figure 4C). After DC calculation, 946 nodes and 43,494 edges with a DC of more than two-fold the median were selected for integration into the subnetwork (Figure 4D). Next, these genes were further screened based on topological property analysis of DC, BC, CC, EC, LAC, and NC using the CytoNCA plugin in the Cytoscape software. The median values showed BC values of 394.718, CC of 0.498, EC of 0.021, NC of 21.359, LAC of 19.469, and the two-fold median DC of 236 (Figure 4E), and finally, a core network including 92 targets together with 1,578 edges was identified for YZQX against AD and MDD comorbidities. The detailed topology parameter of the core network involving 92 targets is presented in the Supplementary Table S3.

According to the ranking of MCODE score, the largest cluster has 25 nodes and 215 edges, with a score of 17.917 points; The smallest module has three nodes and three edges, with a score of three points. And the clusters are presented in Supplementary Figure S1.

To reveal the biological mechanisms of the 92 crucial targets, the cluster profile package in R software was used to perform GO and KEGG enrichment analysis, yielding a total of 1,212 BPs, 95 CCs, 119 MFs, and 117 KEGG pathways (p < 0.05) (Supplementary Tables S4, S5). As indicated in Figure 5, the highly enriched GO terms in BP were involved in DNA-binding transcription factor activity regulation, DNA metabolic process regulation, protein deubiquitination, protein modification by small protein removal, and the DNA biosynthetic process. The main CCs of the 92 core targets were related to nuclear chromatin, focal adhesion, cell-substrate adherens junction, cell-substrate junction, and ubiquitin ligase complex. MFs focused on ubiquitin-protein ligase binding, ubiquitin-like protein ligase binding, phosphatase binding, ubiquitin-like protein transferase activity, and RNA polymerase II transcription factor binding. In addition, based on the KEGG pathway data analysis, we found that the PI3K-Akt signaling pathway, ubiquitin-mediated proteolysis, estrogen signaling pathway, MAPK signaling pathway, and TGF-beta signaling pathway concurrently occupied the process of YZQX treatment of comorbidity with AD and MDD (Figure 5E).

The PPI network indicated that HSP90AA1, ESR1, AKT1, VCAM1, EGFR, CDK1, MAPK1, CDK2, MYC, HSPB1, and HSPA5 are pivotal targets in the core network mentioned above, and these nodes were also overlapping genes of AD and MDD comorbidity. In accordance with key targets, we worked backward to their corresponding bioactive ingredients serving as ligands, thus performing validation of molecular docking. In total, 2 ingredients were docked with AKT1, 2 with CDK1, 3 with CDK2, 6 with ESR1, 11 with HSP90AA1, 1 with EGFR, 1 with HSPA5, 1 with HSPB1, 1 with MAPK1, 1 with MYC, and 1 with VCAM1. The ligand-receptor binding conformation had the lowest binding energy, indicating that this conformation is the most stable, and thus the possibility of interaction is the greatest. The top three docking combinations for the target macromolecules and corresponding ingredients included Fumarine and HSP90AA1, Worenine and ESR1, berberine and HSP90AA1, with binding energies of −6.83, −6.6, and −6.26 kcal/mol, respectively. The other specific docking scores of the ingredients with targets are shown in Supplementary Table S6. If the docking score was less than −6 kcal/mol and ranked as the top 4, PyMOL was used to display the three-dimensional view of the docking mode of the ingredients and the target (Figure 6).