Docking and ADME / T analysis of silibinin as a potential inhibitor of EGFR kinase for ovarian cancer therapy

1 Molecular Modeling & Drug Design Laboratory (MMDDL), Pharmacology Research Division, Bangladesh Council of Scientific & Industrial Research (BCSIR), Chittagong-4220, Bangladesh. 2 Department of Pharmacy, International Islamic University Chittagong, Chittagong-4203, Bangladesh. 3 Department of Pharmacy, State University of Bangladesh (SUB), Dhaka1205, Bangladesh. 4 Department of Pharmacy, Jahangirnagar University, Savar, Dhaka-1342, Bangladesh. 5 Superintendent of Drugs, Directorate General of Drug Administration, Bangladesh. 6 Department of Pharmacy, BGC Trust University Bangladesh, Chittagong-4000, Bangladesh.


INTRODUCTION
Recently, the leading cause of death in women is ovarian cancer amounting to 1,40,000 deaths every year (Jemal et al., 2011).Though most women with epithelial ovarian cancer (EOC) are diagnosed at advanced stages, the poor prediction of this disease causes more death (Slatnik and Duff, 2015).As a result, the clinical response rate at the outset is very high, where platinum based compounds are benchmark first-line agents for ovarian cancer treatment (Herzog, 2004).However, due to the development of acquired chemoresistance in most ovarian cancer, the treatment strategies are now in a big challenge (Agarwal and Kaye, 2003;Rabik and Dolan, 2007).Therefore, in such condition there is an urgent necessity to develop novel and effective therapeutic strategies against ovarian cancer.(Hennessy et al., 2009;Bible et al., 2012).Epidermal Growth Factor Receptor (EGFR), which belongs to the family of Human Epidermal Receptor (HER), is highly expressed in multiple malignancies, mainly in ovarian cancer (Bull Phelps et al., 2008).It's one of the tyrosine kinases and mostly activated by extracellular ligands that trigger receptor auto phosphorylation, which may lead to the activation of downstream pathways involved in proliferation, survival, angiogenesis, and invasion (Scaltriti and Baselga, 2006;Lemmon and Schlessinger, 2010).Recent studies have indicated that the EGFR pathway is recurrently activated in most cancer cell, and its targeted inhibition with a small molecule kinase inhibitor has been successful for lung, breast and ovarian cancer with EGFR mutations (Yarden and Sliwkowski, 2001;Cataldo et al., 2011).In this context, silibinin is a member of flavanoids that has high pharmacological and therapeutical use.
It is one of the most active of three isomers that can be isolated from the seeds of the milk thistle, silyburn Marianum (Ferenci et al., 2008).Several studies have indicated that silibinin possesses chemopreventive and therapeutic properties (Rahman and Sarkar, 2005;Kandala and Srivastava, 2010).Therefore, the present study was aimed to investigate the therapeutic properties of silibinin as EGFR tyrosine kinase inhibitor through various computational analysis, including molecular docking (Dash et al., 2014;Dash et al., 2015b;Raju Dash et al., 2015) and ADME/T analysis.

Protein Preparation
Three dimensional crystal structure of epidermal growth factor receptor kinase domain (PDB id: 1M17) was downloaded in pdb format from the protein data bank (Berman et al., 2000).After that, structure was prepared and refined using the Protein Preparation Wizard of Schrödinger-Maestro v9.4.Charges and bond orders were assigned, hydrogens were added to the heavy atoms, selenomethionines were converted to methionines, and all waters were deleted.Reorientation of certain hydroxyl and thiol groups, amide groups of asparagines, glutamines and imidazole ring of histidines, protonation states of histidines, aspartic acids, and glutamic acids was optimized at neutral pH.Using force field OPLS_2005, minimization was carried out setting maximum heavy atom RMSD to 0.30 Å.

Ligand Preparation
Compound was retrieved from Pubchem databases, i.e. silibinin (CID 31553).The 3D structures for these were built by using Ligprep2.5 in Schrödinger Suite 2013 with an OPLS_2005 force field.Their ionization states were generated at pH7.0±2.0 using Epik2.2 in Schrödinger Suite.Up to 32 possible stereoisomers per ligand were retained.

Receptor grid generation
Receptor grids were calculated for prepared proteins such that various ligand poses bind within the predicted active site during docking.In Glide, grids were generated keeping the default parameters of van der Waals scaling factor 1.00 and charge cutoff 0.25 subjected to OPLS 2001 force field.A cubic box of specific dimensions centred on the centroid of the active site residues (Reference ligand active site) was generated for receptor.The bounding box was set to 14 Å × 14 Å × 14 Å for docking experiments.

Glide Standard Precision (SP) ligand docking
SP flexible ligand docking was carried out in Glide of Schrödinger-Maestro v9.4 (Friesner et al., 2004;Friesner et al., 2006) within which penalties were applied to non-cis/trans amide bonds.Van der Waals scaling factor and partial charge cutoff was selected to be 0.80 and 0.15, respectively for ligand atoms.Final scoring was performed on energy-minimized poses and displayed as Glide score.The best docked pose with lowest Glide score value was recorded for each ligand.

Prime MM-GBSA
Prime MM-GBSA approach was used to calculate ligand binding energies and ligand strain energies for a ligand and a single receptor.MM-GBSA is a method that combines OPLSAA molecular mechanics energies (EMM), an SGB solvation model for polar solvation (GSGB), and a non-polar solvation term (GNP) composed of the non-polar solvent accessible surface area and van der Waals interactions.Here, the Glide pose viewer file of the best conformation chosen was given as the source in Prime MM-GBSA simulation.The total free energy of binding: ΔG bind = G complex -(G protein + G ligand ), where G = EMM + GSGB + GNP

Ligand based ADME/Toxicity prediction
The QikProp module of Schrodinger is a quick, accurate, easy-to-use absorption, distribution, metabolism, and excretion (ADME) prediction program design to produce certain descriptors related to ADME.It predicts both physicochemical significant descriptors and pharmacokinetically relevant properties.ADME properties determine drug-like activity of ligand molecules based on Lipinski's rule of five.ADME/T properties of the compound (Silibinin) was analysed using Qikprop 3.2 module (Natarajan et al., 2015).

Structure-based P450 site of metabolism
The P450 Site of Metabolism (SOM) of selected compound was determined by 'induced fit docking' on CYP2D6 isoform using the 'P450 Site of Metabolism' module of Schrodinger Suite v9.3.For a given atom in a molecule to be a significant site of metabolism by a CYP450 isoform, it must have some degree of 'reactivity' in the absence of the enzyme and also be accessible to the reactive heme iron center.To address both these requirements, the 'P450 Site of Metabolism' workflow which combines induced-fit docking (IFD) for the determination of accessibility to the reactive center, with a rule-based approach to intrinsic reactivity, was used.

RESULT Molecular docking and binding energy analysis
In order to study the interaction of silibinin with EGFR kinase, we performed Glide docking analysis by schrodinger suite, where silibinin processed higher docking score as shown in Table 1.As seen in Figure 1, four hydrogen bonds were formed between silibinin and EGFR kinase, where residues like ASP 813 , ARG 817 , MET 769 and ASN 818 are involved in the interactions.It was also seen that ALA 719 , VAL 702 , LEU 820 were involved in pi-alkyl interactions in the active site of EGFR kinase, when it was visualized in Discovery studio 4.1 for better understanding.Again, in order to observe the overall clarity of docking studies, we introduced Prime MM-GBSA (Adasme-Carreno et al., 2014) approach to calculate ligand binding energies.In our study, the binding energy of silibinin indicates a strong bonding interaction and the accuracy of ligand protein binding (Table 1).

ADME and Toxicity analysis Ligand based ADME/Toxicity prediction
The drug-like activity of the ligand molecule was categorized using ADME properties by QikProp module of Schrodinger.The ADME properties of silibinin were evaluated with QikProp module of Schrodinger, shown in Table 2.The selected properties are known to influence metabolism, cell permeation, and bioavailability.
All the predicted properties of the isolated compound was in the range for 95% of known oral drugs and also satisfy the Lipinski's rule of five to be considered as drug like potential.Total Solvent Accessible Surface Area in using a probe with a 1.4 radius (acceptable range: 300-1000).f Predicted IC50 value for blockage of HERG K+ channels (concern: below −5).g Predicted aqueous solubility, S in mol/dm −3 (acceptable range: −6.5-0.5).h Predicted octanol/water partition coefficient (acceptable range: −2-6.5).i Predicted human oral absorption on 0 to 100% scale (<25% is poor and >80% is high).

Structure-based P450 site of metabolism
In the era of development and use of pharmaceuticals, cosmetics, nutritional supplements and agrochemicals, the prediction of sites and products of metabolism for xenobiotic and endogenous compounds is an important avenue of research.Here, cytochrome P450 enzymes play an integral role in the metabolism pathway of drugs by undergoing a variety of reactions, such as hydroxylation, dealkylation, and double-bond oxidation that result in the degradation of small molecules (Liu et al., 2013).
Regarding this, the metabolism of silibinin were studied by induced fit docking with CYP2D6 and the results were analyzed on the basis of "Fe accessibility", "intrinsic reactivity" and "overall SoM (Site of Metabolism)".Details of these parameters, the metabolic maps of are illustrated in Figure 2.
In Fe accessibility, labelled atoms are considered as the accessibility of the atoms to the iron.Lower of null values of silibinin indicate having no accessibility.The reactivity is predicted with a linear free energy approach based on the Hammett and Taft scheme, where the reactivity of a given atom is the sum of a baseline reactivity rate and a series of perturbations determined by the connectivity (Rydberg and Olsen, 2012).The intrinsic reactivity of silibinin in some position is quite high, where positive values are more reactive and negative values are less reactive.

DISCUSSION
From the last decades, the implication of EGFR pathway in cancer progression has been suggested by various studies (Blume-Jensen and Hunter, 2001;Chan et al., 2005).In case of ovarian cancer treatment, the outcome is often serious due to the development of platinum drug resistance.Briefly, the constitutive activation of STAT3 (signal transducers and activators of transcription 3) is linked with the over expression of EGFR that confers resistance to chemotherapy-induced apoptosis in epithelial malignancies (Coley et al., 2006;Duan et al., 2006;Dash et al., 2015a).Indeed, this protein (STAT3) is phosphorylated and abnormally activated with a high frequency in ovarian solid tumors by the action of EGFR and other tyrosine kinases (Huang et al., 2000;Burke et al., 2001).The main fact is that the activation of STAT3 at tyr-705 leads to the accumulation of antiapoptotoc events, such as Bcl-2, Bcl-xL and PI 3-kinase/akt activity (Bowman et al., 2000;Turkson, 2004;Richardson and Kaye, 2005), Literatures also suggested that, the platinum sensitivity in human ovarian cancer may be improved by the combination therapy like platinum drug with EGFR blocker.Here in this study, we have selected silibinin to investigate the binding interactions mode in the active site of EGFR kinase domains.It can be concluded from the previously outlined result section, silibinin processed a greater attraction towards the EGFR, and inhibition may be resulted by forming a stable protein ligand complex formation.Moreover, from the previously published report, silibinin was found to act as inhibitor of EGFR in breast cancer (Kim et al., 2011a;Kim et al., 2011b), in rat giloma cell line (Qi et al., 2003), in renal cell carcinoma (Li et al., 2008).Furthermore, it was also found that silibinin reduces toxicity and enhances therapeutic efficacy of doxorubicin in lung cancer (Singh et al., 2004).
The ADME and toxicity profiles of silibinin were moderated, where the oral absorption rate is medium.The concerning value of hERG K+ channel (which indicates the potential of a compound for cardiac toxicity) is also high.According to the SOM prediction, sibilinin may undergo aliphatic and aromatic hydroxylation, o-delakylation, due to having higher reactivity of silibinin (Zheng et al., 2009) (Shown in Figure 2).

CONCLUSION
From the perspectives of platinum resistant ovarian cancer, the current study is an attempt to explore out the therapeutic potentiality of silibinin compound.In parallel to the existed literature data, our experimental data not only supports the evidence but also explored out the binding mechanisms.Furthermore, the study is also useful for further clinical based studies and also for the validation of toxicological and pharmacokinetic study.

Fig. 1 :
Fig. 1: Molecular Docking analysis of silibinin and EGFR kinase protein complex obtained from Glide docking.silibinin molecule is shown as stick, protein as cartoon, and H-bonds as pink lines.

Fig. 2 :
Fig. 2: Metabolic maps.(i) Overall SOM, (ii) Fe accessibility and its (iii) Intrinsic reactivity for silibinin .Overall SOM score is the linear combination of the accessibility and the intrinsic reactivity larger scores mean higher reactivity.Fe accessibility (← ): values indicate accessibility of the atoms to the iron center of CYP, full ray represents a unit of accessibility, with the length of the final ray representing the decimal portion.Intrinsic reactivity ( ): Labels the intrinsic reactivity of atoms, positive values are more reactive.The size of the circle indicates the overall SOM score, color of the circle indicates the atomic intrinsic reactivity and a blue ring around the circle indicates that the atom passed through the CYP filtering stage

Table 1 :
Docking results of selected hits in kcal/mol.

Table 2 :
ADME and drug likeness properties of the selected ligands by QikProp.
b HB donor c HB acceptor d SASA e QP log HERG f QP log S g QP log Po/w h % Human Oral Absorption i b Molecular weight (acceptable range: <500).cHydrogen bond donor (acceptable range: ≤5).dHydrogen bond acceptor (acceptable range: ≤10).e