Large-Scale Target Identification of Herbal Medicine Using a Reverse Docking Approach.

Herbal medicine has been used to countermine various diseases for centuries. However, most of the therapeutic targets underlying herbal therapy remain unclear, which largely slow down the novel drug discovery process from natural products. In this study, we developed a novel computational pipeline for assisting de novo identification of protein targets for herbal ingredients. The pipeline involves pharmacophore comparison and reverse ligand-protein docking simulation in a high throughput manner. We evaluated the pipeline using three traditional Chinese medicine ingredients such as acteoside, quercetin, and epigallocatechin gallate as examples. A majority of current known targets of these ingredients were successfully identified by the pipeline. Structural comparative analyses confirmed that the predicted ligand-target interactions used the same binding pockets and binding modes as those of known ligand-target interactions. Furthermore, we illustrated the mechanism of actions of the ingredients by constructing the pharmacological networks on the basis of the predicted target profiles. In summary, we proposed an efficient and economic option for large-scale target exploration in the herb study. This pipeline will be particularly valuable in aiding precise drug discovery and drug repurposing from natural products.

Introduction

Herbal medicines show a wide range of therapeutic effects in countermining cancers, virus infections, inflammations, hypertension, and so on.[1] The power of herbal therapeutics has attracted increasing attentions. Each year, an increasing number of publications described the mechanism of action (MOA) underlying herbs such as ginseng,[2] echinacea,[3] green tea,[4] and ginger.[5] Some herbal ingredients such as pancratistatin[6] and paclitaxel[7] have been applied in cancer therapy.

In addition to the conventional “wet-lab” experimental approaches, computer-aided drug discovery (CADD) serves as an attractive complementary option in the efficient and economic drug design.[8] For instance, docking-based methods have shown their impressive performance in large-scale structure-based virtual screening for candidate drugs.[9] Most of docking approaches screen pool of ligands against a defined protein cavity to identify the best fitting pose. On the other hand, the reverse docking methods search for putative protein targets that bind to a single ligand.[10] It has been acknowledged that biologically active herbal ingredients may exert therapeutic effects through a mode involving multiple targets and low dosage.[11] Hence, the reverse docking methodology meets the requirements of current mechanistic investigation toward herbal medicine. Previous studies have proposed reverse docking as a viable direction in the network exploration of herbal medicine[12] and in systematic understanding of drug pharmacology and toxicology.[13] Up-to-date, only a few tools such as INVDOCK,[12] TarFisDock,[14] and idTarget[15] are available for rapid identification of protein targets for a ligand. These tools are powerful while having several disadvantages: (1) They are either limited by small cavity database or constrained by the long computing time. (2) It is challenging to determine a proper threshold value to filter out those potential false ligand–receptor interactions. (3) They are inconvenient to be implemented or updated locally, which largely hinders their use in practical applications.

Here, we introduce a novel computational pipeline for large-scale identification of protein targets for a given ligand. We evaluated the pipeline efficiency using three well-studied herbal ingredients, namely, acteoside,[16] epigallocatechin gallate (EGCG),[17] and quercetin[18] as examples. Furthermore, we illustrate the possible molecular basis of MOA underlying these herbal ingredients on the basis of their predicted protein target profiles.

Results

Large-Scale Target Search of Herb Ingredients

Overall, the computational pipeline identifies 151, 143, and 128 nonredundant targets for aceteoside, quercetin, and EGCG, respectively. In this study, we further refined the target lists by the vina score threshold of −8.5 kcal/mol for later target analyses. The refined targets were most likely ingredient targets, and they were listed in a descending order of the estimated binding ability in Supporting Information Tables S1–S3. For aceteoside, 38 out of 151 predicted targets owned a vina score stronger than the threshold, including 10 cancer-related targets and 2 asthma-related targets (Supporting Information Table S1). Among the 38 targets, there are 4 known targets for aceteoside and 30 predicted targets. These results are in agreement with current applications of aceteoside in antiallergic (asthmatic) and anticancer therapy.[19,20] In the same way, we refined 20 and 19 high binding affinity targets for quercetin and EGCG, respectively. Some of these targets have been aimed at countermining various diseases such as cancer, Alzheimer disease, and some psychiatric disorders (Supporting Information Tables S1–S3). Among the 20 refined targets of quercetin, there are 7 known targets and 13 predicted targets. Among the 19 refined targets of EGCG, there are 7 known targets and 12 predicted targets.

Target Validation by Flexible Binding Analysis

Because the pipeline preset a large pocket size for high throughput docking simulation, it’s thus necessary to evaluate the capability of the Autodock Vina in locating the correct binding site. For this purpose, we chose the top five predicted protein targets in this study and downloaded the corresponding known ligand–target complexes from the PDB database, if available, as reference to evaluate the binding modes of the predicted ingredient–target complexes. These selected targets consisted of nine acteoside–target complexes, seven quercetin–target complexes, and eight EGCG–target complexes in the predicted list (Supporting Information Tables S1–S3). Of the total 24 predicted ingredient–target complexes, 16 cases (9/9 for acteoside, 6/7 for quercetin, and 5/8 for EGCG, respectively) had highly similar binding modes with the known ligand–target complexes (Table ).

Table 1

Statistic Summary for the Predicted Ingredient–Target Interactions for Acteoside, Quercetin, and EGCG Listed in Supporting Information Tables S1–S3a

name of ingredient	number of predicted targets with binding energy above −8.5 kcal/mol	ratio of known targets in the top five predicted targets	total protein targets in VINA docking simulation	ratio of same binding modes as the known PDD structures in selected five targets
acteoside	34	3/5	151	5/8
quercetin	20	5/5	143	6/7
EGCG	19	4/5	128	5/8

The comparison of protein pockets and docking modes were undertaken on the selected top five protein targets of each ingredient docking list in Supporting Information Tables S1–S3.

As illustrated in Figure , both the original ligands (blue color) and ingredients (red color) were well aligned in the same protein pockets.

Figure 1

Binding comparison between ingredients and known ligands with four predicted targets. The first column of (a–d) illustrates the binding mode similarity. The red one stands for the ingredient, and the blue one stands for known ligand in the original PDB file. The second column of (a–d) illustrates the key residues in the ligand–protein complexes determined by the commercial software Molecular Operating Environment (MOE). The third column of (a–d) shows the key residues in the ingredient–protein complexes.

For instance, quercetin has very similar binding pose with SGFR and cGMP-specific 3,5-cyclic phosphodiesterase (PDE5A) as that of adenosine-5′-diphosphate and 5,7-dihydroxy-2-(4-methoxyphenyl)-8-(3-methylbutyl)-4-oxo-4H-chromen-3-yl 6-deoxy-α-l-mannopyranoside (7CA) (Figure b,d). However, in some cases, the best predicted binding site is not at the known ligand-binding pocket in the PDB structures. For example, the top prediction of quercetin with serine/threonine-protein kinase (Chk1) is not in the known pocket; however, the binding affinity between quercetin and the known pocket is comparable with that of between quercetin and the new pocket. This finding suggests Chk1 may have multiple binding sites, and the quercetin–Chk1 interaction is non-specific.

We have settled the parameters to generate maximum 20 output conformations for each docking, to increase the possibility to sample the correct pocket during blind docking. We have checked whether the docked structure was able to reproduce the X-ray conformation by checking the root-mean-square deviation (rmsd) between top score conformation and native conformation from docking of 64 proteins and its corresponding cocrystallized ligands. There are 23 in 64 cases that the rmsd between top score conformation and native conformation are smaller than 0.6 nm, and 38 in 64 cases that the rmsd between top score conformation and native conformation are smaller than 0.8 nm. It can see the ability of Autodock Vina to find the correct binding site, and conformation in the top score conformation is still challenging, but it is still one of the best method to find the binding conformation efficiently based on the best of my knowledge. Because the docking is a large approximation and the diversity of the vina score training data set is limited, it is still necessary to further validate the result with molecular dynamic (MD) simulation or wet experiment. Because the vina score cannot guarantee that the predicted target is a true target, the MD simulation, which is physical based and more rigid, is important to further validate before doing experimental validation.

Target Refinement Using MD Simulation

The MD simulation is helpful to determine the best binding conformation for the predicted ingredient–target interactions. In almost all cases, the predicted ingredient–target binding remained relatively stable after 100 ns MD simulation (Supporting Information Figure S1). Further replicated MD simulations consolidated the robustness of the simulation outcome (Supporting Information Figure S2). The protein target did not experience large conformation changes after two independent MD simulation, even they started from random initial velocities. We determined the binding free energy after MD for the predicted ingredient–target complexes in the Table using the tool G_mmpmsa.

Table 2

Binding Free Energy of NOS2 (PDBID 4nos) with Acteoside, PDE5A (PDBID 2h44) with Quercetin, and WARS (PDBID 1r6t) with EGCG, Respectivelya

complex name	binding energy (deviation) (kJ/mol)	van der Waal energy (deviation) (kJ/mol)	electrostatic energy (deviation) (kJ/mol)	polar solvation energy (deviation) (kJ/mol)	SASA energy (deviation) (kJ/mol)
acteoside–NOS2	–264.145 (2.185)	–352.715 (1.525)	–205.055 (2.686)	323.767 (1.347)	–30.148 (0.092)
quercetin–PDE5A	–185.516 (3.824)	–134.277 (3.057)	–199.478 (4.765)	163.389 (3.575)	–15.027 (0.295)
EGCG–WARS	–55.807 (4.840)	–92.895 (7.816)	–67.844 (6.139)	117.273 (11.493)	–11.394 (0.958)

The free energy was determined at the last 10 ns simulation trajectory by g_mmpbsa tools with inner dielectric constant 2 and solvent dielectric 80. The energy deviations are shown in the brackets.

Comparatively, the known protein targets were selected in priority (Table ). At the same time, we summarized the critical residues that contribute to free energy in Supporting Information Table S4. The snapshots after 100 ns MD simulation for the three test cases are shown in Figure . The MD simulation analysis provides an ensemble of conformations for the detailed interactions between the ingredients and target proteins with more accurate free binding energy. This is helpful for further precise drug redesign.

Figure 2

Snapshot conformations after 100 ns simulation for acteoside–NOS2 complex, and quercetin–PDE5A complex, and EGCG–WARS complex. Panel (a) shows the conformation of acteoside–NOS2 complex; panel (b) shows the conformation of quercetin–PDE5A complex; panel (c) shows the conformation of EGCG–WARS complex. The residues which have free energy contribution larger than 3 kJ/mol (blue) or smaller than −3 kJ/mol (red) were explicitly shown. The ligands were shown in green color.

Understanding Pharmacology of Herb Ingredients via Target Network

Targets of Acteoside and Their Pharmacological Activities

Of total 151 targets predicted for acteoside, 4 are known targets and 34 are highly potential targets. For instance, the NOS2 (or iNOS) is a known target of several approved or experimental drugs indications such as inflammatory disease, endocrine disorder, and antithrombotic therapy as documented in the DrugBank.[21] For instance, we predicted acteoside–NOS2 interaction with an Autodock Vina score of −10.3 kcal/mol. The binding position and binding mode were similar to that of the known protoporphyrin–NOS2 interaction (PDB_ID: 4nos) (Figure a). The acteoside–NOS2 interaction was also supported by previous work that acteoside could inhibit enzymatic activity of NOS2 and 5-lipoxygenase.[22] The cAMP-specific 3,5-cyclic phosphodiesterase 4B (PDE4) is an known target of acetoside’s analog plantamajoside.[23] It ranked top 2 in the predicted targets of acetoside. The aldose reductase (AKR1B1) is another known therapeutic target, which is related to the acteoside’s anti-diabetic activity.[24] In the predicted target list, the inosine-5-monophosphate dehydrogenase 2 (IMPDH2) owned a high binding affinity with acteoside (vina score: −9.5 kcal/mol). The IMPDH2 is a known target of immunosuppressive agents such as the approved drug mycophenolate mofetil and mycophenolic acid.[25] The predicted acteoside–IMPDH2 interaction may explain the immunomodulatory activity of acteoside in vitro.[26]

Via mapping the 151 predicted targets to the KEGG pathway by the DAVID webserver,[27] we explored the possible synergistic effects of acteoside (Supporting Information Table S5). Of the 151 predicted targets, 53 are related to metabolism, 20 are related to cancer, and 18 are related to the PI3K-Akt signaling pathway. In particular, 10 predicted targets are involved in the prostate cancer (path: hsa05215), including GSK3B and CDK2 that own predicted affinity better than the cutoff −8.5 kcal/mol. These results partially help to understand the cancer cell cycle arresting effects of acteoside.[28] Seventeen predicted targets are related to the biosynthesis of antibiotics, and seven are involved in the measles pathway (path: hsa05162). Interestingly, it is noted the Tibetan herb “Ye-Xin-Ba” which contains acteoside was used to treat measles.[29] Besides, four potential targets are related to the Alzheimer’s disease, which may be responsible for the memory enhancement effect of acteoside.[30] The possible pharmacology network is shown in Figure .

Figure 3

Pharmacology network of acteoside by mapping the potential targets to the KEGG pathways.

During construction of the pharmacological networks, we highlight targets that have predicted binding affinity stronger than the threshold for aceteoside because we care more about the pharmacological effects of the more relevant targets: targets of quercetin and their pharmacological activities. The pipeline predicted 143 targets for quercetin, including 7 known therapeutic targets and 13 highly potential targets. The cGMP-specific 3,5-cyclic phosphodiesterase (PDE5A, PDBID: 2h44) is on the top of potential targets in Supporting Information Table S2. This result is consistent with the previous finding that extraction of red onion peel could inhibit phosphodiesterase 5A (PDE5A), and the activity was mainly caused by flavonoid quercetin.[31] Because PDE5A is strongly related to erectile dysfunction,[31] quercetin may serve as the potential lead in countermining male sexual dysfunction. We also predicted three known therapeutic targets for cancer therapy: the serine/threonine-protein kinase Chk1 whose activity could be decreased by quercetin,[32] serine/threonine-protein kinase pim-1 (PIM1) whose activity was inhibited by quercetin (IC50: 1.1 μm)[33] and proto-oncogene tyrosine-protein kinase Src (SRC) whose activity was inhibited by (64.0 ± 4.24)% using 50 μm quercetin.[34,35] Besides, quercetin owned the phytoestrogen-like activity, which is likely via the interaction with the estrogen receptor β.[36,37]

We mapped the 143 predicted targets into the KEGG pathway database. Many targets are related to the metabolic pathways (43 targets), cancer pathways (23 targets), PI3K-Akt signaling pathway (21 targets), Ras signaling pathway (20 targets), and so on. In particular, 13 protein targets are involved in the tuberculosis (path: hsa05152), including 2 high-affinity binding proteins SRC and JAK2. This finding may partially explain the quercetin’s ability in combating tuberculosis.[38] The possible pharmacology network is shown in Supporting Information Figure S6.

Targets of EGCG and Their Pharmacological Activities

Overall, the pipeline predicted 128 targets for EGCG, including 6 known therapeutic targets and 11 highly potential targets. Of these targets, the β-secretase 1 (BACE1), a Alzheimer disease-related protein, owned high affinity with EGCG (vina score:–9.9 kcal/mol).[39] This result agreed with prior findings that green tea catechins was an inhibitor of BACE1 and may partially explain the potential application of catechins in Alzheimer’s disease treatment.[40] Moreover, the proto-oncogene serine/threonine-protein kinase pim-1 (PIM1),[41] histamine N-methyltransferase, peroxisome proliferator-activated receptor alpha (PPAR-alpha), macrophage metalloelastase (MMP12), and glutathione S-transferase A3 (GSTA3) are all known targets of EGCG.

The tryptophanyl-tRNA synthetase, cytoplasmic (WARS, PDBID: 1r6t) is on the top of predicted targets in Supporting Information Table S7. Previous study found the WARS expression downregulated by EGCG.[42] Literature surveillance also suggested that WARS was likely the target of l-tryptophan,[43] a drug used in promoting serotonin production, helping healthy sleep, relieving depression, strengthen pain tolerance, and controlling weight.[44] Interestingly, green tea can promote healthy sleep, as well as mental and emotional well-being.[45] EGCG, as a major activity ingredient of green tea, may contribute a lot to green tea’s antifatigue effect.[46]

We mapped the 128 predicted targets of EGCG into the KEGG pathway database. Many targets are related to the metabolic pathways (44 targets), cancer pathways (19 targets), PI3K-Akt signaling pathway (18 targets), and so on. The possible pharmacology network is illustrated in Supporting Information Figure S4.

Conclusion

Thousands of years’ practice of herbal medicine has provided valuable clues for safe and better drug development. However, there still is a huge demand of information such as herb ingredient–target interactions for aiding rational drug discovery. The docking method is an economical and practical way to potentially understand MOA of herb by simulation of ligand–receptor binding. In this study, we proposed a novel computational pipeline for high throughput ligand target search against the user-defined structure database. Comparing with the current existing reverse docking packages, the new pipeline has several advantages: (1) the pipeline is flexible in controlling scope of target search by self-defining the protein structure database. This keeps the pipeline unrestrained in dealing with various simulation jobs without largely increasing maintenance cost. (2) The pipeline is comparatively more efficient. Combination of pharmacophore comparison and flexible docking allows large-scale target screening to be done in a reasonable computational cost. (3) The core components of the pipeline are two state-of-art tools, PharmMapper and Autodock Vina, that can be easily accessed on-line or downloaded. In case the PharmMapper fails, we further introduced an alternative target search method on the basis of ligand similarity, whose core is a freeware LiSiCA, to help narrowing down the target pool. Hence, the pipeline can be easily updated accompanying with the component software update. (4) The optional MD simulation on selected docking results further refines the predicted targets, thus largely reducing the false positives.

Using this pipeline, we predicted the potential targets of three herbal ingredients, acteoside, quercetin, and EGCG, in human structural proteome. Upon the predicted targets, we explored the possible MOAs underlying these three ingredients by demonstrating KEGG pathway enrichment analyses. Furthermore, we built the molecular interaction network for each ingredient to aid better understanding of pharmacological activities.

However, the computational pipeline is still far away from perfectness. For instance, the pipeline is constrained in search of protein targets with well-solved structures. The pipeline highly relies on the performance of PharmMapper. Its background pharmacophore database and online interactive job submission have become the bottleneck of the pipeline. Although, we introduced a ligand-based target search method as an alternative choice to PharmMapper, the method is also limited by the availability of ligand–protein complexes, for example, the complete ligand–protein interactome extracted from the PDB. Hopefully, with the dramatical increase in 3D structures of ligand–protein interactions every year and the consistent enrichment of PharmMapper background database,[47] this weakness is expected to be countermined in the future.

Instead of exploring target network of isolated compounds, it would be more interesting to research several compounds within an herb or from the same source of plant. This is because several ingredients from one herbal may exert a synergic effect or adding effect by acting on multiple pathways. In the future, we can search protein targets with our current pipeline for several herbal ingredients from one herb to explore the underlying herb curing mechanisms in a more efficient and systematic way.

Nevertheless, the reverse docking strategy introduced in this study provides an alternative and practical way to rapidly identify the possible protein targets within acceptable computational time. It would be useful in aiding current pharmacological understanding of novel compounds in a systematic manner. Furthermore, it will prompt the modernization process of herbal medicine by identifying potential ingredient–target interactions.

Materials and Methods

Computational Pipeline for Reverse Docking

The computational pipeline consists of two major steps (Figure ): the pharmacophore-based target filtering and reverse docking simulation.

Figure 4

Reverse docking pipeline proposed in this study. The computational pipeline consists of three major stages. In the first stage, the PharmMapper sever is used to narrow down the protein pool by the pharmacophore matching. In the second stage, reverse docking by Autodock Vina is used for selecting target candidates by structure-based virtual screening simulation. Finally, MD simulation is used to refine and validate the candidate targets.

Generation of Potential Target Pool by PharmMapper

The target filter step helps to filter out the unlikely targets from the original protein structure pool by pharmacophore comparison. We used an on-line tool PharmMapper (http://lilab.ecust.edu.cn/pharmmapper/index.php)[48] for rapid prediction of potential targets by comparing the query ligand against the predefined pharmacophore models. Following the guidelines of PharmMapper, we uploaded the chemical structures via the job submission form of the website. The default running parameters were adopted. The server outputs a list of potential ligand targets subject to the top mapping scores. In this study, we selected the top 300 PDB structures of human proteins as potential target pool for further prediction.

As an alternative solution to PharmMapper, we have developed a tool that can be used for inverse protein target searching based on deep neural network model and word2vec techniques.[49] The protein pocket and ligand was converted into vectors as input, and a densely fully connected neural network was used to do the training over the PDBBIND database. The resulted model (we named it DeepBindVec here) shows high efficiency and relative high accuracy (AUC above 0.9 in test set) and indicates that it’s suitable for inverse protein target searching and has potential to replace PharmMapper server in our pipeline. It is independent of protein–ligand binding complex structure; thus, it has the potential to overcome the limitation of LigandScout[50] which is the core of PharmMapper. When using the DeepBindVec as core to do inverse target searching, we named it IVS2vec, deposited in the gitlab (https://gitlab.com/llb671205/densely-fully-connected-neural-network/tree/master). We also proposed a ligand-based target search method to help narrowing down the target pool as another alternative solution to PharmMapper. The principle of the method is that ligands having similar structures may have similar protein binding abilities. The core components of the ligand-based target search method are a freeware LiSiCA[51] (http://insilab.org/lisica/) and user-defined ligand–protein interaction database. The LiSiCA applies the clique algorithm to measure the 2D and 3D similarities between query ligand and target ligands by the Tanimoto coefficients in large scale. The background ligand–protein database can be created by extracting user-interested ligand/cofactor–protein interactions from the protein data bank (PDB), which includes about total 14 761 ligand–protein complexes. Alternatively, the user can download the prebuilt ligand–protein interaction data set, the PROTEUS data set, from insilab server (http://insilab.org/datasets/) for convenience. The PROTEUS data set was created upon the entire PDB data. Underlying the ligand-based target search is a database of known ligand–protein interactions. We converted the ligand from the PDB format into the MOL2 format using the openbabel.[52] Then, we wrote a python script to remove the redundant ligands from the file. By doing so, we built the nonredundant relation database between ligand and protein targets. For user convenience, we deposit the demo package of the method at the GitHub page (https://github.com/haiping1010/LBTS). To be noted, several previous studies utilized the similar ligand-based target search method in prediction of potential drugs for PDE4 and PDE5.[53] In this study, we chose the PharmMapper approach, instead of ligand-based target search method, as the component of the reverse docking pipeline for better illustration.

Target Refinement by Reverse Docking Using Autodock Vina

We used Autodock Vina[54] to simulate ligand–target interaction in large scale. The simulation starts by selecting the receptor mass center as the pocket center using the python module pdb_centermass.py of pdbTools_0.2.1.[55] We set a large cavity volume space (126 Å in x, y, z dimension) as the pocket space. The simulation was initialized by converting the input PDB file format to PDBQT file format using the AutoDockTools.[56] The exhaustiveness was set to 8, the num_modes was set to 9, and the energy_range was set to 3. We adopted the default scoring function and optimization algorithm of Autodock Vina, which have been well discussed in a previous article.[54] The simulation eventually outputted a list of ligand–protein interactions in a descending order of binding affinity. In this study, we selected the most likely targets for further validation by setting a binding energy threshold value of −8.5 kcal/mol. The choice of the proper threshold value is a tradeoff between true positive rate and precisions. In our case, we only focus on the potential targets that can bind with the acquired ligand with relatively high affinity. If the cutoff is very low (such as −6 kcal/mol), then the precisions would be very low, many nonbinding protein targets will be chosen incorrectly. While if the cutoff is very high (such as −10 kcal/mol), the true positive rate would be low, and many binding protein targets will not be identified.

It should be noted that the docking method can be used to do the reverse-docking on all available PDB structures without the import from PharmMapper. There are more than 30 000 human proteins deposited in the PDB database. It would take about 3–4 days to finish all of the blind docking with ∼60 CPU cores. If we only focus on the disease-related therapeutic target (about 2000–3000), it would be completed within 1 day by this reverse docking with 10–20 CPU cores. Also, if we focus on the adverse drug reaction (about several hundred), the time that is needed to finish the reverse docking would be shorter. Thus, the pipeline can be flexible without the importing from pharmMapper.

Prediction of Protein Targets for Herb Ingredients Using the Computational Pipeline

Acquisition of Herb Ingredients

In this study, we used three herb ingredients acteoside, quercetin, and EGCG as examples for pipeline evaluation. Acteoside is a bioactive ingredient of traditional chinese medicine (TCM) drug Herba Cistanches. It also often exists in the species of the Lamiales order.[57] Acteoside has several pharmacological properties including being antioxidant, immunosuppressive, and anti-inflammatory.[58] Quercetin is a widely distributed flavonoid in nature; its potential applications include but are not limited to anticancer, anti-inflammatory, and anti-asthmatic.[17] EGCG is a type of catechin.[59] EGCG is a major active component in green tea and has potential therapeutic activity in many disorders such as cancer and chronic fatigue syndrome.[18] We retrieved structures of these three herbal ingredients from the PubChem database.[60] Before further use, we converted the downloaded *.sdf files to *.pdb format by the UCSF Chimera.[61]

Prediction of Ingredient Targets

We submitted the three herbal ingredients to the PharmMapper server for PDB wide target prediction, achieving a list of 300 potential targets for each ingredient. We downloaded protein structures from the PDB database. We cleaned the structures by removal of cocrystalized ligands and addition of hydrogen bonds using the VMD[62] or the UCSF chimera[61] via TCL or python scripts, respectively. Then, we simulated the ingredient–receptor binding by docking each ingredient against the predefined pockets of the potential protein pool using Autodock Vina. The reverse docking process took about 24 h of computational time under a Linux system with seven CPU cores (1064 Hz/core) for each ingredient. The outcome of docking simulation was a list of PDB IDs, sorted in descending order of binding affinity.

Refinement of Ingredient Targets by Classical MDs

To further refine the ingredient–target interactions for accurate complex conformation, we carried out MD simulation using the Gromacs program with AMBER-14SB force field.[63,64] In this study, we selected three binding complexes for MD simulation as examples: NOS2 (PDBID 4nos) with acteoside, PDE5A (PDBID 2h44) with quercetin, and WARS (PDBID 1r6t) with EGCG. The starting protein structures were downloaded from the PDB structure database.[65] The topology of ligand was generated by ACPYPE,[66] which relies on Antechamber.[67] The partial charges of ligand were calculated with R.E.D development sever using the Gaussian method.[68,69] First, we created a cubic box and put the ligand–receptor complex at the center. A minimum distance from the protein to box edge was set to 1 nm. We filled the cubic box with TIP3P water molecules[70] (the numbers of water molecules are 42 271, 18 513, 27 676, respectively) and then added counter ions (the numbers of sodium ions are 1, 4, 3, respectively) to neutralize the total charge using the Gromacs program tool.[71] Eventually, the three systems contained 140 399, 60 885 and 89 260 total atoms, respectively. We used the particle mesh Ewald method[72] for calculation of electrostatic interactions under the periodic boundary conditions. A cutoff of 14 Å was set in quantifying van der Waals nonbonded interactions. Covalent bonds involving hydrogen atoms were constrained by applying the LINCS algorithm.[73] We performed the MD simulation with parameters of 50 000 energy minimization steps with a step-size of 0.01 nm, 100 ps simulation with isothermal–isovolumetric ensemble (NVT), and 100 ps simulation with an isothermal–isobaric ensemble (NPT) for water equilibrium. This was followed by a 100 ns NPT production run (step size 2 fs) with initial velocities from the previous NPT equilibrium step. The Parrinello–Rahman barostat and modified Berendsen thermostat were used for simulation with a fixed temperature of 308 K and a pressure of 1 atm. For each of the three complexes, we demonstrated two replicated simulations with random initial velocities to check the robustness of MD. We measured the rmsd of final protein conformations of two MD replicates. A small rmsd value indicates the robustness of the MD simulation. Although MD simulation is expensive compared with docking, it is still much cheaper than experimental validation. The MD simulation can also be conveniently carried out on this academic free webserver (http://mmb.irbbarcelona.org/MDWeb/).[74]