Designing, molecular docking, and dynamics simulations studies of 1,2,3-triazole clamped Uracil–Coumarin hybrids against EGFR tyrosine kinase

Since the last decade, hybrid drug strategies have attracted many researchers for their improved anti-cancer potential in comparison to single drug components. Complying to this approach, 28 novel Uracil–Coumarin hybrids with different sized linkers (2–5 carbon atoms) and substituents were designed to occupy the active site of protein epidermal growth factor receptor (EGFR) tyrosine kinase (Protein Data Bank ID: 1M17). Molecular docking studies were performed for all ligands (A1-D7) to identify the potential candidate using Schrödinger software. The relative binding affinity of hybrids toward EGFR was compared with standard Erlotinib on the basis of gScore and Emodel score. Positively, all the hybrids docked inside the cavity and showed significant interactions, compounds A6, A2, and A7 with short-chain linker (two carbon atoms) and halogen substituents were found to have more interactions and better docking score than standard Erlotinib. The visualization results depicted that compound A6 showed the highest affinity and formed the best binding pose to the target EGFR with gScore = −8.891 kcal/mol and Emodel score = −100.744 in comparison to standard Erlotinib (gScore of −8.538 kcal/mol and Emodel score = −80.588). Moreover, a molecular dynamics study also reveals that ligand A6 forms a stable complex with root mean square deviation (RMSD) of 0.3 nm and the plateau phase started just after 10 ns (time). Hence, the present research provides computational insights of Uracil–Coumarin hybrids as potential ligands against EGFR tyrosine kinase and in future in vitro investigations of these hybrids may prove their therapeutic potential against cancer.


INTRODUCTION
Cancer is a diversified disease, mainly recognized by unrestrained cell division that can lead toward mortality by acquiring metastatic form. Cancer has been declared the second most deadly disease that makes one out of eight deaths worldwide. According to the latest global cancer data released by the International Agency for Research on Cancer, nearly18.1 million new cases and 9.6 million deaths were reported in the year 2018 (WHO, 2018). Approximately 13.2 million cancer-related deaths and 21.4 million new cancer cases are expected by 2030 (Siegel et al., 2019). Unfortunately, all currently available chemotherapeutic agents suffer from undesirable off-target effects and frequently develop resistance. Thus, the development of potent and effective novel anti-neoplastic drugs is one of the most intensely persuaded goals of existing medicinal chemistry.
Epidermal growth factor receptor (EGFR) is a transmembrane tyrosine kinase receptor that belongs to the human epidermal growth factor receptor family. EGFR mediates the pathways of cell proliferation, growth, survival, and also involved in drug resistance. Pertaining to such involvement in cancer progression, EGFR has become an imperative target to discover novel chemotherapeutic agents (Sigismund et al., 2018). EGFR is highly expressed in many types of cancer cells; a study reveals that EGFR is over-expressed in 14%-91% of breast cancer cases (Fox et al., 2007;Klijn et al., 1992). Advantageously, EGFR provides a larger hydrophobic pocket in which bigger molecules like Erlotinib (Yun et al., 2007) can be docked easily.
Coumarin (2H-1-benzopyran-2-one) is found in the entire plant kingdom, including essential oils such as cinnamon bark oil and cassia leaf oil, and fruits, green tea, and foods such as chicory (Skehan et al., 1990). Coumarin derivative (Fig. 1), i.e., 7,8 dihydroxycoumarin (daphnetin) has been reported to potentially inhibit EGFR tyrosine kinase (Yang et al., 1999) and effectively bind in a hydrophobic pocket of EGFR with good docking score (Kumar et al., 2018). Moreover, Coumarin hybrids ( Fig. 1) containing 1,2,3-triazole ring have been reported to inhibit cancer cells with good IC 50 values (Pingaew et al., 2014;Singh et al., 2016;. Uracil (pyrimidine 2,4 dione) is a common nucleobase present in the nucleic acid of RNA. In RNA, Uracil makes two hydrogen bonds with adenine and in DNA sequencing thymine replaces Uracil nucleobase. Owing to have DNA intercalation properties, Uracils impart anti-cancer activity (Fadda and Pome 2011). Analogous to Coumarin, Uracil hybrids ( Fig. 2) containing 1,2,3-triazole moiety in their structure have also been proved to actively inhibit various cancer cells with lesser IC 50 values (Kumar et al., 2012;Ruddarraju et al., 2016). Concerning to EGFR inhibition, 5-Fluorouracil (5-FU) is reported to form a complex with epidermal growth factor (EGF) protein. This 5FU-EGF complex when docked with EGFR protein results in more affirmative binding (Kanagathara and Kavitha, 2015). Furthermore, 5-FU in combination with oxaliplatin also found to increase EGFR phosphorylation in human colorectal cancer cell lines (Schaeybroeck et al., 2005).
Furthermore, incorporation of 1,2,3-triazole ring between two pharmacophores is evidently a useful tool to produces bifunctional hybrid molecules with improved efficacy toward cancer target (Huber et al., 2009;Jagasia et al., 2009;Zhang et al., 2009). EGFR docking studies of some hybrids confirmed that the role of the triazole ring is in making water-mediated H-bond with Thr766, Lys721, Arg817, and Met769 of target's amino acid residues (Banerji et al., 2018;Nguyena et al., 2018).
The foremost cytotoxic potential of the moieties Uracil, Coumarin, and 1,2,3-triazole prompted us to design a scaffold of 28 novel bi-functional N 1 , N 3 Uracil hybrids linked with Coumarin via 1,2,3-triazole using different carbon chain linker (n = 1 − 4) to get the synergistic inhibition on EGFR. Docking study was performed to find out the potential ligand amid all designed hybrids, which have the best affinity toward the EGFR tyrosine kinase domain. Virtual screening (2D and 3D) of all potential hybrids was compared with standard Erlotinib, a well-known EGFR inhibitor. The binding efficacy of the ligands has been estimated on the basis of Glide score (gScore) and Emodel score because Glide uses Emodel to pick the "best" pose of a ligand (pose selection). Various interactions between ligand and protein were also studied and compared with standard using premium Schrödinger software. To explore further, molecular dynamics studies on different parameters were performed using Groningen Machine for Chemicals Simulations software to check the stability of the ligand-protein complex and to observe the h-bonding interactions in motion. We expect that this study will provide a useful approach for anticancer drug discovery.

MATERIALS AND METHODS
In this present study, the initial binding study was carried out using molecular docking at standard precision (SP) using Glide (Schrödinger Suites 2017-2, LLC, Cam-bridge, USA). The hardware used for this study was from DELL brand with Intel Core i7-6700 clocked at 3.40 GHz X 8 processor with 8 GB RAM running on 64-bit CENT OS Linux operating system.

Preparation of enzyme
The 3-dimensional (3D) crystallographic model of EGFR tyrosine kinase [Protein Data Bank (PDB) ID: 1M17] was retrieved from Protein Data Bank (www.rcsb.org) using the Protein preparation wizard in Maestro. The protein was preprocessed by performing assign bond orders, removing original hydrogens, adding hydrogens, creating zero-order bonds to metals, creating disulfide bonds, converting selenomethionines to methionines, filling in missing side chains using prime, filling in missing loops  using prime, deleting water molecules beyond 5 Å from hetero groups, and generating het states using Epik at pH: 7.0 ± 2.0. Then, the protein was refined by doing H-bond assignment by sample water orientations, using PROPKA pH: 7.0, and optimization was carried out. The restrained protein minimization was performed using OPLS3 force field.

Receptor grid generation
Receptor grid of size 20 Å was generated after preparing the protein and by selecting the co-crystallized ligand Uracil-Coumarin hybrid inside the cavity of EGFR tyrosine kinase as the centroid of the active site.

Preparation of ligand
All the molecules were built, 3D optimized, and energy minimized using ligprep and confgen function in Maestro. The library was saved in .maegz format.

Molecular docking
Subsequently, molecular docking was performed using Glide at the SP algorithm to evaluate the protein-ligand interaction and their binding affinities. For our analysis, we checked to add Epik state penalties to docking score and enhanced planarity of conjugated pi groups. The validation for molecular docking was done by re-docking the co-crystalized ligand Uracil-Coumarin hybrid in the grid. The resulting poses and interactions of co-crystalized ligand and library were compared.

Molecular dynamics simulation
Based on docking results, the lowest energy and best-posed complex of molecules A6 were selected for the simulation study using Groningen Machine for Chemicals Simulations 2019 package (Pronk et al., 2013) with OPLS2015 force field. Ligand topology files were generated using PRODRG (Schuttelkopf and van Aalten, 2004). The charge of the system was neutralized by the addition of the ions. The energy minimization of the complex (nsteps ¼ 50,000) was executed using the steepest descent approach (1,000 ps). The particle mesh Ewald method was employed for energy calculation and for electrostatic and Van der Waals interactions. The linear constraint solver (Amiri et al., 2007) algorithm was used for covalent bond constraints. Finally, 50 ns molecular dynamics simulation (MSD) with periodic boundary conditions was conducted for respective protein-ligand complexes with nsteps 5,00,000. The root-mean-square deviation and fluctuation (RMSD/F), intramolecular hydrogen bonds, radius of gyration (ROG) (Rg), and thermodynamic parameters were analyzed using Xmgrace (http://plasma-gate.weizmann.ac.il/ Grace/) and VMD 1.9.2 software.

Molecular docking studies
Docking studies against the target protein EGFR (PDB ID: 1M17) that involved in cancer cell proliferation and growth were carried out for 28 novel designed Uracil-Coumarin hybrid compounds (parent structure is given in the following) to check the binding affinity of ligands with the receptors. To validate the applied docking method, corresponding ligands were docked into the tyrosine kinase domain of EGFR using Schrödinger software. The protein structures presented were natural substrates as a co-crystallized ligand, whereas in others, the co-crystallized ligand was a known EGFR inhibitor (i.e., Erlotinib) took as the reference drug (Fig. 3); in both cases, the same docking and scoring validation process were used. Each co-crystallized ligand was previously removed from the respective protein binding site. The predicted docking poses 2D and 3D were compared with the experimental co-crystallized binding pose.
All the data and docking gScore (Glide score is related to ligand binding free energy) and Emodel score (related to electrostatic and van der Waals energies) are provided in Table 1. The hydrogen-bonding interactions, pi-pi stacking interactions, and pi-cation stacking interactions of ligands with the target protein residues were also mentioned. Amid 28 designed hybrid compounds, almost all of them exhibit good interactions with the target protein EGFR tyrosine kinase.
In comparison to standard Erlotinib, compounds A6, A2, and A7 (Fig. 4) were found to show better interactions with target and having a docking score of −8.891, −8.712, −8.658 kcal/mol. The hybrid compound A6 exhibits the most effective binding with target protein EGFR tyrosine kinase (PDB ID: 1M17). In terms of affinity, A6 forms four different hydrogen-bonding interactions; out of them, three H-bonds were with Cys773, Met769, Lys721 amino acid residues and one with water (H 2 O) molecule with gScore −8.891 kcal/mol and Emodel score −100.744. On the other side, reference drug Erlotinib forms three hydrogen bonding interactions; two of them were with Cys773, Met769 and one with (H 2 O) molecule with gScore −8.538 kcal/mol and Emodel score −80.588.
Virtual 2D and 3D screening comparison between most active hybrid ligand (A6) and standard Erlotinib (Figs. 5 and 6) clearly indicates that A6 forms similar kind of interactions and same mechanistic inhibition as of standard drug follows. One additional hydrogen bonding with Lys721 amino acid residue, extra nitrogen heteroatom, and flexibility provided by the linker makes it a superior ligand for the target EGFR.
Various individual hydrogen bonding interactions with virtual bond length in Å of compound A6 with target protein EGFR tyrosine kinase are depicted in Figure 7, and interactions of reference Erlotinib are also shown in Figure 8 for more focused comparison. Referring to Table 1, hybrid compounds A2 and A7 were also reported to have good gScore and Emodel score as compared to the reference drugs. The compound A2 with protein EGFR forms four hydrogen-bonding interactions with amino acid residues Cys773, Lys721,Asp831 and one interaction with water H 2 O molecule. All the interactions of A2 with the target EGFR are shown in Figure 9 in both 2D and 3D dimensions.
Hybrid compound A7 was also found to have hydrogen bonding interactions with Cys773, Lys721, and water (H 2 O). Moreover, compound A7 also has pi-pi stacking and pi-cation stacking interactions with Phe699 and Lys721 residues, respectively. All the interactions are collectively mentioned in Figure 10, whereas additional interactions make it an effective binder in comparison to standard Erlotinib. Hence, on the basis of docking scores of hybrid compounds in Table 1, it can be suggested that among the library of hybrid ligands, Uracil-Coumarin hybrids having two carbon chain length (n = 1) with halogen substituent(X = −I, −Cl, −F) at Uracil moiety A6, A2, and A7 were found to have the best affinity for target EGFR. Although A3 with bromine substituent (X = Br) is found to form prominent interactions with target with comparable docking score (−8.29 kcal/mol) to the reference Erlotinib (−8.538 kcal/mol). Being in the same category, compound A5 showed interaction with poor binding score, which justifies its least affinity toward EGFR and ineffective binding too. Compounds A1, A4, and A5 with substituents other than halogens (X = −H -CH 3 , and −NO 2 ) were found to have a lower affinity and poor docking scores as well.

Molecular dynamic simulations study
MSD is very important to understand structure-tofunction relationships in motion and other adaptive conformational changes by the proteins (Gelpí et al., 2015). Complying to this, an MSD study was performed to evaluate the dynamic behavior of the top-scoring ligand A6 in complex with EGFR tyrosine kinase (PDB ID: 1M17) for a period of 50 ns. The study was analyzed on the basis of parameters such as protein-ligand RMSD, root-meansquare fluctuation (RMSF), ROG, thermodynamic parameters, and protein-ligand interactions formed.

RMSD and RMSF analysis
Analyzing the protein RMSD can give insights into a structural conformation that protein undergoes throughout the simulation. The plots for protein Cα versus time (50 ns) for the simulation are shown in Figure 11, which shows that complex attained an equilibrium value of 0.85 nm at around 30 ns. At Initial stage, a sharp increase in the deviation value from 0.2 to 0.8 nm occurs in the first 3 ns pointing toward a large conformational change in the protein during that period.
Analyzing the ligand backbone RMSD from the plot in Figure 12 indicates how stable the ligand is with respect to the protein and its binding pocket. The plot for ligand RMSD (nm)  versus time (50 ns) shows a plateau value of 0.3 nm for the ligand A6 at around 10 ns. The deviation value observed was within limits, which is quite stable for binding of the ligand with the protein.
The protein RMSF is useful for characterizing local changes along the protein chain. The plot for protein RMSF (nm) versus residue number index shown in Figure 13 describes fluctuations in the range of 0.1-1.4 nm for the complex. The fluctuations in both the terminal ends were higher. The catalytic residue values are no higher than 0.3 nm and are suitable for stable binding.
The ligand RMSF is useful for characterizing adaptive changes in the ligand atom positions (Dong et al., 2018) and can give insights into how ligand fragments interact with the protein and their entropic role in the binding event. Ligand with their atomic index number is shown in Figure 14.
The plot for ligand RMSF (nm) versus ligand atom index shown in Figure 15 describes fluctuations in the range of 0.05-0.33 nm for ligand A6. The central Uracil fragment acts as the core and holds the two triazole-linked coumarins like two arms attached to the core. The Uracil fragment has the least fluctuation which is nearly 0.1 nm. The triazoles attached tend to fluctuate between 0.1 and 0.2 nm, and the coumarins present at the extended length tend to fluctuate above 0.2 nm. The higher value of RMSF (nm) for coumarin fragments may be due to the three atom long linker that attaches them to triazoles.

H-bond analysis
The H-bond interactions formed during the Molecular Dynamics simulations up to 50 ns were also examined. The plot in Figure 16 represents the number of H-bond contacts formed between the ligands and the protein. Ligand A6 shows one stable H-bond contact during the simulation, although occasionally two H-bond contacts, during the period of 26-42 ns, were also formed. From the discussion, it can be summarised that A6 has an efficient binding inside the cavity of EGFR tyrosine kinase.
The H-bonding interaction fraction shown in Figure 17 conveys that stable H-bonding occurs with residue MET769 for around 38% of duration, followed by supplementary H-bond interactions with residues ARG817, LYS704, and TYR703 with interaction for a duration of 8.34%, 6.30%, and 3.90%, respectively.
Results also suggested that irrespective of substituent X on Uracil moiety, the carbon chain linker (n = 1 − 4) played an important role in making an effective binding pose between ligand and the target. Table 1 clearly suggested that as the carbon linker chain increase gScore of the ligand varies continuously from stable −8.891 kcal/mol to unstable −5.53 kcal/mol. The variation in docking scores clearly indicates     that as the linker chain length increases, the flexibility of the molecule increases too, which makes the molecule either unfit to bind in the tyrosine kinase domain of the EGFR or the ligand is unable to make an effective binding pose with protein residues. Evidently, molecular dynamics studies of most effective ligands A6 provided additional proof of forming a stable complex with target EGFR. Ligand A6 forms plateau phase within 10 ns, which retains affirmative up to 50 ns time fraction with little RMSD of 0.3 nm. Both the studies suggested that bigger linker chain hybrids have the property of flexibility and their larger size keep the ligand out of range of interactions, which led to weaker binding interactions and poor docking score. Contrary to this, hybrid molecules with the shortest linker form more     affirmative binding probably due to their less flexibility, higher residue interactions in comparison to standard and their fit size and fit size according to the pocket of the kinase domain of EGFR were proved to be the best ligands.

CONCLUSION
This study emphasizes on the potential interactions between 1,2,3-triazole tethered Uracil-Coumarin hybrids and target protein EGFR tyrosine kinase. Docking results clarified that all hybrid ligands have good H-bonding interactions with the target and few of them also showed pi-pi stacking interactions too. It was also deduced from the results that short carbon chain linker (n = 1) between triazole and Coumarin ring and a halogen substituent (X = −I, −Cl, −F, −Br) at Uracil moiety are two essential parameters for effective binding and to fit properly in the active tyrosine kinase pocket of EGFR. In the library of designed hybrids (A1-D7), compounds A6, A2, and A7 are reported to have good binding interactions as compared to standard Erlotinib. Hybrid compound A6 was found to have the highest binding affinity and form a perfect binding pose with EGFR, justified by gScore and Emodel score. Furthermore, the molecular dynamic study of the most effective ligand A6 has proved the stability of ligand-protein complex  with little deviation and fluctuation up to 50 ns. Hence, it can be concluded that hybrid ligand A6 can serve as the lead compound in the designing of a novel anticancer molecule. Hopefully, the results may help the medicinal chemist in designing and development of a new hybrid drug having a better anti-cancer profile.