Exploring the interaction between epidermal growth factor receptor tyrosine kinase and some of the synthesized inhibitors using combination of in-silico and in-vitro cytotoxicity methods.

Quinazoline derivatives are potent inhibitors of human epidermal growth factor receptor (EGFR) as anticancer agents. In this study, the cytotoxic effects of a new series of synthesized quinazoline derivatives were evaluated using MTT assay against MCF-7 and HT-29 cell lines. Using molecular docking, the binding modes of all compounds were analyzed at the binding site of EGFR. Based on the results, the compounds L1, L2, L4, L5, L6, L7, L10, L15, and L18 may be promising EGFR inhibitors based on docking score and hydrogen bonds. Consistent with the experimental data, Met769 is recognized as a key residue in the binding of potential inhibitors. According to the MTT cytotoxicity assays, Lipinski's rule of five (RO5), absorption, distribution, metabolism, excretion, and toxicity (ADMET) parameters, and docking studies, three compounds L4, L15, and L10 with IC50 values of 80, 60, and 1 μM against the MCF-7 were selected for further comparative assessments. The dynamics of free EGFR, and selected ligand-EGFR complexes were investigated using molecular dynamics (MD) simulation studies. The results indicated that the three compounds bound to EGFR active site in a stable manner during the simulation through the formation of new hydrogen bonds with Phe699, Leu694, Gly700, Lys721, Met769, Arg817, and Asp831 with the superiority of compound L15. These features can promote future drug candidate designing to produce better derivatives in the search for the anticancer agents.

INTRODUCTION

Epidermal growth factor receptor (EGFR) is the cell-surface receptor of protein tyrosine kinase family which is over-expressed in numerous human tumors, containing ovarian, prostate, breast, bladder, lung, and colon(12). Several quinazoline derivatives such as erlotinib, gefitinib, and lapatinib have been synthesized as reversible inhibitors of tyrosine kinase(3). Due to developed resistance to first-generation drugs, irreversible inhibitors have been introduced with limited clinical efficacy(3). So far, numerous studies have been targeted at finding new structures based on the modification of quinazoline as potent inhibitors of EGFR in many anticancer studies(456). In addition, various compounds with substituted quinazolinone structural motif were reported as potent inhibitors of EGFR which were overexpressed in the number of cancer cell lines(78). Thus, quinazoline/quinazolinone-containing compounds represent an attractive scaffold for designing anticancer drugs(45678910).

In two previous studies, our research group have synthesized 4-anilinoquinazoline and 4(3H)-quinazolinones Schiff base derivatives reporting the antimicrobial activity(1112). As explained above, these derivatives may be regarded as new inhibitors of human EGFR(45789).

In the present study, cytotoxic effects of the 19 newly synthesized compounds were investigated against two cancer cell lines containing human breast adenocarcinoma (MCF-7) and human colon adenocarcinoma (HT-29) cell lines by MTT assay. The reason for choosing two types of cell lines was that variable levels of EGFR, human epidermal growth factor receptor 2 (HER-2), and/or HER-4 express by the growth of MCF-7 and HT-29 cell lines(19).

Molecular docking and molecular dynamics (MD) were used to recognize the binding sites on targets, conformational changes of biomolecules, stability, and protein folding(13). Therefore, in order to evaluate whether the synthesized compounds can potentially bind to EGFR an analysis of drug-receptor interaction was carried out by molecular modeling approaches.

To determine the suitability of the newly synthesized compounds, physicochemical properties are very important to proceed for further modifications until they achieve clinical trials(14).

Significant molecular properties for drug's pharmacokinetics in the human body were defined by the Lipinski's rule of five (RO5)(15).

Absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties which establish the pharmacokinetic outline of a drug molecule, are very crucial in evaluating its pharmacodynamics activities(16). As a result, the best compounds were introduced after evaluating its cytotoxic effects, molecular docking, in silico prediction of ADMET properties following Lipinski's RO5. Free EGFR, complexes of EGFR, and three compounds (L4, L15, and L10) were chosen for MD simulation to understand the ligand-receptor possible intermolecular interactions in details.

MATERIALS AND METHODS

Two cell lines were obtained from Pasteur Institute of Iran (Tehran, I.R. Iran). RPMI-1640 was prepared as culture media. 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT), and dimethyl sulfoxide (DMSO) were purchased from Merck, Germany. Absorbance was measured using an ELISA plate reader (ELX 808, USA).

Biological studies

The cytotoxic activities of synthesized compounds against MCF-7 and HT-29 cell lines were studied by MTT as previously described(2). The cell suspension was treated with 20 μL of various concentrations (0.1-100 μM) of synthesized samples in 1% DMSO diluted with culture medium. Erlotinib and DMSO were employed as the positive control and negative control, respectively. The experiments were repeated three times. The percentage of cell viability was calculated by reported equation(2). The half maximal inhibitory concentration (IC50) values in μM against both mentioned cell lines are reported in Table 1.

Table 1

Estimated free energy of binding (▵Gb) and inhibition constants (Ki) for 19 novel (A) 4-anilinoquinazoline, (B) 4(3H)-quinazolinones derivatives, and (C) erlotinib as potential epidermal growth factor receptor inhibitors.

Compliance of synthesized compounds to criteria of the prospective drugs

Lipinski's RO5 is a rule of scan to evaluate drug-likeness or determine whether a chemical compound with a given biological or pharmacological activity has properties that lead to a possible oral drug in humans(15). Besides, other significant conditions and essential elements of pharmacokinetics are ADMET properties(14). Today a lot of online tools and offline software programs are available, which can help in predicting these properties of a drug candidate. The properties displayed in Table 2 are calculated using the SwissADME website(17). The ADMET profiling indications were based on statistically derived from pkCSM(16)or Tox Prediction web server(18). Here, this program was applied to calculate the essential parameters, including Ames toxicity, hepatotoxicity, P-glycoprotein substrate, and median lethal dose (LD50) (mol/kg), etc.

Table 2

Drug-likeness prediction and absorption, distribution, metabolism, excretion, and toxicity prediction.

Docking studies

The theoretical binding mode of synthesized compounds to EGFR using the AutoDock 4.2.2 software(19). The target protein was derived from the RCSB protein data bank (PDB)(20). Crystal's structure of the EGFR tyrosine kinase domain with erlotinib (PDB code 1M17) at resolution 2.6 Å imported into AutoDock(21). The two-dimensional (2D) structures of the compounds were optimized using HyperChem 7.0 software as reported in previous studies(2223).

Docking was performed with the procedure similar to the previous works(2224). Also, to ensure the validity of docking, erlotinib as the well-known 4-anilinoquinazoline inhibitor was re-docked to the binding site with the root mean square deviations (RMSD) value relative to the crystal structure of 0.54 Å. The grid box was centered on Cα of Met769. Grid box dimensions were 50 × 50 × 50 (all in Å) with a 0.375 Å grid point spacing. A Lamarckian genetic algorithm (LGA) program was used to calculate 200 different conformers(1924).

Molecular dynamics simulations

MD simulation was used to explore the characterization of free EGFR and interaction between EGFR and three different inhibitors. Topology file for compounds L4, L15, and L10 were generated by the PRODRG server(25). MD simulations were performed using the GROMACS package 5.3.1(26). The “pdb2gmx” program was used to generate the topology file for protein using the G43a1 united-atom force field(27). The systems inserted into dodecahedron box and was solvated with the single point charge (SPC216) water model(28). The complex 1 (EGFR-L4), complex 2 (EGFR-L15), complex 3 (EGFR-L10), and free protein were neutralized by adding two, and one negatively charged Cl counter ion, respectively. The systems were introduced to energy minimization. The Berendsen algorithm was selected for thermostat and barostat in equilibration phase(29). The lengths of hydrogen-containing bonds were constrained using linear constraint solver (LINCS) algorithm(30).

MD simulations were executed in the NVT and NPT ensemble at a weak temperature coupling (τ = 0.1 ps) and pressure coupling (τ = 1 ps) in association with position restraint procedure. The protein, ligand, solvent, and the Cl ions were separately coupled to the thermostat with a reference temperature of 300 K and at the reference pressure of 1 bar. The thermostat and barostat were the Nosé-Hoover thermostat and the isotropic Parrinello-Rahman barostat for production step(29). The long-range electrostatic interactions were preserved by the particle-mesh ewald (PME) model(31). The MD simulations times were established to 40 ns with a time step of 2 fs in the periodic boundary condition. Temperature, potential and other analyses were deliberated using GROMACS package. The flexibility and stability of receptor and ligands during MD were inspected via the residue root mean square fluctuation (RMSF) and RMSD whose square fluctuation (RMSF) and RMSD whose plots are shown in Fig. 1. Computing the number of hydrogen bonds, solvent accessible surface (SASA), the radius of gyration (Rg) and other calculations were done as mentioned previously(32) (Table 3, Figs. 1 and 2). In order to display the binding mode obtained of ligand's trajectories after MD simulation, 3D and 2D graphical tools such as the visual MD (VMD)(33)the LigPlot software(34)were employed (Fig. 2).

Fig. 1

Analysis plots of the protein backbone and ligands structures for free epidermal growth factor receptor (EGFR), complex 1, 2, and 3 during 40 ns molecular dynamic simulation; A, the root mean square deviations (RMSD) plot; B, root mean square fluctuation (RMSF) plot; C, radius of gyration (Rg) plot; D, solvent accessible surface (SASA) plot.

Table 3

The obtained results from structural analysis of three complexes during molecular dynamic simulation.

Fig. 2

Analysis plots of the complexes 1, 2, and 3 during 40 ns molecular dynamic simulation; A, Minumum distance plot between hinge region residue Met769 and ligand; B, number of contacts plot of the hinge region residue Met769 and ligand; The H-bonding distribution plot vs time for C, L4-epidermal growth factor receptor (EGFR); D, L15-EGFR; and E, L10-EGFR.

RESULTS

MTT assay for cell viability/proliferation

All compounds were evaluated against MCF7 and HT-29 cell lines using MTT assay(2). Based on results presented in Table 1, the most of compounds were active against two cell lines. The cytotoxic activities of the compounds were comparable to erlotinib as a reference compound. Compounds L10 and L13 had potency similar to erlotinib against the MCF-7 cell line. In the other hand, replacement of H and Cl at L10 and L13 with 6,8-dichloro led to to compounds L18 and L16, respectively, which showed lower cytotoxicity than the first against one or two cell lines. Compared to erlotinib, compounds L11 and L12 were in the next priority against the MCF-7 cell line.

Drug-likeness prediction and absorption, distribution, metabolism, excretion, and toxicity prediction

All compounds have followed Lipinski's RO5 (Table 2). The ADMET predictions of some compounds have shown satisfactory results. Compound L15 violated no rule of the Lipinski's RO5 and met all ADMET parameters. In addition, compounds L1, L4, L10, L18, and L19, were not mutagenic based on the ADMET predictions and therefore may not be carcinogenic. Interestingly, all compounds were predicted to have absorption from the human intestine if administered orally based on SwissADME server. Toxic doses are often known as LD50 values in mg/kg body weight. The results of Tox prediction website showed that LD50 and predicted toxicity class of the three selected compounds L4, L15, and L10 may be 1500, 1230, and 1000 mg/kg of class IV with prediction accuracy 54.26%, 69.26%, and 67.38%, respectively. The LD50 was 125 mg/kg of class III for erlotinib. This implies that these compounds may have acceptable ADMET properties.

Molecular docking studies

The conformation with the lowest binding energy (~ -6 – -9 kcal mol-1) supports the idea that some of the compounds are well incorporated in the binding pocket. Dependent on the 4-anilinoquinazoline or 4(3H)-quinazolinone derivatives bearing Schiff base moiety, the interacting residues through hydrogen bonding with the studied compound differs. The binding patterns of different compounds are also slightly different, which may be responsible for the activity variations. It must be noted that inhibitors with 4-anilinoquinazoline scaffold have a common feature that in some cases they formed a hydrogen bond with the backbone NH of Met769 in the Hinge region(57). These compounds also were deeply embedded into EGFR via hydrophobic contacts that are conserved in the majority of the structures. These results were in consistent with the previously studied X-ray crystal structures, indicating the important roles of these(345789).

Confirmation of molecular docking by molecular dynamics simulation

The trajectory stability of the free EGFR and complexes 1-3 were confirmed by the analyses (Table 3, and Figs. 1,2). As shown in the RMSD plots (Fig. 1A), the trajectories were stable during the last 25 ns simulation. It is often considered that small RMSD values of a simulation indicate a stable state of the system (Table 3). The EGFR-complex 1 showed more deviation of RMSD average with regard to the free EGFR which relatively was in agreement with RMSF (Table 3). It may be due to the more interaction of EGFR atoms in the presence of L4 that causes some conformational and structural changes. This highlights the stable binding of the L4 with EGFR leading to instability of the protein. In addition, RMSD average shows primary fluctuations in the magnitude of RMSD of ligand's atoms during MD. Therefore, three compounds obtained an equilibrium state as described by the RMSD profile (Fig. 1A). To achieve a more detailed description of flexibility of the protein residues in the absence and presence of ligands, backbone RMSF values was computed (Table 3). Very few fluctuations were seen beyond 0.5 nm during the MD simulations (Fig. 1B). Only, residues 672 and 964 exposed relatively considerable fluctuations in the vicinity of N-and C-terminus. Higher local fluctuations followed at interface between domains IIA followed at the interface between domains IIA and IIB (the regions are specified with three bars) due to the presence of loops preceding domains (Fig. 1B). In addition, there are a few regions, including residues 850-863, and 888-907 that show higher flexibility than that of the free EGFR. It was found that the complexes 1 and 3 exhibited a little more fluctuation than the free EGFR (Table 3 and Fig. 1B). In other words, complex 2 has very low fluctuations identical to free EGFR, which again emphasis on stability of EGFR after L15 binding. It can be said that the binding of ligands did not affect the proteins' overall conformational diversity significantly, as there was no major change between the RMSD and RMSF values. The plots depicted from Rg of protein displayed that EGFR in four systems had a compact structure (Fig. 1C). The Rg values are compatible upon ligands complexation with respect to free EGFR (Table 3). It can be clearly seen that the Rg value of EGFR decreases slightly upon binding of ligands implying a more compact structure after the MD simulation. Rg value was more deviated in the complex 1 and had a more compact structure with more hydrogen bonds, which confirmed more changes of complex 1 in agreement with other results.

As shown in Table 3, complex 1 exhibited relatively higher deviations (142.63 ± 3.13) of SASA with time due to making more hydrogen bonds, while native EGFR structure showed higher values. However, curve for free EGFR did not differ seriously and preserved SASA, implying that the native conformation of EGFR was mainly conserved during the production time of simulation. These results were in line with Rg analysis. Together with SASA values reduction, Rg values decreased from that for starting structure to average, which along with the interaction mode showed the effective binding of L15 to EGFR.

To further confirm of stability and the effects of ligands binding, the DSSP algorithm was obtained. The main secondary structures of the EGFR maintain rather stable during 40 ns MD simulation. There were no significant changes in structural elements (Table 3). Although complex 1 showed the largest deviation of 179.75 ± 5.37 nm2 (Table 3). Very low levels of structural changes in the protein occurred due to ligands interactions especially for compound L15 that is in agreement with other results. This further confirmed the stability of complex 2. Hence L15 is partially affecting the structural conformation of EGFR.

In addition, the average of the minimum distances between “hinge region residue” Met769 in the active site and showed the superiority of L15 of complex 2, because this was very less than complex 1, and 3 (Table 3). It results in the formation of a stable hydrogen bond through Met769 as key residue with compound L15 (Fig. 2A and Fig. 3B). However, the results of this analysis could be compared with the results of the analysis of the number of contacts for more certainty (Table 3). The number of contacts for L4 in complex 1 is very less than that in complex 2, and 3 about key residue Met769 (Fig. 2B), but it is more for other residues in the active site (data not shown for brevity). More spatial prohibition of L4 is more notable than the interaction between L15 and EGFR due to other different interactions via six new hydrogen bonds (Fig. 3A). These residues are in accordance with other results(345789). Results of stability of hydrogen bonding in MD showed that in most conformers, one/two hydrogen bonds were found for complexes 2 and 3. This value reached to a maximum of 8 and an average of 4 for complex 1 Fig. 2C–2E).

Fig. 3

Schematic interaction resulting from LigPlot and VMD softwares for selected three compounds L4 (A and D), L15 (B and E), and L10 (C and F) after molecular dynamic simulations. In parts A, B, and C, the compounds exposure is blue-highlighted and hydrogen bonding is in green. Hydrophobic interactions are presented by red color. The hinge region residue Met769 is in green in D-F.

Temperature and pressure extended to plateau at 300 K and 1 bar (data not shown). The more negative total energy, potential energy, LJ (SR) energy of the complex 2 indicated that this complex was more stable with more hydrophobic interactions than complex 1, but the coulomb (SR) energy was more positive than complex 1 due to less hydrogen bonds (Table 3). Negative values of LJ (SR) and coulomb (SR) energies indicate the attraction between the receptor and the mentioned compounds. The number of hydrogen bonds between receptor and water was consistent with total energy in Table 3. The small changes of hydrogen bonds in the number of hydrogen bonds between protein and water, and intramolecular hydrogen bonds of protein indicated that there are the stable and soluble structures due to the high number of hydrogen bonds (Table 3).

Comparing the binding mode of three compounds before and after simulation showed that conformational changes of the main chain in Met769 led to the formation other hydrogen bonds (Fig. 3). Some new residues such as Phe699, Leu694, Gly700, Lys721, Met769, Cys773, Arg817, and Asp831 are positioned in proximity of ligands that could participate in formation of new six hydrogen bonds with L4 and one new hydrogen bond with L15 via Met769 at distance 2.68 Å in accordance with other results(345789). It must be noted that Lys 721 maintains a hydrogen bond with L10. Thus the end of MD simulation, new hydrogen bonds were established and previous hydrogen bonds, including Met769 and Thr830 were replaced (Fig. 3). Here the valuable application of MD simulations was indicated after docking of ligands in the binding site. Therefore, the used MD simulation was essential to specify geometry of EGFR in vicinity of water. As illustrated in (Fig. 3), compound L15 was located within this binding pocket. Thus, according to items listed above, L15 can be considered as a compound which maintains the EGFR structure stability throughout the simulation time after binding with the “hinge region residue” Met769. Overall, the analysis of MD simulation showed the superiority of complex 2 over complexes 1 and 3. The results obtained are consistent with the reports from Woods et al(35).

DISCUSSION

Cancer is one of the global health problems. Significant progress has been observed in cancer research in the last ten years(27). Considering the adverse complications and resistance to the current chemotherapeutic agents, development of more efficient anticancer agents with less harm has remained as an important subject in drug design(27). The findings of this study provide an integral part of the pre-clinical investigation of our compounds that can be further developed as anticancer agents against tumors, with excellent oral absorption. Based on The MTT results, showed that compounds L10-L13 substituted with OH, Br, Cl, and NO2 groups may play an important role in growth inhibition of the MCF-7 cell line with IC50 values of 1-10 μM compared to erlotinib. This could be attributed to the withdrawing effects on the heterocyclic ring.

The docking results showed that the residues Met769 and Thr830 play a pivotal role in making of the hydrogen bond with EGFR(57). The same conclusion to our docking results was reported for several other compounds(45789). Compared to erlotinib, some of these compounds properly fit into the ATP binding pocket in EGFR crystal structure, suggesting that they may be potential EGFR inhibitors.

Molecular properties, Lipinski's RO5, and ADMET parameters predicted that some compounds have acceptable properties and lack the sign of mutagenicity effect. Three compounds L4, L10, and L15 violated no rule of the Lipinski's RO5 and also were not mutagenic based on the ADMET predictions and therefore may not be carcinogenic. The activity of both compounds L4 and L15 may not be correlated to multi-drug resistant tumors due to not being P-glycoprotein substrates based on Swiss ADME. Therefore, the MTT cytotoxicity assays, Lipinski's RO5, ADMET parameters, and docking studies demonstrated that the three compounds L4, L15, and L10 were selected as the best compounds for further assessments.

Clearing up of ligand binding mechanisms is the essential step to achieve more selective compounds for a given target. Complex delivered in its natural environment is required to achieve more precise ligand-receptor models(2432). Therefore, MD simulations of complexes of the EGFR and three compounds with the weak, moderate, and strong cytotoxic effect were performed for further assessment of the effects of different ligand binding on the conformation of EGFR The three complexes immersed in water molecules stay in equilibrium during MD. We also found that the MD simulation created an improved and more relaxed structure of enzyme and ligands.

However, compound L15 exhibited nearly appropriate properties comparable to those of the standard drugs, which include good potential interaction with EGFR confirmed by MD simulation above. Compound L15, has strong interaction with ”hing region key residue“ Met769 of EGFR target, which is involved in the anticancer treatment strategies. Moreover, the overall analysis of MD simulation showed the superiority of complex 2 over complexes 1, and 3. As a result, in accordance with other results, compound L15 with quinazolinone structural motif could be a potential lead compound for identification of EGFR inhibitors(78).

CONCLUSION

This paper provided an approach for studying the interactions of protein with new three ligands using molecular modeling techniques. This study demonstrated the utility of other quinazoline/quinazolinone derivatives as an attractive scaffold for the development of potential EGFR inhibitors. More improvements are needed in progress to optimize new quinazoline derivatives as anticancer drug candidates. The prediction methods employed in the current study help to accelerate the analysis of the designed compounds, which were developed in our laboratory before pre-clinical, time-consuming, and high risking experimental methods. It can also be used to select the best compound among some novel proposed derivatives. Hence, it is desired that biological evaluation of the designed compounds for their inhibitory activities on EGFR is conducted using enzymatic specific tests. This research highlights the potential inhibitory activity of selected derivatives for further development.