In silico ADME/T and 3D QSAR analysis of KDR inhibitors

1 Molecular Modeling & Drug Design Laboratory (MMDDL), Pharmacology Research Division, Bangladesh Council of Scientific & Industrial Research (BCSIR), Chittagong-4220, Bangladesh. 2 Industrial Microbiology Research Division, Bangladesh Council of Scientific & Industrial Research (BCSIR), Chittagong-4220, Bangladesh. 3 Institute of Food Science & Technology (IFST), BCSIR, Dhaka. 4 Industrial Botany Research Division, Bangladesh Council of Scientific & Industrial Research (BCSIR), Chittagong-4220, Bangladesh. 5 Department of Microbiology, Faculty of Biological Science, University of Chittagong, Chittagong, Bangladesh. 6 Charles Perkins Center-University of Sydney, NSW 2050, Australia. 7 Department of Nutritional Biochemistry, ICDDR, B, Dhaka 1212, Bangladesh.


INTRODUCTION
New blood vessels formation from existing vasculature is termed as angiogenesis, a crucial pathway for tissue repair, growth and embryogenesis (Cherrington et al., 2000, Ferrara, 2002, Harmange et al., 2008 and associated with different types of pathological processes which include metastasis (Liotta et al., 1991), tumor growth (Folkman, 1972), psoriasis (Detmar, 2000), ocular neovascularization (Aiello et al., 1994), rheumatoid . arthritis (Walsh and Haywood, 2001), inflammation (Fava, 1994). Vascular endothelial growth factor (VEGFa) and its receptor tyrosine kinases VEGFR-2 or kinase insert domain receptor (KDR) and VEGFR-1 were reported as key regulators for angiogenesis (Terman et al., 1992). Briefly, downstream pathways are initiated by the accumulation of receptor dimerization and intercellular autophosphorylation, when VEGF binds with its receptor VEGFRs. As a corollary, vascular permeability, endothelial cell proliferation, survival and migration are increased and resulted in tumor growth and metastasis (Kilari et al., 2013;Quentmeier et al., 2012). Moreover, excessive expression of VEGFR is related to the poor prognosis and aggressiveness of the tumor in most of the cancers and therefore, inhibition of KDR-VEGF binding activity is one of the major strategies in cancer drug development (Li et al., 2016). Some successful strategies to inhibit the angiogenesis have been effectively expressed in preclinical and clinical layout. These approaches include VEGF soluble decoy receptors (Cardones and Banez, 2006), antibodies directed against VEGF (D'Adamo et al., 2005), and small molecules that inhibit KDR (Adjei, 2007;Beebe et al., 2003;Gingrich et al., 2003;Harmange et al., 2008;Hennequin et al., 2002;Ruggeri et al., 2003;Ryan and Wedge, 2005;Sun et al., 2003;Thomas et al., 2003). From these observations and ever increasing demands of these biological active KDR inhibitors, we subjected some novel N-4chlorophenylnaphthamide based KDR inhibitors to study their QSAR and QSTR profiles.
The approach are used based on the popular quantum based drug design methodologies to describe the ADME/T profiles and activity by taking some quantum chemical parameters including dipole moment (D), LUMO energy (E LUMO ), hardness (η) and electrophilicity (χ), along with hydrophobic descriptor SLogP.

Dataset preparation
Thirteen N-4-chlorophenylnaphthamides with KDR inhibitory activity were used in systemic ADME/T and QSAR analysis are represented in Table 1. The molecules were drawn in Symyx Drawer and then subjected to generate 3D molecular structure by energy minimizing them with MMFF force field (Halgren, 1996) followed by considering distance-dependent dielectric constant of 1.0 and convergence criterion of 0.01 kcal/mol A.

Prediction of ADME descriptors and toxicity
The major reason for the failure of most of the drug candidates during clinical trial are poor ADME and high toxicity profile. Thus an important aspect of drug discovery is to evaluate critical physio-chemical as well as toxicity profile in initial stages of drug discovery (Ray et al., 2010). Hence, to select a potential drug candidate, all 13 hypothetical ligands were screened on the basis of ADME and toxicity filter using QikProp 4.3 from Schrodinger 2013 software. QiKProp used Maestro-formatted .mae files (also known as 3D SD files) as input. Output of QikProp showed a number of principal descriptors and ADME properties as shown in Table S1 (Choudhary et al., 2015). Lipinski rule of 5 was applied on the ligands and hits which were having more than 1 violation were rejected. Various other physicochemical properties were also calculated, represented by different descriptors such as molar weight (MW), number of rotatable bonds (NRB), lipophilicity parameter [log P(o/w)], number of hydrogen bond acceptors (HBA), number of hydrogen bond donors (HBD), total polar sur-face area (TPSA), solubility (log s), solvent-accessible surface area (SASA), skin permeability (log Kp), binding to human serum albumin (log Khsa), blood-brain partition coefficient (logBB), apparent MDCK cell permeability (affyPMDCK), apparent Caco-2 cell permeability (affyPCaco), percentage human oral absorption (Dash et al., 2015). Ideal ranges of various descriptors calculated with the reference to 95% of drugs are presented in Table S1.

ADME/T prediction in Discovery Studio
For additional validation purpose, Discovery studio 2.5 (Accelrys, San Diego, CA, USA) was used to describe absorption, distribution, metabolism, elimination, and toxicity (ADME/T) properties by using ADME/T descriptors module. Based on the existing information of drug, six mathematical models are used in the modules to predict quantitatively the properties by a using set of rules/keys ( Table S2) that specify threshold ADME/T characteristics for the chemical structure of the molecules. TOPKAT predictor was also selected, where the models that were used to calculate using TOPKAT are Rat Oral LD 50 (v3.1), Daphnia EC 50 (v3.1). Binding is <90% 1 Binding is ≥90% 2 Binding is ≥95%

Quantum Descriptors calculation and QSAR Analysis
All compounds were subjected to geometry optimization by using semi-empirical SCF-MO method, implemented in MOPAC 7.0. Here, basis set PM3 was used to calculate the quantam descriptors. The following quantum chemical descriptors were considered: the energy of the highest occupied molecular orbital (E HOMO ), the lowest unoccupied molecular orbital (E LUMO ), the band gap energy (∆E), the dipole moment (D), Electrophilicity (χ), Global Hardness (η). Using Koopmans' approximation, and ionization energy (I) and electron affinity (A) can be expressed in terms of the energies of the highest occupied (E HOMO ) and the lowest unoccupied molecular orbital (E LUMO ) as: Ionization potential, I= -E HOMO Electron affinity, A= -E LUMO Electrophilicity index is defined (χ) as a measure of the decrease in energy due to the maximal transfer of electrons from a donor to an acceptor system and is given as: Where, µ and η are the chemical potential and hardness, respectively. Chemical potential and hardness can be expressed in terms of ionization energy (I) and electron affinity (A) as given below: η= [HOMOε-LUMO ε]/2 After that, 2D molecular descriptors, SlogPs (according to various subdivided accessible van der Waals surface area) were calculated by the program MOE (MOE software: Chemical computing group's molecular operating environment (MOE) software, version 2014.0901). After that, PLS QSAR (Geladi andKowalski, 1986, Helland, 1988) analysis was performed to determine the relationship between these molecular descriptors and biological activity of the compounds. The maximum condition number of the principal component transform of the correlation matrix S, the condition limit, was set at 1.0 × 10 6 which is a very high setting. The leave-one-out cross validation (LOO-CV) scheme was used to validate the models and the correlation coefficient (Q 2 ) and root-mean-square error (RMSE) were reported.

ADME and toxicity screening
All 13 hypothetical compounds were screened for "druglikeliness and "non-toxic" profile on the basis of various parameters of the software QikProp. ADME filter lead to 13 drugs like hits. The ADME prediction studies data of 13 hits are shown in Table 2. The range of various important parameters were predicted, like molecular weight lied in range between 374.826 and 502.928, the value of total solvent accessible surface area (SASA) ranged 645.681-802.389, estimated number of hydrogen bonds donated (Donor HB) by the solute to water molecules in an aqueous solutions was found between 1 and 3, estimated number of hydrogen bonds accepted (accept HB) by the solute to water molecules in an aqueous solutions ranged 4-7.25, predicted octanol/gas par-tition coefficient (QPlogPoct) ranged 18.439-27.16, predicted water/gas partition coefficient (QPlogPw) ranged 10.268-16.232, predicted octanol/water partition coefficient (QPlogPo/w) ranged 3.602-6.513, predicted aqueous solubility, log S (QPlogS) ranged -5.684 to -8.177, predicted brain/blood partition coefficient (QPlogBB) ranged -0.144 to -1.429, and Lipinski violations were ≤2. Results of ADME-T studies reveal that compounds from 1 to 13 can be considered best, having druglikliness as well as non-toxicity profile.

ADME/T Prediction by TOPAK
The use of In silico approaches to predict ADME/T properties is considered as a first step in the direction to analyze the new chemical entities to reduce wasting time on lead compounds which would be toxic or metabolized by the body into an inactive form and unable to pass through membranes, and some results for this analysis are herein depicted in Table 3 along with ADME/T descriptors and their rules/keys are tabulated in Table  S2.The ADME/T profile of all the molecules under investigation was forecasted by means of six pre calculated ADME/T models provided by the Discovery Studio 2.5 program. Lipophilicity could be estimated as the log of the partition coefficient between noctanol and water (logP). Though logP is generally used to determine a compound's lipophilicity, the fact that logP is a ratio creates a concern about the use of logP to score hydrophilicity and hydrophobicity (Hughes et al., 2008). Thus the information of Hbonding prominences as obtained by calculating PSA could be taken into consideration along with logP calculation (Egan et al., 2000). Therefore, a model with descriptors AlogP98 and PSA 2D were taken into consideration for the accurate prediction for the cell permeability of compounds. According to the model for a compound to have an optimum cell permeability should follow the evaluation criteria (PSA <140 A˚2 and AlogP98 <5) (Egan et al., 2000). where each descriptor represent one physicochemical properties for "druglikeliness" such as MWmolecular weight; SASA-Solvent Assessable Surface Area; HBD-no. of H bond donor; HBA-no. of H bond acceptor; QPlogPo/w-predicted octanol/water partition coefficient; QPlogS-predicted aqueous solubility; QPlogKp-predicted skin permeability; QPlogKhsa-predicted human serum albumin binding; Affy MDCK-Apparent MDCK cell permeability; AffyCaco-Apparent Caco-2 cell permeability; QPlogBB-predicted brain/blood partition coefficient; QPlogPoctpredicted octanol/gas partition coefficient; QPlogPw-predicted water/gas partition coefficient.%HOA-Percentage Human Oral Absorption; PSA-Total Polar Surface Area; Rule of 5-Lipinski violations. All the compounds showed polar surface area (PSA) <140 A˚2. Considering the AlogP98 criteria, all inhibitors had AlogP98 value <5. Table 2 shows that among 13 compounds, 6 compounds have undefined values for BBB penetration levels (4 as mentioned in Table 2) with the exception of 7 compound having very high and high value (level 0 and 1) BBB penetration level. The aqueous solubility plays an essential role in the bioavailability of the candidate, all of them having low aqueous solubility level shown as Table S2.
All compounds have been predicted to have hepatotoxicity level of 1 which is toxic to liver and lower rate of first pass metabolism. All compounds exhibit as an inhibitors with respect to CYP2D6 liver (with reference to Table S1) except compound 1 shown in Table 3.
This indicates that they are less metabolized in Phase-I metabolism except one. Finally, the ADME/T plasma protein binding property prediction denotes that for compounds 1, 2, 3, 7, 12 binding is ≥90% and rest of them binding is ≥95% respectively, (refer to Table S1), clearly suggesting that all of them have good bioavailability and are not likely to be highly bound to carrier proteins in the blood. The Rat Oral LD 50 values of all the compounds are within the Optimum Prediction Space (OPS). And the LD 50 values are satisfactory. These high LD 50 values indicate higher safety of these compounds ( Table 4). The values obtained from daphnia EC 50 model show remarkable EC 50 values ( Table 5).

QSAR models for inhibitor activity
The quantum chemical descriptors calculation that obtained from geometry optimization like, dipole moment, HOMO energy, LUMO energy, ionization energy, electron affinity, chemical potential, hardness, and electrophilicity of thirteen molecules are represented in Table S3. Parameters like dipole moment, LUMO energy, hardness and electrophilicity were considered to predict model, since these parameters are widely used to describe the stability, toxicity and reactivity of the chemical compounds (Kapur et al., 2000). We also considered the hydrophobicity, as many authors reported the quantitative correlation of toxicity, activity and metabolism of the compounds with the octanol-water distribution coefficient descriptor (Bundy et al., 2001, Ren and Frymier, 2002, Worgan et al., 2003. In order to establish a better predictability, PLS based QSAR study was performed by considering experimentally observed values (pIC 50 ), along with various descriptors, are presented in Tables 1. According to the analysis, the best model included the following predictors: SLogP, dipole moment (D), and electrophilicity (χ). The model is expressed in the Eq. 1.
pIC 50 = -6.21625 + 0.18964(D) + 0.55132(χ) + 0.05674(SlogP_VSA9) + 0.00000(SlogP_VSA8) + 0.01887(SlogP_VSA7) + 0.00000(SlogP_VSA6) + 0.02588(SlogP_VSA5) + 0.20657(SlogP_VSA4) +  The model considers the higher contributions of SLogP, dipole moment (D) and electrophilicity (χ) with inhibitory activity. A plot of the predicted activity versus experimental activity for molecules using training set for model of KDR inhibition is shown in Figure 1. The contribution of global hardness (η) was also found in the correlation of biological activity, along with dipole moment (D) and SlogP, which is shown in eq. (3). The correlations that found in eq. 3, is rendered in Figure 3.  Studies relating the molecular properties of the compounds and their KDR inhibitory activity supported this hypothesis, as the descriptor dipole moment and SLogP show the higher contribution to the activity against KDR in the QSAR model, where contribution of electrophilicity is more significant than other descriptors. The higher LUMO energy is correlated with higher electron affinity. Henceforth, it concludes that increasing of electron affinity can increase the KDR inhibition activity. It is remarkable to note that increase in global hardness of the molecule leads to increase in stability and decrease in reactivity of the species (Ayers and Parr, 2000). From eq. (3), it also confirms that hardness of these molecules influence their inhibitory activity. However as a whole, this result is anticipated as the hydrophobicity and lipophilicity of the chemical compounds mainly govern their biological actions at cellular and molecular levels.

CONCLUSION
The current study is designed to reveal the therapeutic potentiality of some selected N-4-chlorophenylnaphthamide derivatives regarding the perspectives in tyrosine kinase inhibitor. In collateral, our findings agree with the existing evidence with concluding a better pharmacokinetics profile. The results of QSAR and ADME/T studies are validated each other and the developed mathematical models could provide insight into the structural requirements for the synthesis of new potential chemical structure having a better KDR inhibition activity.