In-Silico Evaluation, Chemical Reactivity, and Covalent Docking Study of Various Quinazolines and Pyridopyrimidines as Inhibitors for the Epidermal Growth Factor Receptor

Abstract The epidermal growth factor receptor (EGFR) sensitivity to tyrosine kinase inhibitors gefitinib and erlotinib-activated mutants is common in non-small cell lung cancer. An in-silico screening study of quinazolines and pyridopyrimidines against the wild-type (WT) and mutation of EGFR tyrosine kinase inhibitors was conducted, employing several computational approaches such as covalent docking and molecular dynamics simulation followed by reactivity and the absorption, distribution, metabolism, excretion, and toxicity (ADMET). Twenty-two heterocyclic compounds were screened by covalent docking against WT and mutated proteins (L858R and T790M). Compounds 4 and 19, which contain quinazoline[3,4-d]pyrimidine and pyrido[3,4-d]pyrimidine, respectively, were found to have an affinity toward the wild-type and the mutant protein. In addition, they had good chemical reactivity and kinetic stability toward the WT and mutations and desirable ADMET properties. These findings reveal new, robust, and irreversible tyrosine kinase inhibitors for the WT and its mutant proteins. Graphical Abstract


Introduction
2][3] The EGFR function is to enhance cell growth and division.In the case of mutations, the EGFR protein would not function properly.The uncontrolled activation of EGFR leads to proliferation, differentiation, migration, and angiogenesis, which are associated with human cancers, such as breast cancer, non-small cell lung cancer (NSCLC), and glioblastomas. 4EGFR deregulation may be caused by activating mutations. 5,6The tyrosine kinase (TK) domain of EGFR is encoded by six exons.EGFR kinase activating mutations occur in exons 18-21 and are classified into three categories, the deletion of exon-19, the substitution of a single nucleotide that changes the coded amino acid, and in-frame duplication and insertion of exon 20. 6 The first category includes deletions of leucine-747 to glutamic acid-749 and in the second category leucine to arginine at 858 (L858R). 7EGFR-TK mutant responds to anilinoquinazoline based on the structures of gefitinib and lapatinib inhibitors. 8EGFR wild-type (WT) inhibitors canertinib, 9 afatinib, 10 pelitinib, and neratinib 11 (Figure 1) exhibited significant clinical response in patients with NSCLC and increased the sensitivity of NSCLC cells. 12,13The third category involves changing threonine 790 into methionine. 14This mutation occurs in a hydrophobic pocket at the back of the ATP binding pocket and increases the resistance of first and second-generation EGFR tyrosine kinase inhibitors. 15,16The sensitivity of first-and second-generation inhibitors against EGFR WT, L858R, and T790M mutants enables the development of new inhibitors to block these mutations in cancer cells. 17he computational approach has become a useful tool in drug discovery and development.][20] Covalent drugs typically have irreversible interactions with their targets; hence, they are more potent while maintaining a small molecular size, which is an advantage in the lead optimization step.Covalent engagement with the target protein provides the benefit of a longer-lasting biological effect. 21Electrophilic warheads found in covalent inhibitors (eg, epoxides, aziridines, esters, ketones, a,b-unsaturated carbonyl compounds, acetylenes, and nitriles) react with nucleophilic residues in the active site of proteins (eg, cysteine and threonine) and create a covalent bond. 22,23Frequently, attaching the original non-covalent scaffold to a homology model is sufficient.Mostly to enhance a lead covalent inhibitor to discover more powerful molecules in which structure-activity relationship studies can be used.5][26][27][28] Computational approaches methods are being used in many different fields. 29,30Also they could accelerate the development of new tyrosine kinase inhibitors.
2][33] This study attempts to identify a lead compound containing quinazoline or pyrido [3,4-d]pyrimidine moiety for EGFR wild-type and its mutated isoform (L858R and T790M) by various computational techniques such as reactivity, covalent docking, MD simulation, and ADMET.

Molecule library preparation
A dataset of substituted quinazolines and pyrido [3,4-d]pyrimidines 1-22 as irreversible tyrosine kinase inhibitors of the EGFR family was extracted from the literature. 34The IC 50 of selected compounds 1-22, with values (0.002-0.026 mM), are summarized in Figure 2. The structures were optimized using a density functional theory (DFT) method with the B3LYP/6-31G basis set. 35,36to obtain the most stable conformation, which was also used to calculate the global reactivity descriptors through the Gaussian 09. 37All optimized structures were confirmed to be minima with no negative frequencies. 38The affinity of the optimized structures toward EGFR proteins was investigated using the MOE software. 39

Receptor preparation
The 3D-crystal structures of proteins are presented in Figure 3, which were obtained from the RSCB PDB database, 40 wild-type PDB ID: 1XKK, 41 and for L858R, T790M mutations (mutant protein) PDB ID: 5HG5. 42The enzymes were prepared by removing the cofactors, phosphate ion for the 1XKK, and sulfate ion and glycerol for the 5HG5.4][45] MOE software was used for all the steps of enzyme preparation via the preparation panel, which protonates the charged residues and minimizes the angles along with the dihedral angles. 39The residues of the active sites of each enzyme were obtained using the site finder and are presented in Table 1.Based on the literature, Lys745, Thr790, Met793, and Asp855 residues are the most crucial amino acids for kinase inhibition.While the Lys745, Met793, and Cys797 residues are important for enzyme binding. 46

Covalent docking
The 22 heterocyclic compounds were subjected to covalent docking against the structure of the EGFR WT protein and EGFR mutant protein.The molecular operating environment software performed the covalent docking and scoring calculations (MOE). 39The crystal structure of the Table 1.Binding sites residues that used as input for receptor grid generation during the induced fit docking.

Molecular dynamic simulation
The Schr€ odinger LLC Desmond simulation package was used to run the molecular dynamics simulations. 53Across all runs, the NPT ensemble with a temperature of 300 K and a pressure of 1 bar was used.The simulation's length was 100 ns for all the chosen ligands, and its relaxation time was one ps.All simulations employed the OPLS3 force field settings. 54In Coulomb interactions, the cutoff radius was 9.0 Å.The orthorhombic periodic box boundaries were set 10 Å away from the protein atoms.,The water molecules were explicitly described using the transferable intermolecular potential with three points (TIP3P) model. 55,56Salt concentration was set to 0.15 M NaCl and built using Desmond's System Builder utility. 57The Martyna À Tuckerman À Klein chain coupling scheme with a coupling constant of 2.0 ps was used for the pressure control and the Nos eÀHoover chain coupling scheme for the temperature control. 58,59Nonbonded forces were calculated using a RESPA integrator, where the short-range forces were updated every step, and the long-range details were updated every three steps.The trajectories were saved at 20 ps intervals for analysis.The behavior and interactions between the ligands and protein were analyzed using the Simulation Interaction Diagram tool provided in the Desmond MD package.The stability of MD simulations was monitored by looking at the RMSD of the ligand and protein atom positions in time.

Global reactivity descriptors
The global reactivity indices are derived from the conceptual density functional theory (DFT) 60 and can be used to understand compounds' chemical reactivity and kinetic stability. 61,62The global reactivity descriptors can be described by the energy of the highest occupied molecular orbital (E HOMO ), the power of the lowest unoccupied molecular orbital (E LUMO ), energy gap 3][64] The TCE was taken as a reference for the nucleophilicity value since it presents the lowest EHOMO in a large series of molecules. 65

Computational pharmacokinetics
Most therapeutic compounds failed during the clinical trial due to poor ADMET properties. 66he QSAR models can predict the ADMET properties for drugs at the design stage. 67,68There are useful databases to predict the ADMET properties such as ADMETlab, 69 SwissADME, 68 and admetSAR. 69The ADMETlab online software 58 was used to predict the ADMET properties and provide a more precise prediction than the SwissADME and the admetSAR. 70he ADMET properties were calculated for the selected compounds that resulted from covalent docking as tyrosine kinase inhibitors of the EGFR family.The ADMET was an important step in achieving oral bioavailability.The following parameters were predicated which are important for the absorption (Caco-2 permeability (human colorectal adenocarcinoma cells) > À5.15, Pgp-inhibitor, Pgp-substrate, and HIA (human intestinal absorption)), distribution (PPB: plasma protein binding, VD: volume distribution, and BBB: blood-brain barrier), metabolism CYP450 enzyme (1A2-inhibitor and substrate, 3A4-inhibitor and substrate, 2C9-inhibitor and substrate, 2C19-inhibitor and substrate and 2D6-inhibitor and substrate), excretion (T1/2: half-life and CL: clearance) and for toxicity (hERG, H-HT: human hepatotoxicity, AMES, and LD50.The results of ADMET properties are classified by 1: inhibitor, substrate, or blocker, 0: non-inhibitor, no substrate, or no blocker.

Molecular docking
The binding affinity of heterocyclic compounds 1-22 with the 1XKK (WT) The docking results of 22 compounds containing quinazolines or pyrido [3,4-d]pyrimidines moiety, with the active site of 1XKK protein, are shown in Table 3.The results revealed that heterocycles 4, 6, 18, and 19 have the best RMSD values of 2.01, 2.12, 2.29, and 1.60 Å, respectively.It is worth mentioning that these compounds formed a covalent bond with the Cys797.Meanwhile, compound 18 had the best RMSD value and the best binding score.The residues of the active site that interact with the ligands and their interaction type with the other compounds are reported in Table S1.
The interactions of the most active heterocycles that gave the best docking results are shown in Figure 4. Compounds 4, 6, and 18 have an H-acceptor interaction with methionine amino acid, 71 which is considered an important residue.Several tyrosine kinase inhibitors combined with methionine. 71Also, most heterocycles containing Michael's acceptors interacted with the Cys797 by a covalent bond, while the remaining compounds had a covalent bond with the Thr854 (Figure S1).
The binding affinity of heterocyclic compounds 1-22 with the 5HG5 (mutant protein) Results from docking analysis of 22 compounds quinazolines or pyrido [3,4-d]pyrimidines moiety, with the active site of mutant protein 5HG5, are summarized in Table 4.As can be seen from the results shown in Table 4, compounds 4, 13, 17, and 19 have good energy scores and RMSD values compared with the other ligands.The other type of amino acid residues involved in the protein-ligand interaction, such as hydrogen bonding of other heterocyclic compounds, is summarized in Table S2.
Compound 4 has an p-H interaction with the Leu718 and H-acceptor interaction with Asp800, and compounds 13 and 19 have a covalent bond with the Cys797 residue (Figure 5).Whereas compound 17 has an H-D interaction with the Cys797 residue.The interactions of the remaining compounds are shown in Figure S2.

Analysis of complexes with the 1XKK
To monitor the effect of ligands on the conformational stability of the 1XKK during simulations, RMSD values of Ca atoms were estimated for all the complexes concerning the initial structure (Figure 6).The fluctuation of the proteins was within the acceptable variation with the RMSD (less than 3 Å), indicating the stability of protein conformation.Figure 6 shows that all complexes reach their stable states after 60 ns.The RMSD of the ligands was plotted as a function of simulation time (Figure 7).The RMSD of a ligned and measured ligand was just on its reference conformation within the active site.Figure 7 shows that compounds 4, 6, 18, and 19 were stable within the protein's active site and reached equilibrium in the first ten ns.The average RMSD of these compounds over the last 40 ns were 1.69, 1.02, 2.68, and 1.64 Å, respectively.
The protein-ligand histogram diagram panel, generated with simulation interactions, is implemented in Maestro software.The histogram explains the contacts that occur during the  simulations between ligands and proteins.It can be seen from Figures 8(a) and 9 that hydrogen bonding between 18 and the Met793 and Asp855 residues was maintained most of the time, either directly or through the water bridge H-bond.In addition, hydrophobic binding was observed with the Leu718, Leu792, and Leu844 residues.The Gln791 residue was bound with the ligand through a water bridge (Figure 9). Figure 8(b) illustrates the contacts that occur during the simulations between 19 and protein.Hydrogen bonding between the Asp800 residue and ligand was maintained most of the time, either directly or through the water bridge H-bond, as seen in Figure 9. Also, the ionic binding was observed with the same residue.Other residues, such as Cys797 and Arg841, were able to develop strong hydrophobic interactions most of the time.The interaction histograms and 2D interactions of other ligands are presented in Figure S3.To monitor the protein-18 interactions during simulation, a plot of active site residues was plotted against trajectory frames (Figure 10(a)).Notably, the Met793 and Asp855 were in contact with 18 most of the time.While the Gln791 was in contact with 18 in the first 10 n, then lost contact till 50 ns, and contacted again during the simulation time.In addition, the other residues, such as Leu718, L792, and L844, were in contact with 18 most of the time.
To monitor the protein-19 interactions during simulation, a plot of active site residues was also plotted against trajectory frames (Figure 10(b)).Notably, the Cys797, Asp800, and Arg841 were mostly in contact with compound 19.The Arg803 had a strong interaction with 19 at the beginning of simulation time, and then it was lost at about 20 ns.The interaction trajectory frames of other ligands are presented in Figure S4.

Analysis of complexes with the 5HG5
To monitor the effect of the compounds on the conformational stability of the 5HG5 during simulation, the RMSD values of Ca atoms were estimated for all the complexes concerning the initial structure.Figure 11 shows that all complexes tend to reach their stable states after 10 ns; the fluctuation of the proteins was within acceptable variation (RMSD less than 3 Å), indicating the stability of protein conformation.For the last 20 ns, the RMSD value of Ca for 4/5HG5 was more than 3 Å because of the fluctuation in the C and N-terminal tails of the protein (Figure 12).
Figure 13 shows that heterocyclic 4, 13, and 17 were stable within the protein's active site and reached equilibrium within the first 10 ns.The average RMSD of these compounds were 1.67,  1.28, and 2.53 Å, respectively.Compound 13 showed stability within 10 ns as well and moved by around 1.5 Å from its original position and hold on to the new position toward the end of the simulation.Also, heterocycle 17 showed stability at 10 ns, moved by around 2.5 Å from its original position, and held on to the new position toward the end of the simulation.Finally, 19 showed the lowest stability within the active site of the mutant protein 5HG5.It held its position from the beginning of the calculation till around 40 ns, where it reorients itself inside the active site and moved by approximately 3.5 Å, concerning its initial position.
Figure 14(a) shows the interactions between 17 and protein.Hydrogen bonding between the Asn842 residue and ligand was observed most of the time.Figures 14(a) and 15 show that the Asp837 interacted with compound 17 throughout ionic and hydrogen bonding.Also, the hydrogen bonding between the Lys879 residue and the ligand was maintained most of the time, either directly or through the water bridge.The other residues, such as Phe723, Cys797, Arg841, Phe856, and Arg858, were able to develop strong hydrophobic interactions most of the time.
The hydrogen bonding between the Asp800 and Arg841 residues and ligand 19 was maintained most of the time, either directly or through the water bridge H-bond (Figures 14(b) and 15).The other residues, such as Cys797 and Phe856, were able to develop strong hydro-phobic interactions most of the time.The interaction histograms for the other ligands are presented in Figure S5.
The protein-17 interactions during simulation are shown in Figure 16(a).The Phe723, Cys797, Asp800, Arg841, Asn842, Phe856, Arg858, and Lys879 were in contact with 17 most of the time.The Asp837 has a strong interaction with 19 from 30 ns.
The interactions between protein and 19 during simulation are shown in Figure 16(b).The Cys797, Asp800, Phe856, and Arg841 were in contact with 19 most of the time.The Phe723 residue has a strong interaction with 19 from 40 ns.While the Ala723, Asp837, and Asn842 had a strong interaction with 19 at the beginning of the simulation time, and then it was lost at about 40 ns.The other interaction trajectory frames are presented in Figure S6.

Free binding energy analysis
The free energy of compounds 4, 6, 18, and 19 binding with the 1XKK and 5HG5 were computed using the MM/GBSA approach by the Schrodinger program.The MM/GBSA technique can advise designers of more powerful ligands about important protein-ligand interactions, according to the findings and debate.Table 5 lists the MM/GBSA findings for the investigated ligands interacting with various target proteins.
Free energy calculations are widely used and commonly accepted to estimate the ligandbinding affinities in the protein system.All the binding energy values (DG Bind ) were negative,   indicating that all the compound's hits with protein were favorable except compound 13 with 5HG5, which gave an unexpected positive value.

Global reactivity descriptors
The calculated global reactivity descriptors for the compounds that have the best results from the covalent docking study are summarized in Table 6.Analysis of the DFT descriptors reveals more information about the compounds' stability, electrophilicity, and nucleophilicity under investigation.From the results, E HOMO and E LUMO have negative values, 72 which refer to the strength of investigated complexes.The energy gap (DE) measures a molecule's chemical reactivity and kinetic stability.A large DE means high kinetic stability, low reactivity, and poorly polarizable molecules. 73Compound 18 has the lowest energy gap (DE ¼ 3.49 eV), followed by 4 and 19.This indicates that these compounds have a charge transfer within the molecules and reflect their biological activity. 74,75Additionally, 18 was the hardest (g ¼ 1.75 eV) and softest (S ¼ 0.29) compound, followed by 4 and 19.Compound 18 was the hardest according to the presence of (dimethylamino) at the end of the chain and structure base pyrido [3,4-d]pyrimidine. 76According to the above analysis, compounds 4, 18, and 19 are more stable and reactive than those reinforced in the literature. 77,78ll coordination processes are spontaneous, according to all molecules' negative chemical potential values. 79Compounds 4 have the highest electronic chemical potential (m ¼ À3.55), so it can exchange the electron density with the environment efficiently. 80Organic molecules can be classified as strong (N > 3 eV), moderate (2.0 eV N 3.0 eV), and marginal nucleophiles (N < 2.0 eV) in polar organic reactions. 65All compounds under investigation were classified as moderate nucleophiles except 4, 6, and 13, which were classified as strong nucleophiles.
The electrophilicity (x) provides important information about the reactivity of organic compounds that participate easily in polar reactions and should be more than 2.0 eV. 63,81Compounds 18 (5.33 eV), and 19 (4.80 eV), are more electrophilic than other compounds.In addition, heterocycles 4, 18, and 19 have good DE, g, and S.

Pharmacokinetics properties
The ADMET properties for quinazolines and pyrido [3,4-d]pyrimidines are reported in Table 7.
Compound 18 has a Caco-2 value > À5.15.All compounds are listed as P-gp inhibitors, passed human intestinal absorption HIA, and none were substrate P-gp.Compounds 18 and 19 passed the F 30% bioavailability.All compounds have a PPB of more than 90% (high PPB-bound).All Compounds passed the BBB criteria, and for the VD distribution value: < 0.07 L/kg highly hydrophilic, 0.07-0.70evenly distributed, and > 0.7 highly lipophilic.Compound 18 is highly hydrophilic, and the others are evenly distributed.All compounds were not considered inhibitors for the CYP450-1A2 residue.They were inhibitors and substrates of the CYP450-3A4.The compounds were inhibitors for the CYP450-2C9 except for 18 and 19.Compounds 18 and 19 were substrates of the CYP-2C9.All compounds were inhibitors for the CYP-2C19 except for 17.All compounds were substrates of the CYP450-2C19 except 13.All compounds were inhibitors for the CYP450-2D6, whereas only 4, 13, 17, and 18 were not substrate of the CYP450-2D9.
All compounds have a T1/2 value of > 0. 5h and CL < 5 mL/min/kg as canartinib value and act as hERG blockers.Compound 15 is an H-HT human hepatotoxic.Compounds 4, 11, 13, and 17-19 have no Ames mutagenicity.The LD50 acute toxicity value should be > 500 mg/kg for low toxic, 50-500 mg/kg for toxic, and < 50 for highly toxic compounds.The compounds under investigation have low toxicity.

Conclusions
The inhibitory activity of 22 heterocyclic compounds containing quinazoline or pyrido [3,4-d]pyrimidine ring system against the WT growth factor receptor and its mutant protein was evaluated.The docking assessment with the 1XKK and 5HG5 revealed that heterocyclics 4 and 19 had promising results against the WT growth factor; also, they were the best to interact with the mutant protein.The interactions were further confirmed by molecular dynamics simulation.The ligand-protein complexes were stable for the last 50 ns.Heterocycle 4, which contains quinazoline [3,4-d]pyrimidine moiety, was a good inhibitor for the WT.While compound 19, which contains pyrido [3,4-d]pyrimidine, was a good inhibitor for the WT and the mutant proteins.Heterocycle 19 had good global reactivity descriptors.Heterocycle 4 was considered a strong nucleophile with N ¼ 3.27 eV, while 19 was a moderate nucleophile with N ¼ 2.58 eV.The heterocycles under investigation have passed the ADME property and have no toxicity.The use of a combined computational approach that includes global reactivity descriptors, covalent docking, and molecular dynamic simulation could provide an alternative way to features of the binding mechanism for quinazolines and pyrido [3,4-d]pyrimidines as good inhibitors of growth factor receptor mutation.

Figure 1 .
Figure 1.The chemical structures of first and second-generation tyrosine kinase inhibitors.

Figure 4 .
Figure 4. Compounds binding with the WT epidermal growth factor receptor; PDB ID: 1XKK.

Figure 6 .
Figure 6.Plots of RMSD for Ca atoms (Å) with respect to the initial structure vs simulation time (ns) for All the complexes with the 1XKK.

Figure 7 .
Figure 7. Plots of RMSD for ligand atoms (Å) with respect to the initial structure vs simulation time (ns) for All the complexes with the 1XKK.

Figure 10 .
Figure 10.The (a) 18/1XKK and 19/1XKK interactions shown by the active site amino acids in each trajectory frame.Note: White refers to zero interaction and deep color indicates more interactions.

Figure 11 .
Figure 11.The plots of the RMSD for the Ca atoms (Å) with respect to the initial structure vs simulation time (ns) for the complexes of compounds 4, 13, 17, and 19.

Figure 12 .
Figure 12.The overlap of snapshots of 4/5HG5 complex at 1, 25, 50, 75, and 100 ns shows the fluctuation in the C and N-terminal tails (red).

Figure 13 .
Figure 13.The plots of RMSD for ligand atoms (Å) with respect to the initial structure vs simulation time (ns) for the complexes of compounds 4, 13, 17, and 19.

Figure 16 .
Figure 16.The (a) 17/5HG5 and (b) 19/5HG5 interactions shown by the active site amino acids in each trajectory frame.White refers to zero interaction and deep color indicates more interactions.

Table 2 .
Crystallization, data collection and refinement statistics of the WT and mutant protein.

Table 3 .
The docking results of lead compounds with the 1XKK.

Table 4 .
The docking results of lead compounds with the 5HG5.

Table 5 .
The calculated MM-GBSA binding energies (kcal/mol) for the selected compounds against the 1XKK and 5HG5 proteins over MD simulations.