Identification of potent EGFR inhibitors from TCM Database@Taiwan.

Overexpression of epidermal growth factor receptor (EGFR) has been associated with cancer. Targeted inhibition of the EGFR pathway has been shown to limit proliferation of cancerous cells. Hence, we employed Traditional Chinese Medicine Database (TCM Database@Taiwan) (http://tcm.cmu.edu.tw) to identify potential EGFR inhibitor. Multiple Linear Regression (MLR), Support Vector Machine (SVM), Comparative Molecular Field Analysis (CoMFA), and Comparative Molecular Similarities Indices Analysis (CoMSIA) models were generated using a training set of EGFR ligands of known inhibitory activities. The top four TCM candidates based on DockScore were 2-O-caffeoyl tartaric acid, Emitine, Rosmaricine, and 2-O-feruloyl tartaric acid, and all had higher binding affinities than the control Iressa®. The TCM candidates had interactions with Asp855, Lys716, and Lys728, all which are residues of the protein kinase binding site. Validated MLR (r² = 0.7858) and SVM (r² = 0.8754) models predicted good bioactivity for the TCM candidates. In addition, the TCM candidates contoured well to the 3D-Quantitative Structure-Activity Relationship (3D-QSAR) map derived from the CoMFA (q² = 0.721, r² = 0.986) and CoMSIA (q² = 0.662, r² = 0.988) models. The steric field, hydrophobic field, and H-bond of the 3D-QSAR map were well matched by each TCM candidate. Molecular docking indicated that all TCM candidates formed H-bonds within the EGFR protein kinase domain. Based on the different structures, H-bonds were formed at either Asp855 or Lys716/Lys728. The compounds remained stable throughout molecular dynamics (MD) simulation. Based on the results of this study, 2-O-caffeoyl tartaric acid, Emitine, Rosmaricine, and 2-O-feruloyl tartaric acid are suggested to be potential EGFR inhibitors.

Introduction

Target-specific therapies have generated much attention in addition to conventional cancer treatments [1]–[3]. By targeting key molecules essential for cellular function, replication, or tumorigenesis, such therapies may exert cytostatic or cytotoxic effects on tumors while minimizing nonspecific toxicities associated with chemotherapy or irradiation [4].

The epidermal growth factor receptor (EGFR) signaling pathway is one of the most important pathways in mammalian cells [5]. Specific ligands, such as epidermal growth factor (EGF) and transforming growth factor alpha (TGFα), bind and activate EGFR, triggering autophosphorylation of the intracytoplasmic EGFR tyrosine kinase domain [6], [7]. The phosphorylated tyrosine kinase residues serve as binding sites for signal transducers and activators of intracellular substrates, which then stimulate intracellular signal transduction cascades that upregulate biological processes such as gene expression, proliferation, angiogenesis, and inhibition of apoptosis [8]. EGFR overexpression has been shown to activate downstream signaling pathways, resulting in cells that have aggressive growth and invasive characteristics [9]. Tumor cell motility, adhesion, metastasis, and angiogenesis have also been associated with stimulated EGFR pathways [10]–[12]. Since EGFR over-expression often differentiates tumor cells from normal cells, it is possible for EGFR inhibitory molecules to act on tumor cells and attenuate their proliferation rates [4].

Several tyrosine kinase inhibitors were approved for clinical use. Iressa® (gefitinib) is highly selective for EGFR tyrosine kinase and is commonly used for treating lung cancer [13]. EGFR downstream signaling is competitively inhibited by Iressa® at its ATP binding site [14]. Other therapeutic agents with inhibitory mechanisms similar to Iressa® include Erlotinib (Tarceva®) against non-small cell lung cancer (NSCLC) and pancreatic cancer [15], [16], and Vandetanib (Zactima®) against late stage medullary thyroid cancer [17]. Lapatinib (Tykerb®) is a dual inhibitor of EGFR and HER2 tyrosine kinases approved for metastatic breast cancer [18], [19]. Though the effect of Iressa® on lung cancer has been well established, severe side effects has also been reported [20]. Adverse reactions listed under Iressa® product information include diarrhea, skin rash and dryness, nausea, vomiting, haemorrhage, anorexia, asthenia, and in some cases, interstitial lung disease with fatal outcomes [21]. The adverse effects of available treatments necessitate continuous search efforts for alternatives with less toxicity.

Computational predictions in biology and biomedicine are of significant importance for generating useful data which otherwise be time-consuming and costly through experiments alone [3], [22]–[27]. Computational predictions, combined with information derived from structural bioinformatics analysis, can provide useful insights and timely information for both basic research and drug development [28], [29]. Much cutting-edge cancer drug development has been conducted through the use of computational bioinformatics and modeling [30]–[37]. The powerful ability of modern computational prediction and bioinformatics were adopted in this research to search for novel EGFR inhibitors.

Traditional Chinese medicines (TCM) are natural substances with therapeutic effects on many diseases [38]–[40]. The vast number of TCM represents a rich resource that can be explored for drug development. We had investigated kinase inhibitor candidates from TCM targeting HER2 and HSP90 receptors before [28], [41]–[42]. Though EGFR kinase inhibitors have been investigated through different screening and modeling scenarios [43]–[47], none from TCM compounds has been reported to date. This study utilized the world’s largest TCM Database@Taiwan [48] to screen for potential EGFR inhibitors from TCM compounds and applied structure- and ligand-based methods to evaluate the suitability of candidates as EGFR inhibitors.

Materials and Methods

A useful predictor for a biological system should include the following steps [49]: (i) selection of a valid dataset to train and test the predictor; (ii) formulate samples with an effective mathematical expression that reflects intrinsic correlation with the attribute to be predicted; (iii) develop a powerful algorithm to operate the prediction; (iv) objectively evaluate accuracy of the predictor through cross-validation tests. The experimental design of the current study follows these guidelines and details are explained in the following sections.

EGFR Protein Sequence, Structure, and Characteristics

The EGFR protein sequence (EGFR_HUMAN, P00533) used in this study was obtained from Swiss-Prot [50], and the 3D structure (PDB: 2ITY) [51] used for analyses was downloaded from Protein Data Bank. The tyrosine kinase was encoded from Phe712-Leu979, and the ATP binding site was located at Leu718–Val726.

Candidate Screening and Docking Studies

The Traditional Chinese Medicine (TCM Database@Taiwan) database (http://tcm.cmu.edu.tw) was used to screen for potential EGFR inhibitors from more than 20,000 compounds within the database. All compounds were operated using the Prepare Ligands module with Lipinski’s rule of five using Discovery Studio 2.5 (DS 2.5; Accelrys Inc., San Diego, CA). Iressa® was selected as the control. The LigandFit program (DS 2.5) was used to locate the best docking pose for different confirmations under the Chemistry at HARvard Macromolecular Mechanics (CHARMm) force field [52]. Results for the docking studies were ranked according to Dock Score.

Descriptor Generation Using Genetic Function Approximation (GFA) Algorithm

Twenty ligands with demonstrated inhibition against EGFR were used in this study (Table S1) [53]. Descriptors for each ligand were identified using the Calculate Molecular Properties program in DS 2.5. Predictive models containing five optimum descriptors were generated using the Genetic Function Approximation (GFA) algorithm. Descriptors in the model with the highest r2 value were used to construct ligand activity prediction models.

Ligand Activity Predictions Using Multiple Linear Regression (MLR) and Support Vector Machines (SVM)

A MLR model using the descriptors from the top GFA algorithm was constructed using Matlab Statistics Toolbox (MathWorks, Natick, MA) and validated using MLR Leave-One-Out validation [54]. The MLR model was trained with 17 randomly selected ligands with EGFR inhibitory activity (Table S1) and used to predict the activity (pIC50) of the control and TCM candidates.

The identical descriptors were normalized to the range of [−1,+1] and plugged into the libSVM program to generate a SVM prediction model[55]. Following model training with the 17-ligand training set, the SVM model was used to predict the activity of the control and TCM candidates.

3D Quantitative Structure-Activity Relationship (QSAR) Model

Ligands used in the previous sections were also used for 3D-QSAR analysis. The 2-dimensional (2D) and 3-dimensional (3D) ligand structures were drawn with ChemBioOffice 2008 (PerkinElmer Inc., Cambridge, MA) under a Molecular Mechanics 2 (MM2) force field. Following ligand alignment, Comparative Molecular Field Analysis (CoMFA) and Comparative Molecular Similarities Indices Analysis (CoMSIA) models were constructed using partial least squares statistical method (PLS). Cross-Validated (CV) correlation coefficient (q2) and non-cross-validation correlation coefficient (r2) were used to validate the models. Biological activities of Iressa® and TCM candidate compounds were predicted using the generated 3D-QSAR contour map.

Molecular Dynamics Simulation

Molecular dynamics (MD) of Iressa® and the TCM candidates were simulated using DS2.5 Standard Dynamics Cascade and Dynamics package. Sample preparation was conducted under the following parameters: [minimization] steepest descent and conjugate gradient: each with maximum steps of 500, [heating time] 50 ps, [equilibration time] 200 ps. The simulations were produced with a total production time of 20 ns with NVT, constant temperature dynamics of Berendsen weak coupling method, a temperature decay time of 0.4 ps, and a target temperature of 310K. Root mean square deviations (RMSD) of protein-ligand complex and individual ligands, total energy of protein-ligand complex, hydrogen bond (H-bond), and H-bond distance were analyzed using the Analyze Trajectory function following MD simulation.

Results/Discussion

The top four TCM candidates with the highest Dock Score were 2-O-caffeoyl tartaric acid, Emitine, Rosmaricine, and 2-O-feruloyl tartaric acid (Table 1). Corresponding scaffolds of the top TCM candidates are illustrated in Figure 1. Iressa®, Emitine, and Rosmaricine had amine groups available for H bonding whereas 2-O-Caffeoyl tartaric acid and 2-O-feruloyl tartaric acid had carbonyl groups. The different residues available for H bonding resulted in different binding poses (Figure 2). Binding of Iressa® (Figure 2a), Emitine (Figure 2c), and Rosmaricine (Figure 2e) to tyrosine kinase were located within the pocket, with H-bonds formed between the amine group of the ligand compounds and the carboxyl group of Asp855. 2-O-Caffeoyl tartaric acid (Figure 2b) and 2-O-feruloyl tartaric acid (Figure 2e) docked outside the tyrosine kinase pocket and formed multiple H-bonds through their carboxyl groups with Lys716 and Lys728. The binding location of 2-O-caffeoyl tartaric acid and 2-O-feruloyl tartaric acid directly blocks the ATP binding site of tyrosine kinase located from Leu718–Val726. Dock scores of each TCM candidate is given in Table 1. All candidates have higher dock scores than Iressa®, indicating higher binding affinities to the tyrosine kinase receptor than Iressa®.

Table 1

Docking score and biological activity predictions of top TCM candidates in comparison with the control.

		Predicted Activity (pIC50)
Compounds	Dock Score	MLR	SVM
2-O-Caffeoyl tartaric acid	118.065	8.386	7.041
Emetine	113.424	5.548	7.140
Rosmaricine	105.149	4.543	7.463
2-O-Feruloyl tartaric acid	104.284	8.359	7.242
Iressa*	60.791	6.715	5.110

*control.

Figure 1

Scaffolds of the control and TCM compounds.

(A) Iressa, (B) 2-O-Caffeoyl tartaric acid, (C) Emetine, (D) Rosmaricine, and (E) 2-O-Feruloyl tartaric acid.

Figure 2

Docking pose of different compounds in EGFR using LigandFit.

(A) Iressa, (B) 2-O-Caffeoyl tartaric acid, (C) Emetine, (D) Rosmaricine, and (E) 2-O-Feruloyl tartaric acid. Binding site amino acids are shown in red (negatively charged amino acids) and blue (positively charged amino acids). Hydrogen bonds are color coded based on bond distance: 1.6–1.8Å (orange), 1.8–2.0Å (yellow), 2.0–2.2Å (light green), 2.2–2.4Å (neon green), and 2.4–2.5Å (dark green).

Ligand Activity Predictions Using MLR and SVM

Representative descriptors from the top GFA algorithm include: Num_H_Acceptors_Lipinski (equivalent of N+O count), Molecular_SurfaceArea (the total surface area for each molecule using a 2D approximation), Kappa_1 (Kappa Shape Indices), PMI_Y (principle moment of inertia Y-component), and Shadow Xlength (length of molecule in the X dimension). The descriptors were validated using Leave-One-Out method which is the most objective of all available cross-validation methods [56]. The MLR model established with the determined descriptors was:

The SVM model was also established with the five identified descriptors using libSVM.

Correlation between the predicted and observed pIC50 activities on EGFR ligands of known activity using the constructed MLR and SVM models were illustrated in Figure 3a and 3b, respectively. Correlation coefficients based on the training set were 0.7858 for the MLR model and 0.8754 for the SVM model. Activity predictions of Iressa® and the TCM candidates using MLR and SVM were summarized in Table 1. Both models indicate that Iressa and the TCM candidates are compounds with acceptable predicted activities. Predicted activities (pIC50) of Iressa by the trained MLR and SVM models were 6.715 and 5.110, respectively. The Iressa activity predicted by SVM was closer to experimentally determined Iressa activities (pIC50) between 4.76–5.96 [57], thus SVM values may be more accurate predictions of the actual activity.

Figure 3

Correlation of observed and predicted activity (pIC50) using 2D-QSAR models.

(A) MLR and (B) SVM.

3D QSAR Model

The results of CoMFA and CoMSIA model generation are detailed in Table 2. Steric field was the sole factor in the CoMFA model since the electrostatic field value was zero. Cross-validated (q2) and non-cross-validated (r2) correlation coefficient values of 0.721 and 0.986, respectively, indicated a high level of confidence for this model. The small standard error of estimates (SEE) and large F-test value further supported the reliability of this model. In contrast, CoMSIA models were influenced by multiple factors including steric field, hydrophobic region, and hydrogen bond donor/acceptors. Among all generated versions of the CoMSIA model, CoMSIA_SHD had the highest r2 (0.988), lowest SEE (0.133), and highest F value (134.272), thus was selected as the optimum CoMSIA model for use in this study. The pIC50 of 20 ligands predicted by the constructed CoMFA and CoMSIA models were compared with observed pIC50 reported by Fidanze et al. [53] in Table 3. In general, both models gave similar predicted values and were close to the experimentally determined activities. Correlations between predicted and observed pIC50 using CoMFA and CoMSIA models are summarized in Figure 4a and 4b, respectively. High correlation coefficients validated the reliability of the constructed CoMFA (r2 = 0.9860) and CoMSIA(r2 = 0.9877) models.

Table 2

CoMFA and CoMSIA models as a factor of various fractions and the corresponding validation values.

	Cross validation		Non-cross validation			Fraction
	ONC	q2 cv	r2	SEE	F	S	H	D	A
CoMFA	6	0.721	0.986	0.142	117.843	1.00	-	-	-
CoMFA
S	6	0.764	0.975	0.189	65.257	1.00	-	-	-
H	6	0.331	0.980	0.168	83.249	-	1.00	-	-
D	6	0.236	0.945	0.281	28.780	-	-	1.00
A	6	−0.344	0.784	0.558	6.034	-	-	-	1.00
SH	6	0.541	0.986	0.141	118.218	0.37	0.63	-	-
SD	6	0.615	0.982	0.159	92.965	0.48	-	0.52	-
SA	6	0.718	0.977	0.180	72.008	0.53	-	-	0.47
HD	6	0.569	0.982	0.160	91.804	-	0.56	0.44	-
HA	6	0.265	0.979	0.176	76.040	-	0.60	-	0.40
DA	6	0.003	0.956	0.250	36.635	-	-	0.60	0.40
SHD*	6	0.662	0.988	0.133	134.272	0.26	0.40	0.35	-
SHA	6	0.441	0.983	0.156	96.519	0.26	0.42	-	0.33
SDA	6	0.664	0.979	0.173	78.179	0.33	-	0.38	0.30
HDA	6	0.543	0.979	0.175	76.649	-	0.42	0.32	0.26
SHDA	6	0.665	0.985	0.148	107.106	0.20	0.31	0.27	0.22

ONC: Optimal number of components.

SEE: Standard error of estimate.

F: F-test value.

PLS: partial least squares.

S: Steric.

H: Hydrophobic.

D: Hydrogen bond donor.

A: Hydrogen bond acceptor.

*: CoMISA model selected for 3D-QSAR simulation.

Table 3

Observed and predicted activities of EGFR ligands using the constructed CoMFA and CoMSIA models.

		CoMFA		CoMSIA
Comp.	Observed pIC50 *	PredictedpIC50	Residual	PredictedpIC50	Residual
1	6.620	6.571	0.049	6.600	0.020
2	7.081	7.192	−0.111	7.230	−0.149
3	7.260	7.234	0.026	7.147	0.113
4	6.638	6.394	0.244	6.522	0.116
5	8.102	8.337	−0.235	8.275	−0.173
6	8.721	8.508	0.213	8.493	0.228
7	6.060	5.940	0.120	6.012	0.048
8	6.180	6.237	−0.057	6.247	−0.067
9	7.000	6.893	0.107	6.952	0.048
10	6.721	6.828	−0.107	6.717	0.004
11	7.201	7.293	−0.092	7.322	−0.121
12	8.208	8.149	0.059	7.806	0.402
13	9.108	9.077	0.031	9.167	−0.059
14	9.018	9.059	−0.041	9.023	−0.005
15	8.638	8.563	0.075	8.566	0.072
16	7.252	6.377	0.875	6.012	1.240
17	7.244	7.210	0.034	7.159	0.085
18	7.796	7.790	0.006	7.710	0.086
19	7.620	7.744	−0.124	7.729	−0.109
20	8.194	8.089	0.105	8.216	−0.022

*: Experimental values of ligand bioactivity adapted from Ref [53].

Figure 4

Correlation of observed and predicted activity (pIC50) using 3D-QSAR models.

(A) CoMFA and (B) CoMSIA.

Ligand activities of Iressa® and the TCM candidates can be predicted based on structural conformation to the 3D-QSAR feature map, including features in steric field, hydrophobic field, and H-bond donor/acceptor characteristics. As illustrated in Figure 5, Iressa and the TCM candidates were able to match the generated 3D-QSAR model features. The benzene in Iressa® favored steric and hydrophobic fields, and H-bond was favored between its amine group and Asp855. In 2-O-Caffeoyl tartaric acid, the benzene structure favored steric and hydrophobic fields, and the carboxyl group favored H-bond formation with Lys716 and Lys728. The carbon chain structure in Emetine contoured to the steric and hydrophobic fields, and the amine group favored H-bond formation with Asp855. Rosmaricine had benzene and isopropyl structures that favored steric and hydrophobic fields, and an amine group that favored H-bond with Asp 855. The benzene structure in 2-O-feruloyl tartaric acid favored steric fields and the carboxyl group favored H-bond formations with Lys716 and Lys728. Iressa® and the TCM candidates have structural components that contour to the features of the 3D-QSAR model, thus were likely to be biologically active.

Figure 5

Structural contouring of different compounds to 3D-QSAR mapping.

(A) Iressa, (B) 2-O-Caffeoyl tartaric acid, (C) Emetine, (D) Rosmaricine, and (E) 2-O-Feruloyl tartaric acid. 3D-QSAR features are represented by the following: steric field favor/disfavor (green/yellow), hydrophobic field favor/disfavor (cyan/white), and hydrogen bond donor avor/disfavor (magenta/orange). Bond distances are shown in yellow.

Binding stability of the control and TCM candidates was validated using MD simulation. RMSDs of protein-ligand complex (Figure 6a) and individual ligand (Figure 6b) stabilized after 10 ns. The RMSDs of the protein-ligand complexes stabilized at approximately 1.6Å. With regard to individual ligands, the RMSDs of Iressa and 2-O-caffeoyl tartaric acid was 2.0 and 1.6Å, respectively. All other compounds registered RMSD values of approximately 1.0Å. The lower RMSD values of the TCM candidates suggest more stability within the receptor compared to Iressa. The energy trajectory of each compound is shown in Figure 6c. Complexes formed by Rosmaricine and 2-O-feruloyl tartaric acid had the lowest total energy (<−14,800 kcal/mol), followed by Iressa® and Emetine (approximately −14,700 kcal/mol), and 2-O-caffeoyl tartaric acid (−14,600 kcal/mol). Stabilization of total energy in ligand-protein complexes was achieved after 12 ns.

Figure 6

RMSD and total energy during MD simulation.

(A) Protein-ligand complex RMSD (Å), (B) ligand RMSD (Å) and (C) total energy of protein-ligand complex.

H-bond distance profiles in the EGFR receptor were summarized in Figure 7. A single H-bond between the amine group on Iressa® and the carboxyl group on Asp855 was formed after 9.74 ns and stabilized after 20 ns (Figure 7a). Two H-bonds were formed between the carboxyl group of 2-O-caffeoyl tartaric acid and Lys716 and Lys728 of the EGFR receptor (Figure 7b). The formation of two H-bonds contributed to a higher stability between 2-O-caffeoyl tartaric acid and the EGFR receptor. However, an increase in H-bond distance was observed towards the end of the 20 ns simulation period, suggesting a weakening of the H-bond at Lys728. Emetine formed a total of four H-bonds with the receptor, two with Asp722 and two with Ala855 (Figure 7c). Bond distances stabilized after 10 ns for Ala722 and 4 ns for Asp855. Rosmaricine formed three H-bonds each at Asp841 and Arg855 (Figure 7d). The multiple H-bonds enabled Rosmaricine to remain in a stable state within the protein. 2-O-Feruloyl tartaric acid also formed multiple H-bond at Lys716 and Lys728, enhancing its stability within the receptor site (Figure 7e). However, similar to 2-O-caffeoyl tartaric acid, an increase in H-bond distance was also observed at Lys728 for 2-O-feruloyl tartaric acid. These observations suggest that the bond at Lys728 weakens throughout the MD simulation process, and that the H-bond at Lys716 may be the primary bond for 2-O-caffeoyl tartaric acid and 2-O-feruloyl tartaric acid. In addition, periodic fluctuations in H-bond distances were observed in 2-O-caffeoyl tartaric acid, Rosmaricine, and 2-O-feruloyl tartaric acid. These phenomena can be attributed to the rotation of the amine group where the H-bond is formed. These MD results support our docking findings which identify Asp855, Lys716, and Lys 728 as key residues for docking.

Figure 7

Hydrogen bond distance profile during MD simulation.

(A) Iressa, (B) 2-O-Caffeoyl tartaric acid, (C) Emetine, (D) Rosmaricine, (E) 2-O-Feruloyl tartaric acid.

As determined in the CoMSIA model, hydrophobic interactions were key factors contributing to ligand bioactivity. Toward the final 20 ns of analysis, hydrophobic amino acids surrounding the docking region were Leu718, Val726, Ala743, Cys775, Phe795, Cys797, and Leu844. The hydrophobic subgroups of Iressa®, Emetine, and Rosmaricine were surrounded by Val726, Cys797, and Leu844 (Figure 8a). Hydrophobic groups of 2-O-caffeoyl tartaric acid were also surrounded Val726, Cys797, and Leu844 (Figure 8b). The hydrophobic region of 2-O-feruloyl tartaric acid was attracted to the Phe795 on EGFR (Figure 8b). The significance of matching the hydrophobic region of the ligand to that of the receptor may be to increase stability of the ligand-protein complex, and contribute to the bioactivity of the activated ligand. Our results indicate that Iressa® and the TCM candidates remained stable within the EFGR hydrophobic area following MD simulations.

Figure 8

Hydrophobic interactions of different compounds in EGFR.

(A) Iressa (orange), Emetine (green), and Rosmaricine (violet) (B) 2-O-Caffeoyl tartaric acid (blue), and 2-O-feruloyl tartaric acid (green). The hydrophobic regions in the binding site are depicted in red, and specific hydrophobic amino acids close to TCM compounds are indicated in red.

Conclusion

Structural and ligand based methods supported 2-O-caffeoyl tartaric acid, Emetine, Rosmaricine, and 2-O-feruloyl tartaric acid as potential EGFR inhibitors. Structurally, the TCM candidates were capable of forming H-bonds with key residues Asp855, Lys716, and Lys728 and matched hydrophobic regions of the receptor. Bioactivity of the candidates were evaluated using validated MLR, SVM, CoMFA, and CoMSIA models. All models indicated that the TCM candidates have good predicted bioactivity. Molecular simulation results further supported the high potential for the TCM candidates in drug development. Iressa®, the drug currently used clinically, bound to the ERGF receptor through a single H-bond at Asp855. In comparison, multiple H-bonds formed at Asp855 and additional H-bonds formed at Ala722 and Arg841 increase the stability of Emetine and Rosmaricine, respectively. The ability of carboxyl groups in 2-O-caffeoyl tartaric acid and 2-O-feruloyl tartaric acid to form multiple H-bond networks that directly blocked the ATP binding site was also a unique characteristic worthwhile of further investigation. Contour to hydrophobic regions of the TCM candidates within the receptor site provides additional support for the stability of the protein-ligand complex. In summary, using different simulation and validation methods, we have identified four TCM compounds that may have potential as novel EGFR inhibitors. As the four TCM compounds have two distinctive types of binding locations and bond formation within the EGFR binding site, we suggest exploring the possibility of connecting Emetine/Rosmaricine with 2-O-caffeoyl tartaric acid/2-O-feruloyl tartaric acid through a spacer. The connection could allow more of points of attachment, which in turn would contribute to more stable binding within the tyrosine kinase site.

Molecular structures and biological activities of ligands used for model training. Structural details of the 20 ligands adopted for ligand-based studies are listed within this table. Experimental bioactivity values for each ligand were adapted from Ref [53].

(DOC)

Click here for additional data file.