{'fraction_non_alphanumeric': 0.028318727685816295, 'fraction_numerical': 0.027669587796170073, 'mean_word_length': 4.831322450910811, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational modeling of novel quinazoline derivatives as potent epidermal growth factor receptor inhibitors

QSAR modelling on Thirty (34) novel quinazoline derivatives (EGFRWT inhibitors) as non-small cell lung cancer (NSCLC) agents was performed to develop a model with good predictive power that can predict the activities of newly designed compounds that have not been synthesised. The EGFRWT inhibitors were optimized at B3LYP/6-31G* level of theory using Density Functional Theory (DFT) method. Multi-Linear Regression using Genetic Function Approximation (GFA) method was adopted in building the models. The best one among the models built was selected and reported because it was found to have passed the minimum requirement for the assessment of QSAR models with the following assessment parameters: R2 of 0.965901, R2adj of 0.893733, Qcv2 of 0.940744, R2test of 0.818991 and LOF of 0.076739. The high predicted power, reliability, robustness of the reported model was verified further by subjecting it to other assessments such VIF, Y-scrambling test and applicability domain. Molecular docking was also employed to elucidate the binding mode of some selected EGFRWT inhibitors against EGFR receptor (4ZAU) and found that molecule 17 have the highest binding affinity of -9.5 kcal/mol. It was observed that the ligand interacted with the receptor via hydrogen bond, hydrophobic bond, halogen bond, electrostatic bond and others which might me the reason why it has the highest binding affinity. Also, the ADME properties of these selected molecules were predicted and only one molecule (34) was found not orally bioavailable because it violated more than the permissible limit set by Lipinski's rule of five filters. This findings proposed a guidance for designing new potents EGFRWT inhibitors against their target enzyme.


Introduction
Among the foremost cancer issues that results in loss of lives each year in the globe which was estimated for almost one-third of the entire cancer deaths is lung cancer. Non-small cell lung cancer (NSCLC) is the principal subset of lung cancers that estimates for about 85% of the problems raised above [1]. The most common cause of NSCLCs recognised was EGFR kinase. It was found in about 10-15% and 30-40% of the population of patients in Caucasia and Asia. It mostly affects women and cigarette smokers in general [1].
Development of inhibitors for mutant-selective kinase is among the difficulties faced in medicinal chemistry and is the principal concern for EGFR tyrosine kinase inhibitors [2]. The remedy of epidermal growth factor receptor (EGFR) to managed non-small cell lung cancers with the T790M resistance mutation prevails a vital medical necessity [3].
In patients with stimulating modifications of EGFR, EGFR inhibitors show a very high response rate. EGFR inhibitors are categorised into two classes: First generation EGFR inhibitors which are referred to reversible inhibitors and include gefitinib and erlotinib. The second class which consist the second and third generation EGFR inhibitors. The second and third generation EGFR inhibitors are referred to as irreversible inhibitors (examples are afatinib and osimertinib). All these classes of drugs were designed to mitigate the problem of NSCLC most especially the EGFR L858R mutations (First generation EGFR inhibitors were designed to treat this type of mutation), EGFR T790M mutations (while Second generation EGFR inhibitors were designed for the treatment of this type of mutation) and EGFR T790M/L790M double mutations (third generation EGFR inhibitors) were designed to treat this type of mutation [2,4,5,6].
QSAR modeling is a molecular modeling method which quantitatively correlate response variable (biological activities) and molecular descriptors (physicochemical properties) of a molecule [7]. In addition, the QSAR technique of computer-aided drug design plays a significant role in predicting the biological activities of small molecules that have not been synthesised [8]. Another virtual screening method applied in computer aided drug design is molecular docking which give an overview of how the ligand and the receptor interact with one another using their individual 3D structures [9]. To have an insight on how body response to the administration of drugs there is need to study the ADME and drug likeness of the drugs before it reaches the final (clinical) stage [10].
The aim of this work is to develop a model with good predictive power which could be used to predict the inhibitory activities of newly designed compounds using QSAR technique, study the mode of binding interactions between some selected EGFR WT Inhibitors and EGFR enzyme via docking and also to predict the ADME properties of these selected EGFR WT Inhibitors.

Dataset source
Thirty four (34) quinazoline derivatives bearing various 6-benzamide moieties as potent EGFR WT inhibitors with their inhibitory activities (IC 50 ) in nM were selected from the work of Hou et al., for this research [11]. The inhibitory activities (IC 50 ) of all the dataset were then converted to their corresponding negative logarithms (pIC 50 ) using Eq. (1) [12]. Table 1 presents the structures, IC 50 and pIC 50 for all the data set  used in this research. pIC 50 ¼ -log IC 50 Â 10 À9 (1)

Sketching of structures and optimum structure generations
After data collection, the sketching of the 2D structures of the studied molecules was achieved using Chemdraw software version 12.0.2 [13]. After sketching the 2D-structures of the dataset, Spartan 14 software was used to convert the 2D-structures to 3D-strucutres before energy minimization. Energy minimizing was performed to reduce constrain in the structures before geometry optimization. Geometry optimization is a process of finding the most optimum structure of a molecule on potential energy surface and this was performed by utilizing Spartan 14 software. DFT at B3LYP/6-311G* level of theory was used in finding the most optimum structures of all the studied molecules on global minima on the potential energy surface (PES) [14].

Descriptors computation, data pre-treatment and daataset splitting
In order to compute the independent variables (descriptors), the most optimum structures obtained in 2.2 above were saved in SDF a file format that is been recognized by the software used in computing the descriptors, PaDEL descriptor tool kit. PaDEL descriptor tool kit was used to compute both Fragment count descriptors, Topological descriptors and Geometrical descriptors [15].
To eliminate redundant and constant descriptors, data pre-treatment was performed manually in this regard.
After pre-treating the data, Data division software was further used in splitting the data into model building and validation set utilizing Kennard-Stone algorithm [16]. The model building which comprise 24 molecules (70%) were used for the generation of the models as the name implies and the validation set which contain 10 molecules (30%) were used for the assessment of the generated models [17].

Building of the model
The models were built utilizing Genetic Function Approximation (GFA) method with the descriptors as independent parameter and the actual pIC 50 as the response parameter. In the case of variable selection, the GFA creates an original population of descriptor sets and determines the most suitable set from it by utilizing evolutionary crossover and mutation speculators which generates a succeeding derivative population of descriptor sets. GFA select most highly correlated descriptors to develop so many models which is one of the distinct characteristic of GFA [18]. The MLR-GFA equation for the model is shown below: where X's are the descriptors, y's are the co-efficient of the corresponding descriptors and z is the regression constant.

Assessment of the model built
The assessment parameters used in evaluating or validating the quality of a QSAR model are the; Squared correlation coefficient of the training set (R 2 training ), Adjusted R 2 (R 2 adj ), Cross-validation coefficient (Q cv 2 ), and Squared correlation coefficient of the test set (R 2 test ) [19,20,21].
The large value of these parameters seem to be important but not enough [22]. In this regard, the multi-collinearity between descriptors can be identified using their variation inflation factors (VIF), to identify whether these descriptors correlate with each other or not. If the estimated VIF values are equal to 1 it means there is no correlation between them; if it happens to be between 1-5, there is high chance of accepting the model; and if it is greater than 10, the model cannot be accepted is therefore rejected [23]. It can be calculated using the equation below: The assessment of importance and participation of each descriptor to the selected model is made using the value of the mean effect (ME) of each descriptor. The equation used in calculating the ME is shown below: where ME represents the mean effect of a descriptor j in a model, the coefficient of the descriptor J is represented by β j in the model and the value of the descriptor in the data matrix for each molecule in the model building set is d ij , n is the number of molecules in the model building set and m is the number of descriptor that appear in the model [24]. Y-Scrambling test was performed to assure the robustness of a model and also the model was not achieved by chance correlation. It is done by reshuffling the actual activities and holding the descriptors fixed to generate new QSAR models for many trials, the new built QSAR models were anticipated to give low Q 2 and R 2 value. The validation parameter for this test is cRp (cR 2 p > 0.5) [25].

Applicability domain
The applicability domain (AD) of a model was carried out to determine whether a model can be regarded valid and void if the model can make a good prediction of new activities of the training and test molecules. As such, the model is subjected to AD to find out whether there are influential or outliers molecules in the studied ones [26]. Leverage approach is among the methods used in assessing the AD of QSAR models and thus is given as h i : where the model building set matrix I is given by x i , n Â k descriptor matrix of the model building set is represented by X and X T is the transpose matrix X used in generating the model. The thresh-hold for the value of X is the warning threshold (h*) which is presented in the equation below: where the number of chemicals of the model building set is given by q, and the number of the descriptors in the model under evaluation is represented by x.

Molecular docking analysis
To elucidate the mode of binding interactions between the active site of EGFR enzyme and some selected EGFR WT inhibitors (ligands), A Dell Latitude E6520 computer system, with the following specification: Intel ® Core™ i7 Dual CPU,M330 @2.75 GHz 2.75GHz, 8GB of RAM was  utilized with the help of Pyrex virtual screening software, Chimera, PyMOL and Discovery studio.

Ligands and EGFR enzyme preparation for the molecular docking computational analysis
The first thing to do in any molecular docking analysis is ligands preparation. The preparation of the ligands was adopted from the optimized structures in 2.2 above saved in pdb file format using Spartan'14 wave software [27]. The next thing ought to be done is the retrieval of 3D structure of the EGFR enzyme to be used in this study. The EGFR enzyme with pdb code: 4zau was downloaded from the Protein Data Bank (RSCPDB). Discovery Studio Visualizer was utilized in preparing the EGFR enzyme for the docking analysis, in the course of the preparation, hydrogen was added, water molecule, heteroatoms and co-ligands present on the crystal structure were completely eliminated and saved in pdb file.

Execution of the molecular docking computational analysis
Autodock vina of Pyrex software was used for the docking of the ligands to the active site of EGFR enzyme (pdb ID: 4zau) [28]. Re-coupling of the ligand-receptor (complexes) for further investigation was done with the help of Chimera software [29]. The elucidation of the binding mode interactions of the complexes was achieved using PyMOL and Discovery studio visualizer [30,31].

ADME and drug-likeness properties prediction
SwissADME a free online web tool used in evaluating ADME and drug-likeness properties of small molecules was used to predict the ADME and drug-likeness properties of some selected EGFR WT inhibitors among the data set [32]. SMILES is the input file for SwissADME which contains a molecule per line separated by a space with a name (optional). Molecules can be inserted in SMILES format or pasted, or drawn using the molecular sketcher available in the web tool. If the molecule is ready, the calculations can be setup by clicking on the "Run" button [32].
The Lipinski's rule of five filter is very useful at pre-clinical stage of drug discovery which state that if any compound violate more than 2 of these criteria (Molecular weight < 500, Number of hydrogen bond donors 5, Number of hydrogen bond acceptors 10, Calculated Log p 5 and Polar surface area (PSA) <140 Å 2 ), the compound is said to be impermeable or badly absorbed [33].

QSAR modeling
The reported model was observed to have excelled the limit for the evaluation of a good model with the following evaluation parameters: R 2    of 0.965901, R 2 adj of 0.893733, Q cv 2 of 0.940744, R 2 test of 0.818991 and LOF of 0.076739 as reported by [21] (Table 2).
pIC 50 ¼ 1.496601581 * ATSC6m -1.098227666 * ATSC8e þ 0.786519460 * MATS7m þ 0.410257932 * SpMax3_Bhp -1.755609854 * SpMax5_Bhs þ 2.392960763 * maxHBint10 þ 5.590744728 The details of the descriptors in the reported model were presented in Table 3. The negative coefficients of these descriptors (ATSC8e and SpMax5_Bhs) highlighted their negative correlation to the inhibitory activities of the quinazoline derivatives (EGFR WT inhibitors). The lesser the number of these descriptors in the structures of these EGFR WT inhibitors, the higher the potency of these EGFR WT inhibitors toward their target EGFR WT enzyme. On the other hand, the positive co-efficient of ATSC6m, MATS7m, SpMax3_Bhp and maxHBint10 descriptors in the reported model gives the positive correlation of these descriptors to the inhibitory activities of EGFR WT inhibitors. That is, the more the presence of these types of descriptors in the structures of these EGFR WT inhibitors the more the inhibitory activities of the EGFR WT inhibitors toward their target enzyme.   3.1.1. Description of the descriptors that appear in the reported model ATSC6m and ATSC8e are Moreau-Broto autocorrelation of a Topological Structure, ATS (The ATS descriptor is a graph invariant describing how the property considered is distributed along the topological structure). These descriptors can be seen as a special case in which other types of descriptors can also be derived from [34]. This is the most known spatial autocorrelation defined on a molecular graph G as MATS7m is a Moran autocorrelation which if applied to a molecular graph. Moran coefficient usually takes value in the interval [-1, þ1]. Positive autocorrelation corresponds to positive values of the coefficient whereas negative autocorrelation produces negative values [34]. It can be defined as SpMax3_Bhp and SpMax5_Bhs are the maximum absolute eigenvalue of Burden modified matrix -n 3/and -n 5/weighted by relative Istate and relative polarizabilities, called leading eigenvalue or spectral radius, SpMaxA is the maximum absolute value of the spectrum. These kinds of functions were called by Ivanciuc matrix spectrum operators. This eigenvalue has been suggested as an index of molecular branching, this descriptor talks about branching in molecules [35]. As seen from the regression equation and ME values (Table 5), SpMax5_Bhs contributes negatively to the inhibitory activities of the studied molecules. It suggests that reducing chain branching in the studied molecules will improve the inhibitory activities of the studied molecules toward their target enzyme.
MaxHBint10 is a maximum E-State descriptors of strength for potential hydrogen bonds of path length 10. Based on the regression equation and ME values (Table 5), this descriptor gave the highest contribution toward the inhibitory activities of the studied molecules. Increasing the number of hydrogen bond in the molecules might increase their potency against their target protein.
The XY (Scatter) plot of predicted activities of both the test and training sets against the actual pIC 50 was shown in Figure 1A & 1B. The significance of the reported model was confirmed by the distribution of the values around the straight line. Also, the R 2 values from the plots agree with those of the training and test for the internal and external assessment.
On the other hand, the XY (Scatter) plot of actual pIC 50 against the residuals of both the model building and validation sets was shown in Figure 2. The unusual occurrence of these residuals on either side of zero on the plot shows the non-existence of methodological error in the reported model.
The pIC 50 , Predicted pIC 50 and the residual values for all the studied molecules were presented in Table 4. The low residual values noted in the table verified the reliability of the reported model.
The correlation statistical analysis of the descriptors in the reported model was performed (Table 5) and the descriptors were found not to correlate with one another. This shows the high performance of the descriptors utilized in generating the reported model. To further confirm whether there is a similarity or not between the descriptors in the reported model, The VIF values of these descriptors in the model building set were estimated and realized to be less than 2 (Table 5) indicating the applicability of the reported model and thus the descriptors were independent of one another. The ME value (Table 5) gives the contribution of a descriptor in opposition to other ones in the reported model. The signs point the various directions of either increase or decrease in the values of these descriptors which will improve the inhibitory activities of the studied molecules. It is observed that from the model and ME values (Table 5), maxHBint10 descriptor gives the highest contribution.
The Y-scrambling test was presented in Table 6 for the 10 randomly generated models and the R 2 and Q 2 values for the newly generated random models were determined to be very low. This has affirmed the obtainability of the reported model was not by chance and further confirm its robustness.
The plot of leverages against standardized residuals of both the model building and validation sets (Williams plot) presented in Figure 3 identified two (2) influential compounds from which were all in the validation set. It is very paramount to decipher that these molecules (influential compounds) with leverage value greater than the threshold h*(h* ¼ 0.875) are not put into consideration when designing new EGFR WT inhibitors. These molecules might be structurally different from those used to generate the reported model and, thus may have a different mechanism of action.

Molecular docking analysis
The mode of binding interactions between the active site of EGFR receptor (4zau) and some selected EGFR WT inhibitors (ligands) was elucidated through molecular docking (Table 7). From Table 7, Complex 17 was identified to have the highest binding affinity of -9.5 kcal/ mol. With the help of discovery studio visualizer, the ligand was clearly observed to have interacted with the active site of EGFR receptor via Hydrogen bond with the following amino acid residues MET793, MET793, THR854 and ASP855 with bond distances of 2.61394 (Å), 2.18464 (Å), 2.57601 (Å) and 2.68794 (Å). The interaction was not only via hydrogen bond, it also interacted with the active site of the EGFR receptor via halogen bond (GLN791), hydrophobic bond (LEU718, CYS797, LYS745, ALA743, ALA743, and VAL726), electrostatic bond (LYS745) and others (MET766) which might be the reason why it has the highest binding affinity. The next one identified with good binding affinity after the one mentioned above (complex 17) is complex 34. It interacted in the active site of the receptor through hydrogen bond with ARG841, ASN842, LYS745, UNK1 residues with bond distances of 2.7122 (Å), 2.07811 (Å), 2.54982 (Å), and 3.62317 (Å). It also interacted with the active site of the EGFR receptor through halogen bond (LEU788), hydrophobic bond (VAL726, LYS745, LEU718, ALA743 and LEU844), electrostatic bond (ASP855) and others (CYS797). The rest other three complexes interacted in the active site of the receptor through hydrogen bond, halogen bond, electrostatic and hydrophobic bond as shown in Table 7. Figure 4 & Figure 5 showed the 3D and 2D structures of the complexes. Based on the molecular docking results, the most common amino acids to all of the examined compounds were MET793, LEU718, LYS745 and VAL726 ( Table 7). The most important amino acids that might be responsible for the higher binding affinity were CYS797, GLN 791 and ASP855 due to their interaction with the molecules with higher binding affinity (Table 7). Halo substituted molecules (Complex 17) were found to fit better in the active site of the receptor than those with bulkier substituents (Complex 34) as shown in Figure 5.
On comparing the QSAR and docking results, the molecule with the highest activity was among those having higher binding affinity. This means that there is little correlation between the QSAR and the molecular docking studies.

ADME properties prediction
The ADME properties of these selected EGFR WT inhibitors were predicted and presented in Table 8. From Table 8, it can be observed that only one among these molecules violated more than the maximum permissible limit of the criteria stated by Lipinski's rule of five, it means there is a high tendency all of these molecules might be pharmacologically active except molecule 34 which has more than 3 violation. In a null shell the remaining four (4) molecules are said have good absorption, low toxicity level, orally bioavailable and permeable properties. The Bioavailability Radar gives an overview of the drug-likeness of all the selected molecule (Figure 6.). The painted pink area shows the range for each properties (Lipophilicity: XLOGP3 between À0.7 and þ5.0, size: MW between 150 and 500 g/mol, polarity: TPSA between 20 and 130 A 2 , solubility: log S not higher than 6, saturation: fraction of carbons in the sp3 hybridization not less than 0.25, and flexibility: no more than 9 rotatable bonds). Based on this criteria, all the molecules are said to be orally bioavailable except molecule 34 which is too Flexible, Polar, Lipophilic and Insoluble. The plot of WLOGP against TPSA (Boiled-egg plot) to predict gastrointestinal absorption and brain penetration of the selected molecules was shown in Figure 7. It can be seen from the plot that none of the molecules possess the BBB permeant but they are within the GI absorption region.

Conclusion
A very high predictive model was developed using QSAR modelling technique on some EGFR WT inhibitors. The reported model was selected and reported because of its fitness with the following assessment parameters: R 2 trng ¼ 0.919035, R 2 adj ¼ 0.893733, Q 2 cv ¼ 0.866475, R 2 test ¼ 0.636217, and LOF ¼ 0.215884. The high predicted power, reliability, robustness of the reported model was verified by other assessments such as AD, Y-scrambling test and found to be statistically fit. The molecular docking results of the examined compounds showed that CYS797, GLN 791 and ASP855 amino acids might be responsible for the higher binding affinity of molecule 17 (-9.5 kcal/mol) and 34 (-9.4 kcal/mol). Also, the   results of ADME properties predicted indicated that only molecules 34 among others was nor orally bioavailable as it has violated more than the maximum permissible limit for the orally bioavailability of drugs set by Lipinski's rule of five. This research was able to identified compound 17 as a lead among the studied compounds and proposed when designing new compounds it should be used as template for structural modification.

Author contribution statement
Adamu Uzairu: Conceived and designed the experiments. Muhammad Tukur Ibrahim: Performed the experiments; Wrote the paper.
Sani Uba, Gideon Adamu Shallangwa: Analyzed and interpreted the data.

Funding statement
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.