Discovery of Anilino-1,4-naphthoquinones as Potent EGFR Tyrosine Kinase Inhibitors: Synthesis, Biological Evaluation, and Comprehensive Molecular Modeling

Epidermal growth factor receptor (EGFR) has been recognized as one of the attractive targets for anticancer drug development. Herein, a set of anilino-1,4-naphthoquinone derivatives (3–18) was synthesized and investigated for their anticancer and EGFR inhibitory potentials. Among all tested compounds, three derivatives (3, 8, and 10) were selected for studying EGFR inhibitory activity (in vitro and in silico) due to their most potent cytotoxic activities against six tested cancer cell lines (i.e., HuCCA-1, HepG2, A549, MOLT-3, MDA-MB-231, and T47D; IC50 values = 1.75–27.91 μM), high selectivity index (>20), and good predicted drug-like properties. The experimental results showed that these three promising compounds are potent EGFR inhibitors with nanomolar IC50 values (3.96–18.64 nM). Interestingly, the most potent compound 3 bearing 4-methyl substituent on the phenyl ring displayed 4-fold higher potency than the known EGFR inhibitor, erlotinib. Molecular docking, molecular dynamics simulation, and MM/GBSA-based free energy calculation revealed that van der Waals force played a major role in the accommodations of compound 3 within the ATP-binding pocket of EGFR. Additionally, the 4-CH3 moiety of the compound was noted to be a key chemical feature contributing to the highly potent EGFR inhibitory activity via its formations of alkyl interactions with A743, K745, M766, and L788 residues as well as additional interactions with M766 and T790.


■ INTRODUCTION
Cancer is the second most common serious public health problem ranked behind cardiovascular diseases. 1 The burdens of cancer incidence and mortality are rapidly increasing worldwide, 2 which have rendered the discovery of novel effective anticancer drugs a research hotspot. Receptor tyrosine kinases are well known as key protein signaling regulators for many cellular events including cancer-related processes. 3 Currently, many tyrosine kinase inhibitors are approved as anticancer drugs and have taken a promising place in clinical cancer management. 3 Epidermal growth factor receptor (EGFR) is a member of the ErbB family of the receptor tyrosine kinases, 4 in which its upregulations are found in many types of cancers such as breast cancer, non-small cell lung cancer, cholangiocarcinoma, ovarian cancer, prostate cancer, leukemia, and colon cancer. 5−9 An overexpression of EGFR consequently activates prooncogenic downstream signaling pathways including PI3K/AKT/mTOR and RAS/RAF/MEK/ ERK. 10 These pathways play key roles in oncogenesis-related events, including cell survival, cell proliferation, cell differentiation, cellular apoptosis, and metastasis. 11,12 Accordingly, the EGFR is considered one of the attractive targets for cancer treatment. The first-generation tyrosine kinase inhibitors (i.e., erlotinib and gefitinib) were developed as reversible EGFR inhibitors. 13 However, these drugs are concerning for their acquired drug resistance 14 and considerable side effects. 15,16 Thus, the discovery of a novel class of EGFR inhibitors with improved efficacy but minimized toxicity is an area of urgent need. 17,18 Naphthoquinone, a derivative of naphthalene bearing two carbonyl groups, is a privileged scaffold found in many classes of naturally occurring bioactive compounds. Some of naphthoquinone-based compounds have been used as anticancer drugs, including doxorubicin and daunorubicin. 19,20 Various molecular mechanisms underlying the anticancer properties of the naphthoquinone-based compounds have been reported. 21−23 Among these, inhibition of the EGFR signaling pathway was noted for anticancer action against several types of cancer cells. Examples of naphthoquinones acting as EGFR inhibitors are shown in Figure 1. Furano-1,2-naphthoquinone I 24 and synthetic shikonin derivatives containing benzoylacrylic acid II 25 ( Figure 1) were reported to inhibit EGFR and induce cell cycle arrest in the G2/M phase, leading to cellular apoptosis. Some plant-derived 1,4-naphthoquinone derivatives with EGFR inhibitory effect were also reported. For instance, shikonin, a compound found in Lithospermum erythrorhizon, was reported to inhibit EGFR in human epidermoid carcinoma, 26 glioblastoma, 27 and non-small cell lung cancer cells. 28 Likewise, plumbagin, isolated from the roots of Plumbago zeylanica L., was reported to exert the same effect against in vitro pancreatic cancer cell line as well as in vivo study. 29 Notably, 1,4-naphthoquinones embedded with an amino group or substituted amino groups on the quinone ring have been extensively reported to possess anticancer properties. 30−37 However, their underlying mechanisms and molecular targets for anticancer effects are still rarely understood.
Understanding of target−ligand binding modes and interactions as well as drug-likeness (or pharmacokinetic properties) is essential for successful drug design and development. Computational approaches are well known as facilitating tools in drug discovery (i.e., increasing success rate, saving time and cost, and gaining beneficial key findings toward effective structural design and optimization). 38 Molecular docking and molecular dynamics (MD) simulation have been widely utilized for elucidating possible molecular targets and binding interactions. 39 Poor drug-likeness and severe toxicities are highly concerned as causative factors of late-stage failure in the development pipeline. Currently, in silico drug-likeness predictions have been used for prioritizing compounds to facilitate successful clinical development. 40 Examples of computer-aided development of novel anticancer agents were reported by our group 41,42 as well as others. 43 Based on the structural features of available EGFR inhibitors (i.e., naphthoquinones and erlotinib in Figure 1), it is anticipated that the naphthoquinone skeleton could mimic the quinaxoline ring of erlotinib. Modification of the core skeleton with various substituents possessing different hydrophobic, electronic, and H-bonding properties would contribute to the formations of additional interactions and enhanced binding with the target protein. Thus, a set of new naphthoquinone-based EFGR inhibitors was rationally designed by introducing another suitable functionality to the 1,4naphthoquinone system as follows: (i) using anilino-1,4naphthoquinone as a core structure, (ii) varying the substituents (R = electron donating and electron withdrawing groups) on the aniline ring, and (iii) varying the substituents (X = H, Cl, and Br) on the quinone ring. A set of amino-1,4naphthoquinone derivatives 3−18 (Scheme 1) was synthesized, and their anticancer activities were evaluated against six cancer cell lines. Cytotoxicity of the compounds against the normal cell line (MRC-5) was also performed, and the selectivity index (SI) was calculated. Anticancer activities and SI values of the compounds were primarily considered to prioritize a set of compounds for further in vitro investigation on EGFR inhibitory activity. Drug-likeness of the selected compounds was also predicted using in silico tool to ensure their potential of successful development. Subsequently, molecular docking and MD simulation were carried out to elucidate possible binding modalities and interactions. Taken together, the findings from this study would be beneficial for further design and development of anilino-1,4-naphthoquinone-based compounds as novel anticancer drugs targeting the EGFR tyrosine kinase.
In overview, six compounds (4, 8−9, and 12−14) exhibited anticancer effects on all types of tested cancer cell lines (with varied potency), whereas two compounds (5 and 18) were inactive. However, all active compounds were less potent than the reference drug, doxorubicin.
For the A549 cell line, most of the compounds displayed weak activity (IC 50 values greater than 40 μM), except for  Table S1) for further investigations. The most potent compounds are shown in bold. NC, noncytotoxic. IC 50 > 50 μg/mL denoted as noncytotoxic; ND, not determined. For the MOLT-3 cell line, most compounds in series I (3− 4, 7−10, and 12) and series II (13−15) were active anticancer agents with IC 50 values of <10 μM, whereas all compounds in series III (16−18) were inactive. Enhanced cytotoxic effects were observed for the compounds in series I with electronwithdrawing groups (i.e., R = nitro and cyano moieties) as shown by the lower IC 50 values (<3 μM). Promisingly, cyano derivative 9 was the most potent compound among all (IC 50 = 1.75 μM). In contrast, the replacement of the methyl (R) group in compound 3 with halogen and OMe groups leads to analogues with lesser potencies (bromo derivative 5 was inactive, iodo derivative 6: IC 50 = 42.58 μM, and OMe derivative 11 IC 50 = 17.59 μM). Similar potency impairments were observed for fluoro (compound 4), trifluoromethyl (compound 7), and hydroxyl (compound 12) derivatives but with lesser degree. Notably, the anticancer effects were enhanced when the chloro (X) group on the naphthoquinone core in series I (compounds 4 and 5) was replaced with X = bromo to provide derivatives 14 and 15 (series II). This suggested that the substituents (X = Cl, Br) on the naphthoquinone core are essential for anticancer activity since the cytotoxic effect was lost for the compounds without these substituents as observed for inactive compounds of series III (16−18).
For the HepG2 cell line, the most potent compound was a bromo derivative (X) of series II compound 15 (IC 50 = 4.94 μM). It should be noted that the presence of a bromo (X) group of 15 on the naphthoquinone core may be essential for potent anticancer activity against HepG2. This effect was seen when the absence of activity was observed for the inactive compound 5, which contains chloro group (X) with the same Br substituent (R). A similar effect was noted when comparing two fluoro (R) derivatives of series I and II (compound 4: IC 50 = 131.89 μM and compound 14: IC 50 = 19.07 μM). Additionally, other active compounds in series II (13−14) and III (16−17) showed comparable activities (19.07 μM < IC 50 < 25.07 μM), whereas those of series I, except for compound 12, showed lower activity (IC 50 > 73.37 μM). Unlike MOTL-3, it was suggested that the substitution of the chloro group on the naphthoquinone core gave the  compounds weak/inactive anticancer activity against the HepG2 cell line.
Considering the cytotoxic effects on breast cancer cell lines, the most promising compound against MDA-MB231 was compound 12 (R = OH of series I, IC 50 = 4.77 μM) followed by two compounds from series II (i.e., compounds 14 and 13 with IC 50 = 7.54 and 9.48 μM, respectively). In the case of T47D cells, the most potent agent was compound 13 (IC 50 = 10.12 μM), followed by two compounds of series I (compounds 9 and 8 with IC 50 = 11.27 and 11.68 μM, respectively).
Cytotoxicity of all synthesized compounds (3−18) was investigated against the normal embryonic lung (MRC-5) cell line (Table 1), and their SIs were calculated (Table S1). Most of the synthesized compounds (5−7, 11, 15−16, and 18) were noncytotoxic against the normal cell line MRC-5; however, these compounds showed relatively weak or inactive activities against most of the tested cancer cell lines. Potent anticancer activities along with the SI values were simultaneously used for selecting a set of compounds for further studies. As a result, three compounds from series I (i.e., compounds 3, 8, and 10) with promising anticancer activities and high SI values were selected for further studies.
Drug-likeness Prediction. To ensure the potential of successful development, the drug-likeness of three selected compounds (3, 8, and 10) along with the known EGFR inhibitor, erlotinib, was predicted using in silico SwissADME web tool. 44 The investigated properties include molecular weight (MW), number of hydrogen bond donors and acceptors (HBD and HBA), number of rotatable bonds (RB), topological polar surface area (TPSA), and lipophilicity (log P). It was found that the predicted values of all investigated compounds (3, 8, and 10) were within the range of the Lipinski's rule of five criteria: (i) MW ≤ 500 Da, (ii) HBD ≤ 5 and HBA ≤ 10, (iii) RB ≤ 10, (iv) TPSA ≤ 140 Å 2 , and (v) log P ≤ 5, as summarized in Table 2. 45 The findings suggested that these compounds are drug-like compounds, which could likely be further developed as anticancer drugs.
It should also be noted that although many naphthoquinone-based compounds were proposed to inhibit EGFR activity via inhibition of EGFR-mediated signaling pathways, 24,26−29 the mechanism underlying the inhibitory effect on the EGFR protein has still not been fully explored in vitro. Thus, this is the first study to shed light on the potentiality of anilino-1,4-naphthoquinone derivatives as promising EGFR tyrosine kinase inhibitors.
Computational Studies. Molecular Docking. Molecular docking was performed to reveal the possible binding modes of compounds 3, 8, and 10 against the EGFR tyrosine kinase (PDB ID: 1M17). 56 Results showed that all naphthoquinones could accommodate within the adenosine triphosphate (ATP)binding pocket of EGFR in the same manner as erlotinib, a cocrystallized ligand (Figure 3). The naphthoquinone core of all compounds (3, 8, and 10) was positioned near the hinge region domain (red) 56 to mimic the quinaxoline ring of erlotinib, whereas their phenyl ring of all naphthoquinones pointed toward the Cα-helix domain (green) of the EGFR (Figure 3).
The 2D protein−ligand interaction diagrams were generated for an in-depth understanding of ligand−protein interactions  ( Figure 4). It was revealed that compound 3 could form similar interactions as erlotinib. The key interactions included (i) the formation of H-bonding between its carbonyl group and M793 residue at the hinge region 57 and (ii) the formation of alkyl interactions between its 4-CH 3 moiety and K745/M766/L788 residues ( Figure 4A, pink) to mimic the terminal alkyne group of erlotinib ( Figure 4D). It was also observed that the 4-CH 3 substituted aromatic ring of compound 3 played essential roles in the formations of additional interactions with the target protein, which could be a reason of its higher potency than the reference inhibitor, erlotinib (Table 1). These additional interactions were noted as (i) π−sulfur interaction with M766, (ii) π−lone pair interaction with T790, and (iii) alkyl interaction with A743. For both nitro compounds 8 and 10, the nitro phenyl groups formed H-bond and electrostatic interactions with the polar K745, E762, and T790 residues.
In overview, the results showed that all focused compounds (3, 8, and 10) could occupy the ATP-binding pocket of EGFR via the formations of van der Waals (vdW) force, H-bonding, π−π interaction, alkyl interaction, and attractive charge. Among these interactions, the vdW force seems to be the major force for protein−ligand complexation. It was shown that all compounds (3, 8, and 10) could form vdW interactions with residues similar to those formed with erlotinib. These key binding residues were also reported as key binding residues for many classes of small-molecule EGFR inhibitors such as vinyl sulfones, 49 thiazoles, 47 thiophenes, 57 chalcones, 58 and 6arylureido-4-anilinoquinazolines. 7 Based on the EGFR inhibitory activity (Table 3), only the most potent compound 3 was selected for further investigation on its atomistic binding mechanism against the EGFR tyrosine kinase in an aqueous environment using MD simulation and free energy calculation.
MD Simulation and Free Energy Calculation. System Stability. The stability and dynamic behavior along the simulation times of the most potent compound 3 in complex with the EGFR tyrosine kinase domain were determined using root-mean-square displacement (RMSD), number of atomic contacts (#contacts), and radius of gyration (R g ) calculations.
Results showed that the RMSD values of the EGFR/3 complex were slightly fluctuated at the first 20 ns, which was supported by the adaptation of atomic contacts within the 5 Å sphere of compound 3 during this period, and then reached the equilibrium state after 20 ns with an average RMSD value of ∼3 Å ( Figure 5). In the case of ligand 3, the RMSD values were increased at the first 5 ns and then were persistent with an average RMSD value of ∼1.5 Å until the last stage of MD simulation. According to the R g calculation, the compactness of the EGFR tyrosine kinase domain was found to relatively remain constant at ∼20−21 Å from the beginning till the end of the simulation, reflecting the stable closed conformation of the EGFR/3 complex ( Figure 6). Altogether, the results indicated that the molecular complexation between compound 3 and EGFR was highly stable in an aqueous environment.
Binding Hotspot. To investigate the key interacting amino acid residues upon the binding of compound 3 to the target EGFR, the per-residue decomposition free energy (ΔG bind,residue ) based on the molecular mechanics/generalized Born surface area (MM/GBSA) method was calculated using 100 snapshots taken from the last 20 ns of the MD simulation. Among residues 696−987, only residues 700−870 are shown in Figure 7, where the binding orientation of compound 3 inside the ATP-binding pocket of EGFR tyrosine kinase is illustrated in the right panel. The contributing amino acids and vdW contributions are colored according to their ΔG bind,residue and energy contribution values, respectively. Notably, only the residues exhibiting an energy stabilization value of ≤−0.5 kcal/ mol were considered.
Results revealed that six residues (i.e., L718, V726, A743, K745, G796, and L844) were associated with the binding interactions of compound 3 and the EGFR ( Figure 7A). It should be noted that these residues, except G796, along with M793 were previously characterized as the top five contact frequencies forming a core hydrophobic binding pocket, which was noted for the accommodation of the adenine base of ATP or the ATP-competitive inhibitors. 59 Among all, L718 exhibited the highest vdW contribution to the binding of compound 3 ( Figure 7B). The L718 residue, located within the p-loop of the EGFR, was reported to play key roles in ligand−protein binding of many small-molecule EGFR inhibitors. Its mutations (e.g., L718A and L718Q) were also noted to contribute to drug resistance. 60−62 Overall, the main energy contribution for protein−ligand complexation was derived from vdW forces rather than electrostatic attraction, which is similar to that of other reported EGFR tyrosine kinase inhibitors, including erlotinib, gefitinib, lapatinib, 63 vinyl sulfones, 49 thiazoles, 47 thiophenes, 57  . Notably, compound 3 was the most promising compound with 4-fold more potent inhibitory effect than the known EGFR inhibitor, erlotinib. Molecular docking revealed that the most promising compound 3 could occupy the ATP-binding pocket of the target protein EGFR in a similar manner as erlotinib. MD simulation and MM/GBSAbased free energy calculation showed that the molecular protein−ligand complexation was mainly driven by vdW forces and was highly stable in an aqueous environment. Additionally, the key amino acid residues associated with compound−target binding were revealed, including L718, V726, A743, K745, G796, and L844 residues. In summary, this study demonstrated the comprehensive utilization of chemical synthesis, in vitro, and in silico approaches for the discovery of new 1,4naphthoquinone-based EGFR inhibitors for potential development as anticancer drugs. The key findings regarding key chemical structural features and key protein residues participating in the ligand−target bindings would also be beneficial for the future discovery and design of the compounds in related classes as small-molecule EGFR inhibitors.

■ EXPERIMENTAL SECTION
Chemistry. Column chromatography was carried out using silica gel 60 (70−230 mesh ASTM). Analytical thin-layer chromatography (TLC) was performed on silica gel 60 F 254 aluminum sheets. 1 H-and 13 C NMR spectra were recorded on a Bruker AVANCE 300 NMR spectrometer or a Bruker AVANCE NEO 500 NMR spectrometer. Fourier-transform infrared spectra were obtained using universal attenuated total reflectance attached on a PerkinElmer Spectrum One spectrometer. HRMS data were recorded on a Bruker Daltonics (microTOF) instrument. Melting points were determined using a Griffin melting point apparatus and were uncorrected.
General Procedure for the Synthesis of Naphthoquinones (3−21). A mixture of 1,4-naphthoquinone 1 (2.4 mmol) and the appropriate aniline derivative 2 (2.0 mmol) in absolute  ethanol (20 mL) was stirred under reflux until completion of the reaction, as monitored by TLC, and then concentrated under reduced pressure. The crude product was purified using silica gel column chromatography to afford the pure product.   Cytotoxic Assay. Cells suspended in the corresponding culture medium were inoculated in 96-well microtiter plates (Corning Inc., NY, USA) at a density of 10,000−20,000 cells per well and incubated for 24 h at 37°C in a humidified atmosphere with 95% air and 5% CO 2 . An equal volume of additional medium containing either serial dilutions of the test compounds, positive control (etoposide and/or doxorubicin), or negative control dimethyl sulfoxide (DMSO) was added to the desired final concentrations, and the microtiter plates were incubated for an additional 48 h. The number of surviving cells in each well was determined using the MTT assay 73,74 (for  HuCCA-1, HepG2, A549, MDA-MB-231, T47D, and MRC-5 cells) and the XTT assay (for MOLT-3 cells). The IC 50 value is defined as the drug (or compound) concentration that inhibits cell growth by 50% (relative to negative control). The compounds exhibiting IC 50 >50 μg/mL were considered as noncytotoxic.
Drug-likeness Prediction. Drug-likeness of potent naphthoquinones as well as erlotinib, the known EGFR inhibitor, was evaluated using the SwissADME web tool. 44 ADP-Glo Kinase Assay. The EGFR kinase inhibition activity of the target compounds was evaluated using the ADP-Glo kinase assay (Promega, Wisconsin, USA). First, 8 μL of buffer (40 mM Tris-HCl pH 7.5, 20 mM MgCl 2 , and 0.1 mg/ mL bovine serum albumin) was added to a 384-well plate. Second, 5 μL of either wild-type EGFR or L858R/T790M mutant EGFR (1.25 ng/μL, Sigma-Aldrich) and 2 μL of inhibitors were added, followed by 10 μL of a mixture of 25 μM ATP and 12.5 μg/mL poly(glu-tyr), and incubated for 1 h at room temperature. Third, 5 μL of the ADP-Glo reagent was added and incubated for 40 min to terminate the kinase reaction. Finally, 10 μL of kinase detection reagent was added and incubated at room temperature for 30 min. The luminescence was detected using a microplate reader (Infinite M200 microplate reader, Tecan, Mannedorf, Switzerland). All assays were performed in triplicate. The relative inhibition (%) of inhibitors was calculated in comparison to the control with no inhibitor as shown in eq 1. Computational Studies. System Preparation. The crystal structure of EGFR complexed with erlotinib (PDB ID: 1M17) 56 was downloaded from Protein Data Bank (PDB). The protonation state of all ionizable amino acids was determined at pH 7.0 using PROPKA 3.0. 75 The 3D structures of 3, 8, and 10 were created and fully optimized using the Gaussian 09 program 76 with the HF/6−31d basis set. 47,77,78 The protonation state of all studied ligands was determined at pH 7.0 using MarvinSketch program. 79 The site of erlotinib in the EGFR crystal structure was defined as the docking sphere (15 Å) for molecular docking analysis using the CDOCKER module 80 implemented in Discovery Studio 2.5 program. The binding mode between EGFR and inhibitor(s) was visualized using the Discovery Studio Visualizer. The docked complexes with the lowest CDOCKER interaction energy of each system were chosen as the initial structure for MD simulation. The electrostatic potential (ESP) charges were calculated with the HF/6-31d basis set and then converted to restrained ESP charges using the antechamber module of AMBER20. The AMBER ff14SB 81 and the general AMBER force field version 2 (GAFF2) 82,83 force fields were applied for protein and ligand(s), respectively. All missing hydrogen atoms of the protein were added using the LEaP module. Then, each system was neutralized by the chloride ions and solvated using the TIP3P water model. 84 The added hydrogen atoms and water molecules were subsequently minimized using 1000 steps of the steepest descent followed by 2500 steps of conjugated gradient methods. Finally, the whole system was minimized using the same minimization process.
MD Simulation and Free-Energy Calculation. The MD simulations were carried out by a time step of 2 fs using AMBER20 program. The short-range cutoff for nonbonded interactions was set as 10 Å, while the particle mesh Ewald summation approach was applied to treat long-range electrostatic interactions. 85,86 Temperature and pressure were controlled by a Langevin thermostat and the Berendsen weak coupling algorithm, respectively. The SHAKE algorithm was used to constrain all covalent bonds involving hydrogen atoms. 87 In the relaxation phase, all of the models were gradually heated up from 10 to 310 K for 100 ps with the application of a harmonic restraint of 30.0 kcal/mol Å 2 to the EGFR−ligand complex. In the next equilibrium phase, each complex was subjected to restrained MD simulations at 310 K with the harmonic restraint of 30, 20, 10, 5, and 2.5 kcal/mol Å 2 for 500 ps in total, followed by unrestrained MD at 310 K for 500 ps. Subsequently, MD simulations were performed under the NPT ensemble (310 K and 1 atm) until reaching 100 ns. The CPPTRAJ module 88 of AMBER20 was used to compute the structural information, including RMSD, #contacts, and R g . The ΔG bind,residue was calculated using the MM/GBSA method on 100 MD snapshots extracted from the last 20 ns of the MD production phase.
■ ASSOCIATED CONTENT * sı Supporting Information