Molecular docking and network connections of active compounds from the classical herbal formula Ding Chuan Tang

DCT-compounds identification

DCT-compounds were identified by searching the keywords DCT and each of the nine DCT herb names (Chinese, English (common), pharmaceutical and Latin) within five English databases (PubMed, Embase, Science direct, HIT database and CINAHL) and two reference books (Traditional Chinese medicines: Molecular structures, natural sources and applications (Zhou, Yan & Xie, 2003), Chinese herbal medicine: Materia medica (Bensky, Clavey & Stoger, 2004)). The three-dimensional chemical structures of the DCT-compounds were downloaded from PubChem, and DS 4.5 was used for further optimization (step 1 in Fig. 1).

Pharmacophore screening with DCT-compounds

Pharmacophore models are constructed from multiple pharmacological features and targets and are used in computational drug discovery to recognize the molecular features of one or more compounds (ligand) required for a lock and key fit with macromolecules (proteins) with the same or similar biological activity (Thangapandian et al., 2011; Wieder et al., 2017). The DS 4.5 software contains a Ligand pharmacophore profiler (LPP) protocol, allowing for virtual screening of multiple chemical compounds against a pharmacophore database (PharmaDB) containing 117,423 pharmacophore models and constructed from 7,028 protein databank (PDB) protein-ligand X-ray crystal structures. This protocol investigates the potential biological actions for the ligands using the pharmacophore models, with the results for the docking within the LPP calculated in fit-values, as a high fit-value indicates the likely comparability of small molecules and target protein binding (Ming, Yi & Yang, 2018). The identified 396 DCT-compounds were virtually screened with the PharmaDB database and the inclusion criteria was set as fit-value 0.9 and above (1.0 fit-value being perfect prediction), to identify the target proteins with the high comparability score (step 2 in Fig. 1).

Identification and screening of genes for asthma

Genes associated with the asthma disease pathways were identified using “asthma” as a keyword to search the following four databases: Disease gene search engine with evidence (Kim et al., 2013), Genetic Association Database (Becker et al., 2004), Pubmed_Gene (Benson et al., 2013) and GeneCards (Stelzer et al., 2016) (step 3 in Fig. 1). The identified asthma genes were further screened with the pharmacophore and DCT-compound docking results to identify mutual target proteins (step 4 in Fig. 1).

Molecular docking of DCT-compounds and potential target proteins

DCT-compounds that resulted with fit-value scores ≥0.9 were docked using AutoDock Vina software (Trott & Olson, 2010) with the target proteins identified in the asthma gene and pharmacophore screening (step 5 in Fig. 1). Molecular docking included three phases: (1) crystal complex identification, (2) crystal complex preparation, and (3) molecular docking.

Establishment of networks

Cytoscape 3.6.0 software (Shannon et al., 2003) was used to establish the networks between docked DCT-compounds and the identified target proteins (step 6 in Fig. 1). The DCT-compounds and target proteins were defined as the nodes of the network, while the affinity binding energy value was the edge of the network. The Network analysis function in Cytoscape 3.6.0 analyzed the features of the network structure through the yfiles plug-in module to draw the network image. This further optimized the layout of the network and identified the most suitable protein targets for DCT-compounds. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa et al., 2017) includes an advanced mapping tool (Search and Color pathway) linked to Cytoscape 3.6.0 software and accessing the affinity binding results from DCT-compounds and protein docking, identified the signaling pathways targeted by DCT.

Data provided from the Protein Data Bank database (PDB) aided in the establishment of network connections between the DCT-compounds and the target proteins post the molecular docking process. Furthermore, the KEGG database (Kanehisa et al., 2017), UniProt database (UniProt Consortium, 2018), Wikipathways (Slenter et al., 2018) and Cell signaling technology pathways (Cell Signaling Technology, 2018) were searched to define and interpret the identified network pathways for the potential mechanisms of action and relevance to the management of asthma.

DCT-compound identification

Following the literature search strategies, the original search identified 644 DCT-compounds; however, 109 duplicates and one toxic compound (cyanide) were removed and a further 138 were excluded due to lack of compound identification data such as the molecular structure, molecular formula, simplified molecular-input line-entry system data and physicochemical property data. This resulted in 396 DCT-compounds included for analysis, with 40 from Ginkgo biloba L. (Bai Guo), 40 from Ephedra sinica Stapf (Ma Huang), 30 from Perilla frutescens (L.) Britton (Su Zi), 100 from Glycyrrhiza uralensis Fisch (Gan Cao), 15 from Tussilago farfara L. (Kuan Dong Hua), 17 from Prunus armeniaca L. (Xing Ren), 43 from Morus alba L. (Sang Bai Pi), 33 from Scutellaria baicalensis Georgi (Huang Qin) and 78 from Pinellia ternata (Thunb.) Makino (Ban Xia).

Pharmacophore screening with DCT-compounds

The virtual screening of the 396 DCT-compounds with the 117,423 pharmacophore models within the PharmaDB resulted in 159 pharmacophore models (potential target proteins). Following this initial hit list, this step resulted in two further sets of outcomes. The first outcome was the reduction of the original hit list of 159 potential target proteins down to 53 due to the inclusion parameter set at fit-value 0.9 and above (Table S1). The second outcome was the reduction of 396 DCT-compounds to 243 DCT-compounds, also due to the 0.9 and above fit-value inclusion parameter (Table S2). This second outcome narrowed down the number of compounds compared during the molecular docking phase to only those DCT-compounds with high affinity binding; as a high fit-value indicates high comparability of small molecules and target protein binding (Ming, Yi & Yang, 2018).

Identification of genes for asthma

Asthma is classified as a “complex” heritable disease with multiple genes implicated and contributing to its development (Global Initiative for Asthma, 2018). A total of 234 asthma genes were identified, with the results listed in Table S3.

Screening for common target proteins for DCT target proteins and asthma genes

The virtual screening of 53 DCT potential target proteins identified from pharmacophore docking and the 234 identified asthma genes, resulted in seven mutual targets as centrally displayed in Fig. 2, with the list of all proteins provided in Tables S1 and S3. These mutual protein targets included Adenosine A_2A receptor (ADORA2A), Estrogen receptors alpha (ESR1), Kinase insert domain receptor (KDR), Leukotriene A4 hydrolase (LTA4H), Macrophage migration inhibitory factor (MIF), cAMP-specific 3′,5′-cyclic phosphodiesterase 4D (PDE4D) and Peroxisome proliferator-activated receptor gamma (PPARG).

In Fig. 2, the group on the left-hand side (Group 1) is the 53 DCT protein targets resulting from pharmacophore virtual screening. Group 2: On the right are all the identified genes associated with asthma. In the center highlighted are the seven proteins both groups target: ADORA2A, ESR1, KDR, LTA4H, MIF, PDE4D and PPARG. At the time of this investigation, the active binding sites, required for Autodock Vina, could not be identified for ADORA2A and MIF; therefore, they were not included in further examination. Five target proteins (ESR1, KDR, LTA4H, PDE4D and PPARG) were included in the next step: Molecular docking with DCT-compounds.

Molecular docking of active compounds and potential target proteins

Establish network connections

Figure 3 displays the Cytoscape graphic network connections from the molecular docking results between the 243 DCT-compounds and the five target proteins. The Cytoscape 3.6.0 yfiles plug-in module was used to draw and optimize the layout of the network image (Shannon et al., 2003). All five target proteins resulted in 155 network connections with compounds from the nine DCT herbs including: Ginkgo biloba L. (Bai Guo) (n = 15), Ephedra sinica Stapf (Ma Huang) (n = 5), Perilla frutescens (L.) Britton (Su Zi) (n = 13), Glycyrrhiza uralensis Fisch (Gan Cao) (n = 63), Tussilago farfara L. (Kuan Dong Hua) (n = 5), Prunus armeniaca L. (Xing Ren) (n = 1), Morus alba L. (Sang Bai Pi) (n = 24), Scutellaria baicalensis Georgi (Huang Qin) (n = 20) and Pinellia ternata (Thunb.) Makino (Ban Xia) (n = 9) (Fig. 2). Table 2 contains the compound names for each connection node code seen in Fig. 3.

Interpreting network connections

The KEGG database resource is based on current experimental knowledge for “understanding high-level functions and utilities of the biological system” (Kanehisa et al., 2017). The KEGG asthma pathway (Fig. 4) displays the known molecular/cellular interactions and signaling pathways involved in chronic inflammation leading to asthma symptoms. As seen in Fig. 4, there are four signaling pathways involved in the KEGG asthma pathways: (1) T-Cell receptor signaling pathway; (2) B-Cell receptor signaling pathway; (3) the JAK/STAT signaling pathway and (4) the FcεRI signaling pathway. Within this discussion, these pathways will be referred to as the primary pathways. Following on, the pathways further involved within the four primary signaling pathways are referred to as the secondary pathways.

Each of the five identified target proteins, ESR1, KDR, LTA4H, PDE4D and PPARG, have established potential to play significant roles involving both the primary and secondary asthma pathways through many signaling cascade pathways.

ESR1

Nuclear hormone receptor, ESR1, in respect to the pathology of asthma, is involved in cell proliferation in tissues within the respiratory system, increasing and maintaining the number of cells in the target tissues through an increase in cell division and altering the process of cell death (UniProt). ESR1 participates in two primary KEGG asthma pathways: the B-Cell signaling pathway and the JAK-STAT signaling pathway, through the Prolactin signaling pathway (KEGG). ESR1 acts downstream of both the B-Cell signaling pathway and the JAK-STAT signaling pathway and modulates the immune response through regulation of macrophages, T and B lymphocytes and fibroblasts (KEGG). The Prolactin signaling pathway, targeted by ESR1, activates pathways involved in the inflammatory response including the KEGG secondary MAPK pathway and PI3K pathway and is also involved in the secondary KEGG asthma NF-κB signaling pathway; with ESR1 again acting downstream (KEGG). A protein–protein interaction with ESR1 and NF-κB inhibits NF-κB signaling, resulting in a decrease in the inflammatory response through a reduction of pro-inflammatory cytokines and a dampening of inflammatory gene expression in smooth muscle cells (Kovats, 2015). ESR1 also acts downstream within the Thyroid hormone signaling pathway, involved in the secondary Calcium signaling pathways and the MAPK signaling pathway through moderation of gene transcription and cell proliferation (KEGG). ESR1 participates and acts upstream in the Estrogen signaling pathway involved in the secondary KEGG asthma PI3K-Akt signaling pathway, Calcium signaling pathway and MAPK signaling pathway; targeting immune genes, regulating cell cycles and cell adhesion molecules (KEGG). As ESR1 acts upstream of the Estrogen signaling pathway, there is potential that targeting this gene may inhibit or reduce the effects of those three secondary asthma pathways.

Docking results for ESR1 included two crystal complexes classified as antagonists of ESR1 (2fja and 3erd) and one agonist of ESR1 (2yja) (Table 1). Cyclomulberrin is the single compound binding with multiple ESR1 crystals, suggesting both agonist and antagonist functions. The selective estrogen receptor modulator drug tamoxifen is a known ESR1 agonist and antagonist of ESR1, supporting the findings of cyclomulberrin (Gallo & Kaufman, 1997). The research community recognizes the disparity of inflammatory and anti-inflammatory effects of ESR1, hypothesizing the effects are dependent on the available levels of circulating estrogens or the tissue the ESR1 is present (Yung, Fuseini & Newcomb, 2018; Keselman & Heller, 2015; Dijkstra et al., 2006). Inhibiting ESR1 has been found to (1) reduce airway hyperresponsiveness and thymic development, (2) regulate and decrease pro-inflammatory mediators and (3) hinder the development and advancement of airway changes and remodeling, while promoting ESR1 has been found to impede pro-inflammatory responses in both human and animal studies (Yung, Fuseini & Newcomb, 2018; Keselman & Heller, 2015; Dijkstra et al., 2006).

KDR

KDR is classified as a type III receptor tyrosine kinase, with alternative names including vascular endothelial growth factor receptor 2 and Fetal Liver Kinase 1. KDR acts as a cell surface receptor for vascular endothelial growth factors A, C and D, all involved in the migration, growth and survival of endothelial cells (Duffy, Bouchier-Hayes & Harmey, 2013; UniProt). KDR is expressed throughout human tissues, acting upstream to the primary KEGG asthma pathway T-cell receptor signaling pathway through the Ras signaling pathway and effects cell growth, survival, cell cycle progression and gene expression (Pate et al. 2010; KEGG). KDR is involved in three of the secondary KEGG asthma pathways including PI3K-Akt signaling pathway, Calcium signaling pathway and MAPK signaling pathway; with the Ras signaling pathway additionally involved in these (KEGG). Two cascade pathways include KDR acting upstream of these three secondary asthma pathways potentially inhibiting the inflammatory outcomes; the VEGF signaling pathway- resulting in cell survival, migration and proliferation, and the RAP1 signaling pathway- involving cell adhesion, migration, proliferation, survival and gene activation (KEGG). KDR acts upstream from MAPK pathway in an additional pathway associated with focal adhesion, with effects including inflammatory cell proliferation and survival (KEGG).

All four KDR crystal complexes are classified as antagonists of KDR functions (KEGG); therefore, the DCT-compounds that docked with all four crystals provide strong evidence as KDR antagonists (Table 1). Evidence supporting the inhibition of KDR as beneficial for asthma involves a variety of pathological pathways including a dampening of the inflammatory asthma response through a reduction in recruitment of pro-inflammatory cells, a decrease circulating cytokines and cell adhesion molecules and the reduction of airway remodeling through smooth muscle hyperplasia (PDB; Uniprot).

LTA4H

The protein that the human gene LTA4H encodes is an enzyme containing both: (a) hydrolase, with the function of converting leukotriene A4 (LTA4) into leukotriene B4 (LTB4), a proinflammatory mediator and; (b) aminopeptidase, recognized as a biological characteristic of chronic obstructive pulmonary disease (Uniprot; PubMed_Gene). LTB4 promotes the inflammatory response as a mediator of leukocytes including monocytes, eosinophils, macrophages and dendritic cells, along with specific chemotactic function for neutrophils and regulating differentiation of T-cells (Di Gennaro & Haeggstrom, 2014; Liu & Yokomizo, 2015). The primary KEGG asthma FcεRI signaling pathway includes LTA4H through its function in AA metabolism (KEGG). The Eicosanoids pathway instigates both the immune and the inflammatory response through utilizing LTA4H downstream for both the FcεRI signaling pathway and the secondary AA metabolism pathway (KEGG). Pro-inflammatory cells including leukocytes, platelets and fibroblasts, in addition to target tissue of the respiratory system, widely express LTA4H (Uniprot).

Both LTA4H crystal complexes are classified as antagonists of LTA4H, therefore the 38 DCT-compounds that docked with both crystals of LTA4H (Table 1) potentially reduce asthma symptoms such as increased mucous and persistent airway inflammation, through blocking the conversion of LTA4 to LTB4 (PDB; Uniprot).

PDE4D

Nine isoforms have been identified for the PDE4D gene, each encoding hydrolyzing proteins that break down the second messenger Cyclic adenosine monophosphate (cAMP) (UniProt). The immune effects of cAMP regulate both the innate and adaptive immune responses, with many cells within the respiratory system including leukocytes, macrophages, B-cells and T-cells affected and modulated (Billington et al., 2013; Raker, Becker & Steinbrink, 2016). Increased levels of intracellular cAMP inhibit asthma pathways leading to broncho-constriction through the contraction of respiratory smooth muscle, while as a drug target, it amplifies the anti-inflammatory response and simultaneously inhibits pro-inflammatory immune and cellular factors (Billington et al., 2013; Raker, Becker & Steinbrink, 2016). The secondary KEGG asthma disease pathway PI3K-Akt signaling pathway, involves PDE4D through the cAMP signaling pathway, leading to cell growth, survival and proliferation (KEGG). PDE4D acts upstream of the PI3K-Akt signaling pathway; effectively inhibiting the outcomes and thus reducing inflammatory mediated symptoms (KEGG; Billington et al., 2013).

The four investigated PDE4D crystal complexes have all been shown to be antagonists of PDE4D functions, thus the DCT-compounds docked with them (Table 1) present the potential to reduce bronchoconstriction and reduce the asthma inflammatory cellular mediated symptoms (Billington et al., 2013; Raker, Becker & Steinbrink, 2016; UniProt). PDE4D inhibitors have been used successfully as drug targets for asthma, such as theophylline since the 1930s (Spina, 2008) and Ibudilast, an oral asthma drug approved in Asia (Rolan et al., 2008). These anti-asthma drugs present strong evidence for DCT-compounds as potential future drug candidates through computational demonstration to target PDE4D.

PPARG

The type II nuclear receptor PPARG is known primarily for the function of regulating the fatty acid pathways and glucose metabolism (UniProt). Asthma development implicates PPARG in regard to the reduction of respiratory inflammation and the permanent remodeling that occurs due to the chronic asthma response (Oh et al., 2009). The activation of PPARG moderates the initiation of pro-inflammatory Th2 cytokines, eosinophils and macrophages, while allowing the expression of anti-inflammatory cellular factors within the lung (Huang & Glass, 2010; Oh et al., 2009). PPARG is involved and acts downstream from the secondary KEGG asthma disease pathway PI3K-Akt signaling pathway through the AMP-activated protein kinase (AMPK) signaling pathway, to inhibit cell growth, cell survival and protein synthesis (KEGG).

All five PPARG crystal complexes are classified as PPARG agonists, suggesting all 108 docked DCT-compounds (Table 1) promote the reduction of pro-inflammatory Th2 cytokines, eosinophils and macrophages, along with lessening asthma chronicity (PDB; Uniprot). Rosiglitazone, a strong PPARG agonist, has shown to significantly improve lung function and decrease IL-5 and IgE in a murine model (El-Naa et al., 2015). Furthermore, a current review of rosiglitazone defines the PPARG agonist as a key regulator of airway inflammation, with strong evidence as an asthma therapeutic target (Banno et al., 2018).

ADORA2A and MIF

Although these two target proteins were excluded, they are still worth discussion due to virtual screening results. The excluded target protein ADORA2A is a receptor for adenosine, a signaling transductor and second messenger (Uniprot). Adenosine mediates G proteins and in turn activates adenylyl cyclase (an enzyme that synthesizes cAMP), associated with immune function and inflammatory pathways (Uniprot). Adenosine has long been associated with the inhibition of the inflammatory response through the effects on neutrophil and macrophage regulation (Hasko & Cronstein, 2013). The biological processes associated with ADORA2A demonstrate involvement in cell–cell signaling, the cellular defense response and the inflammatory response (Uniprot). Multiple experimental studies have shown agonist ligands of ADORA2A reduce respiratory inflammation (Van Waarde et al., 2018). This research suggests the prospect of the DCT-compounds successfully docked with ADORA2A to display similar agonistic biological effects for asthma.

MIF is a pro-inflammatory cytokine with a regulatory function of the innate immune response (Uniprot). The identified molecular functions of MIF identified and associated to asthma are (1) macrophage mediator and regulator, (2) chemokine attraction and (3) signaling receptor binding (Uniprot). The biological processes of MIF include (1) the inflammatory response, (2) cell proliferation, (3) cell surface receptor signaling pathways, (4) macrophage regulation, (5) leukocyte regulation and (6) B-cell maturation regulation and (7) the innate immune response (Uniprot). Through the direct activation of immune cells and the contribution to airway remodeling in asthma experimental studies (Lan et al., 2018; Hoi, Iskander & Morand, 2007), inhibitors of MIF prove to be a strategic therapeutic target for asthma; thus indicating the potential favorable benefits of the docked DCT-compounds.

Pharmacological contributions to asthma management

Single pharmaceutical compounds, such as the majority of asthma medications, fundamentally target single molecular mechanisms and as a result may be less efficient at treating complex diseases that have multi-cellular/tissue pathology (Kumar & Zhang, 2018; Aggarwal et al., 2007). It has been proposed that pharmaceuticals targeting molecular signaling pathways are the future for drug discovery as they target multiple protein targets and pathways (Kumar & Zhang, 2018; Aggarwal et al., 2007). The network connection results of this investigation are in line with this theory, as the pharmacological and mechanistic predictions provide evidence, albeit computational, of 155 strong affinity binding connections (≥−9.0 kcal/mol) between DCT-compounds and five proteins associated with multiple asthma KEGG disease pathways (Fig. 3).

Specific DCT-compound insights include the compounds that targeted all crystal complexes of KDR, LTA4H, PDE4D or PPARG (Table 3), as each of the crystal complexes for these proteins are in the same classification; agonist or antagonist. These results suggest that these compounds may promote or inhibit the biological actions associated with each protein through accumulation and synergism. Furthermore, the identified DCT-compounds that bound with all investigated crystal complexes of PDE4D and PPARG are significant justifications for the historical clinical evidence of DCT reducing the symptoms of bronchoconstriction and airway inflammation, as we know these are the outcomes of registered asthma pharmaceutical compounds theophylline (PDE4D antagonist), ibudilast (PDE4D antagonist) and rosiglitazone (PPARG agonist). Within these DCT-compounds there are 12 that bound to all PDE4D crystal complexes and all PPARG crystal complexes and are especially highlighted for further investigation as PDE4D antagonists or PPARG agonists. These included four from G. uralensis Fisch (Gan Cao) (apioside-Gc013, gancaonin F-Gc024, paratocarpin B-Gc088 and semilicoisoflavone B-Gc090), five from M. alba L. (Sang Bai Pi) (cyclomulberrin-Sb005, kuwanon S-Sb011, moracin E-Sb016, morusin-Sb020 and sanggenon N-Sb041) and three from S. baicalensis Georgi (Huang Qin) (chrysin-7-glucuronide-Hq020, dihydrobaicalin-Hq021 and scutellarin-Hq030).