A novel series of Benzothiazepine derivatives were designed and evaluated for epidermal growth factor receptor (EFGR) inhibitory potential by using various computational tools along with MTT assay based cell line studies towards anti-cancer activity. Our docking studies evidenced that all the present investigated twenty compounds have the potential to dock inside the active binding pocket of the EGFR Kinase domain with a binding energy in a range of -7.94 to -9.71 Kcal/mol. On the other hand, MTT assay based studies against DU145, MCF7 and HT29 cell lines showed good correlating results to our In-silico predictions. Moreover, all the designed compounds were predicted to be having promising ADMET parameters and found to be well in compliance with Lipinski’s rule of five along with no toxicology profile expect compound 7 based on Osiris property explorer predictions. Among all the twenty compounds tested, compound 9 is the best lead like molecule with -9.71 kcal/mol of binding energy with predicted IC50 value of 76.70 nano molar along with experimental IC50 values of 16±1, 27±1 and 28±1 in µg/ml for DU145, MCF7 and HT29 cell lines respectively. Molecular dynamic simulation studies for this compound 9 in complex with EGFR kinase domain has elucidated several interesting molecular level protein-ligand interactions with some of the important amino acid residues present at the active binding site of EGFR Kinase domain. Conclusively, novel designed compound 9 of the present study have shown promising anti- cancer potential worth considering for further evaluations.
Keywords: EGFR kinase, MTT assay, docking, ADMET, MD simulations, Benzothiazepine, anti-cancer.
|Citation: CH.M.M.Prasada Rao et.al. (2018) Novel series of 1, 5 Benzothiazepine skeleton based compounds as anti-cancer agents – In silico and MTT assay based study. Journal of PeerScientist 1(2): e1000008. doi : 10.5281/zenodo.3372436|
|Received: October 12, 2018; Accepted: December 18, 2018; Published: December 30, 2018|
|Copyright: © 2018 CH.M.M.Prasada Rao et.al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.|
|Data Availability: All relevant data are within the paper and its Supporting Information files.|
|Funding: This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.|
|Competing interests: The authors have declared that no competing interests exist.|
|* E-mail: email@example.com, firstname.lastname@example.org | Phone: +91 9177247605|
Epidermal growth factor receptor (EGFR) is one among the most important receptors involved in angiogenesis, one of the crucial requirements of cancer spread [1-3]. This receptor plays an important role in the progression of many solid tumors [4-5]. Particularly, the EGFR family of tyrosine kinases involvement of the cancer proliferation suggests that a small molecule such as an inhibitor to the receptors related to EGFR family could have significant therapeutic potential . The benzothiazepines  (figure 1a,b) are 7 membered heterocyclic compounds. Characteristically, this class of compounds contains sulfur and nitrogen and posses varied bioactivities useful in drug discovery [8-15].
Benzologs of 1,4-thiazepine (figure 1c) are among the important representatives of 1,5-Benzothiazepines and one among the three plausible benzo-condensed derivatives, viz. 1,4-(figure 1d), 4,1- (figure 1e) and 1,5-benzothiazepines [16-19]. Because of the high target specific binding potency of 1,5- benzothiazepine derivatives, they are always been of special interest in lead like drug candidates discovery [20-34]. Diltiazem is the first small molecule derivative of 1, 5- benzothiazepine which is been used clinically (figure 1f), followed by clentiazem (figure 1g), against cardiovascular arrests. Whereas, thiazesim (figure 1h), clothiapine (figure 1i) and quetiapine (figure 1j). Based on the wide range of applications of these compound derivatives in curing many diseases; in the present study, we have designed some novel derivatives (table 1) and investigated their potential to act as anti-cancer lead like drug candidates targeting EGFR kinase domain.
Figure 1: Structures of representative benzothiazepines compounds: 1) and 2) 1,5-benzothiazepines 3) 1,4-thiazepine 4) 1,4-benzothiazepines 5) 4,1-benzothiazepines 6) diltiazem 7) clentiazem 8) thiazesim 9) clothiapine and 10) quetiapine.
Docking of the compounds with EGFR Kinase domain active site:
We have performed the docking studies for the present studied twenty compounds with the EGFR Kinase domain protein targeting its active binding site in order to know the binding energy involved in this complex formation and to know the molecular interactions responsible for this target specific inhibition. Docking results are tabulated in table 2. All the twenty compounds studied in this present work have shown to be successfully docking inside the active site of EGFR kinase domain with a binding energy in a range of -7.94 to -9.71 Kcal/mol. As per the molecular docking results, it was revealed that Compound 9 has the best estimated -9.71 Kcal/mol of binding energy for the EGFR kinase domain inhibited complex formation by forming a hydrogen bond with LYS723 along with a pi-pi stacking with PHE699. Apart from hydrogen bonds and pi-pi stackings, this compound was also found to be forming Electrostatic interactions with PRO770; GLY772; MET769; LEU764; ILE765; ALA719; THR830 and GLU738 along with Van der waals interactions with ASP831; LEU768; LEU694; LEU820; THR766; ILE720 and VAL702 residues (Figure 2).
Figure 2: Docking snapshot of the best compound (compound 9) in complex with EGFR kinase domain.
Half maximal (IC50) inhibitory concentration is a good indicative parameter towards evaluating the plausible experimental activity of the compounds. Table 2 shows the predicted IC50 value for the compounds 1-20 were in a range of 76.70 to 1510.00 nano molar. Among which the compound 9 has shown the best possible inhibitory potential with 76.70 nano molar, whereas compound 12 with least predicted IC50 value of 1510.00 nano molar. Promisingly high IC50 values have been obtained for the present studied compounds clearly demonstrating their inhibitory potential with kinase domain of EGFR.
MTT assay based Cytotoxicity evaluation on various cell lines:
Cytotoxicity of the present investigated compounds has been evaluated on three different cell lines HT-29; MCF-F and DU-145 as per the default protocol. Results have been presented in table 2. Of all the compounds tested against HT-29 cell lines, the compound BP9 having p-nitro moiety in its structure showed maximum activity with a IC50 value of 28 µg/mL. This is followed by compounds, BP11 having a nitro methyl moiety (IC50 36 µg/mL), BP4 and BP6 having fluorophenyl and dichlorophenyl moieties respectively (IC50 42 µg/mL), BP8 having an o-nitro moiety (IC50 55 µg/mL) and BP1 having a methyl moiety (IC50 56 µg/mL). The other compounds also showed activity but at a higher IC50 values compared to control methotrexate which has shown IC50 11 µg/mL.
Among the compounds tested for cytotoxicity on MCF-7 cell lines, the compound BP9 showed maximum activity (IC50 27 µg/mL). This was followed by compounds, BP11 (IC50 28 µg/mL), BP4 (IC50 42 µg/mL) and BP6 (IC50 48 µg/mL). All the other compounds showed cytotoxicity at higher values compared to control methotrexate which has shown IC50 62 µg/mL.
Among the compounds tested for cytotoxicity on DU-145 cell lines, the compounds, BP9 and BP11 showed maximum activity (IC50 16 µg/mL). This was followed by compounds, BP4 (IC50 33 µg/mL), BP6 (IC50 46 µg/mL), BP8 (IC50 52 µg/mL) and BP1 (IC50 56 µg/mL), whereas control methotrexate has shown IC50 56 µg/mL. It was also observed that among all the compounds tested on these three cell lines, most of the compounds showed maximum activity on prostate cancer cell lines (DU-145) and the compound shows better cytotoxic activity on these three cell lines.
Prediction of pharmacological properties:
Pharmacological properties of the present studied compounds have been evaluated using Osiris Property Explorer online server according to Lipinski’s Rule of Five  for Oral Bioavailability (supplementary table 1). For the tested compounds it was inferred that all the present studied compounds are compiling according to Lipinski’s rule of five for good bioavailability. The result of toxicity analysis of all the analyzed compounds is shown in (supplementary table 2). “None” means low toxic tendency, “medium” means the midcore and “high” means high tendency of toxicity. Supplementary table 1 shows that all the analyzed twenty compounds has low or no toxic tendency however compound 7 has predicted to be having high tumorigenic toxicity.
Docking and drug likeliness based screening of best compound:
Keeping in view of binding energies, IC50 values, ADMET parameters along with predicted toxicology profile of the present investigated compounds it was found that compound 9 has the promising anti cancer drug like properties based on its ΔG binding energy and IC50 value. Based on Pharmacological properties, all the twenty compounds showed good pharmacological attributes. These compounds were found to comply with Veber’s rule, Lipinski’s rule and oral bioavailability parameters. Whereas, compound 9 showed best pharmacological attributes with well below toxicity threshold along with highest binding affinity and least half inhibitory potential.
MD simulations of EGFR Kinase domain protein in complex with compound 9:
In order to analyze the stability and conformational changes along with the aim of revealing their underlying molecular interactions at atomic level of EGFR Kinase domain in complex with our best predicted drug like compound (compound 9) among the twenty present investigated compounds. Before starting the analysis part of the MD simulations trajectories, firstly we have analyzed the simulation quality parameters, which is an important issue to deal with to conform that the molecular dynamic simulations were carried out under our given temperature, pressure and simulation box volume conditions throughout the simulated time. For the molecular dynamic simulations, input parameters given such that the temperature is at 300 k; pressure is at 0 atmospheric and volume of the simulations water box is kept at 365000 Å3 throughout the simulated time. As evident with supplementary figure 1, the given input simulation parameters were maintained thoroughly throughout the simulated time of 10 ns each, thus conforming the quality of the simulations carried out.
Simulation event analysis: After confirming the simulation quality parameters, we have studied the Root mean square deviation (RMSD) [Figure 3a], Root mean square fluctuations (RMSF) [Figure 3b], Radius of gyration (ROG) [Figure 3c], intra molecular hydrogen bonds [Figure 3d]; and total energy [Figure 3e] contributions with time dependant function of MD simulations in order to understand the stability and conformational changes of the EGFR Kinase domain protein in complex with compound 9 at each frame of the MD simulations trajectory. The statistical data for the analysis of the MD trajectories has been presented in table 3.
Protein’s backbone RMSD (figure 3a) was shown to be fluctuating between 0.8 and 2. 5 Å with an average of 1.66 Å; which was comparatively minimized than the RMSD of apo protein backbone. Whereas, cα, sidechains and heavy atoms in the protein were found to be fluctuation upto 3.3 Å. When, the ligand RMSD corresponds to protein’s backbone along with its own structure, it was found to be fluctuation around 6.4 Å with correspondence to the protein structure but found to be quite stable to its own structure by showing well below 0.8 Å deviations throughout the simulation time.
Each residue fluctuations (RMSF) (figure 3b) were calculated and averaged for the entire 10ns simulated timescale in order to identify the higher flexible regions in the protein. When the RMSF graph of EGFR kinase domain in complex with compound 9 was analyzed high peaks of fluctuations were observed at the residues positioned at the first 40, 175 and 245 positions apart from the start and ending of the protein structure. On the other hand, above mentioned residues might be plausibly involved towards forming the minimized conformational changes in the protein in presence of this compound 9.
Protein structure compactness can be predicted using radius of gyration (ROG) . We have analyzed the EGFR Kinase domain in complex with compound 9 was measured in order to estimate the effect of this compound on the overall compactness and stability of the protein and found out that the ROG of the protein in complex with this compound is in a range of 19.7 to 20.3 with an average of 20.0 Å statistically (figure 3c). Total number of intra molecular hydrogen bonds is a valuable parameter to the underlying forces contributing for the stability of the protein structure; higher intra molecular hydrogen bonds mean increased rigidity in the protein.
From the graphs in (figure 3d) it is evident that the EGFR protein structure in complex with compound 9 is strengthened by 216 intra molecular hydrogen bonds in a range of 191 to 243 table 3. This data clearly evidences that our proposed compound 9 has the potential to inhibit the EGFR kinase domain protein by causing contraction in the proteins overall conformation resulting in increased intra molecular hydrogen bonds contributing for the proteins rigidity. Increased rigidity in the protein is a direct indication of reduced activity of its functionality.
Finally, the total energy of EGFR protein in complex with compound 9 was determined (figure 3e). It was observed that EGFR protein in complex with compound 9 has an averaged -9016 Kcal/mol of energy in a range of -10023 to -8188 Kcal/mol. Above mentioned energies are calculated for the protein alone in complex with compound 9. An increase in the overall energy can be seen after the first 2.5ns of simulated time, however found to be lowering towards the end of the simulation. We have also monitored the total secondary structure elements (SSE) present in the protein like alpha helices and beta strands throughout the simulation trajectory. From the analysis it was revealed that the protein EGFR in complex with compound 9 was maintaining an average of 47% of SSE composition of helices and strands over simulated time. Most of the protein was stabilized with strands (blue) and helices (red) followed by loops (white) (supplementary figure 2). Ligand properties profile in complex with EGFR Kinase domain during MD simulations: Simulation interaction diagram program integrated within Desmond software has a novel feature which reports ligand properties such as RMSD, radius of gyration of the ligand; intra molecular H-Bonds; molecular surface area; polar surface area and solvent accessible surface area of the ligand during the simulated length of time scale. As evident with supplementary figure 3, ligand RMSD was found to be fluctuating between 0.5 to 1.2 with an average of 0.6 Å which is good enough to state that the compound is quite stable inside the active site of the protein throughout the simulated time. As the simulation progresses the ligand was found to be contracting at initial stage in its size as evident with its radius of gyration map declining from 36 to 25 Å, however stabilized at an average of 35 Å. Interestingly no intra molecular H-Bonds were found to be responsible for this contraction of the ligand. Molecular surface (MolSA) area, solvent accessible surface area (SASA) along with polar surface area was found to be maintaining an average of 336, 100 and 112 Å2 respectively.
Figure 3: a) Root mean square Deviation (RMSD) b) Root mean square Fluctuations (RMSF) c) Radius of Gyration (ROG) d) Total intra molecular hydrogen bonds and e) Total Energy graphs of EGFR Kinase domain protein in complex with compound 9.
Molecular interactions of EGFR Kinase domain-compound 2 complex during MD simulations:
Detailed molecular interactions between EGFR Kinase domain and compound 9 were studied using Desmond’s simulation interaction diagram program. There were about 21 contacts found in between EGFR Kinase domain and compound 9 in total among which one hydrogen bond was observed with MET769 residue and most of other contacts were found to be hydrophobic contacts followed by water bridging and ionic bonds (Figure 4).
Figure 4: Protein-ligand interactions (or ‘contacts’) of EGFR Kinase domain- Compound 9 complex. The stacked bar charts are normalized over the course of the trajectory: for example, a value of 0.7 suggests that 70% of the simulation time the specific interaction is maintained. A schematic detailed compound 9 atom interactions with the EGFR Kinase domain protein residues as observed during simulation.
Table 4 shows the molecular interactions profile of compound 2 with EGFR kinase domain. Hydrophobic contact were found to forming with residues LEU692; LEU694; SER696; PHE699; VAL702; ALA719; LYS721; MET742; LEU768; PHE771; CYS773; ASP776; ARG817 and LEU820. Several hydrogen bonds were observed facilitated by water brides with the residues LYS704; LYS721; MET742; THR766; GLN767; PRO770; GLY772 and ASP831. Apart from above molecular interactions PHE699 residue was found to be playing a crucial role in binding compound 9 via pi-pi stacking towards stabilizing the protein-ligand complex in its active site (figure 4). Finally, we have analyzed the rotatable bonds present in the ligand. For the compound 9, a total of three rotatable bonds have been observed all of which are were found at the linker chains. From Figure 5 it is clear that the rotatable bonds in compound 9 are consuming energy of 5.70; 0.58 and 5.60 Kcal/mol of energy.
Figure 5: The ligand torsions of compound 9 plot summarizing the conformational evolution of every rotatable bond (RB) in the ligand throughout the simulation trajectory (0.00 through 10.00 nsec). The top panel shows the 2d schematic of a ligand with color-coded rotatable bonds. Each rotatable bond torsion is accompanied by a dial plot and bar plots of the same color. Dial (or radial) plots describe the conformation of the torsion throughout the course of the simulation. The beginning of the simulation is in the center of the radial plot and the time evolution is plotted radially outwards. The bar plots summarize the data on the dial plots, by showing the probability density of the torsion. If torsional potential information is available, the plot also shows the potential of the rotatable bond (by summing the potential of the related torsions). The values of the potential are on the left Y-axis of the chart in Kcal/mol.
Benzothiazepine analogs designing and their cytotoxicity potential was evaluated targeting EFGR receptor in computational approaches along with cell line based studies. All these twenty derivatives showed potent EFGR receptor potential. Our In silico studies provided a rationalization to the ability of present studied novel compounds as a valuable small inhibitor compound with high binding affinity towards EGFR Kinase domain consolidating their complex’s thermodynamic stability for plausible anti-cancer activity. Moreover our hypothesis further substantiated by the predicted IC50 values of these compounds that they have the plausibility to inhibit EGFR Kinase domain. Further, de novo simulations for 10 nano seconds revealed ligand interaction with the residues of EGFR Kinase domain, all or some of which fall under catalytic active site important residues for its structurally stability and/or functionality. Knowledge gained through this study would be of high value towards enhance the discovery of EGFR Kinase domain target specific drug molecules. The present study further evidence that EGFR kinase domain might be the drug target for the present investigated compounds for their anti-cancer property. Moreover, promising binding energies with correlating best activity predictions of the cell line based studies for the present compounds especially compound 9 substantiates the need of further evaluating this compounds ability to inhibit cancer.
Software and program:
Chemsktech and Accelrys Discovery Studio (Version 4.0, Accelrys Software Inc.)  were utilized to draw the ligand compounds and visualize the protein-ligand interactions and to render images respectively. Autodock 4.0  is used to perform the semi-flexible protein ligand docking studies as per the protocol followed elsewhere [39-44]. ADMET properties of the compounds were determined using Molinspiration and Orissis property explorer online programs. Schrodinger’s Desmond module v3.6 [45-46] was used for molecular dynamic simulation studies between EGFR kinase domain crystal structure [PDB: 1M17] and compound as per the standard protocol explained in detail elsewhere [47-48]. Periodic boundary conditions were used to determine the specific size and shape of the water box buffered at 10 Å distances and box volume was calculated as 365000 Å3. Seven Cl- ions for EGFR-compound 9 compound complex were added to neutralize the system electrically and placed randomly in the solvated system.
MTT assay method for cytotoxic studies:
The compounds BP1-BP20 compounds were evaluated for their cytotoxic studies on different cell lines using the MTT assay method [49-52]. The results were compared with standard drug Methothrexate. The cell lines used in this present study are HT-29 (colon cancer), MCF-7 (breast cancer) and DU-145 (prostate cancer), obtained from National Centre for Cell Science (NCCS), Pune, India. Cell lines were grown as adherent in DMEM media and the culture was maintained in a humidified atmosphere with 5% CO2. DMEM (Dulbeccos Modified Eagels Medium), MEM (Minimum Essential Media Eagle), MTT [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide], Trypsin, EDTA were purchased from Sigma chemicals (St.Louis,MO). Fetal bovine serum (FBS) was purchased from Arrow Labs, 96 well flat bottom tissue culture plates were purchased from Tarson. Stock solutions of test compounds (BP1 to BP20) were prepared (10 mg/mL) in DMSO and from that various dilutions were made with sterile water to get the final drug concentrations. The cells are then solubilized with DMSO and the released, solubilized formazan reagent is measured spectrophotometrically at 570 nm. Since reduction of MTT can only occur in metabolically active cells, the level of activity is a measure of the viability of the cells.
Supplementary table 1: Structures of novel designed 3, 5-dinitrophenyl clubbed azoles used in the present study.
Supplementary table 2: Toxicity of sample compounds based on Osiris Property Explorer predictions.