Int J Pharm Pharm Sci, Vol 8, Issue 7, 262-268Original Article


IDENTIFICATION OF NOVEL INHIBITORS FOR MITOGEN-ACTIVATED PROTEIN KINASE KINASE 4 BY VIRTUAL SCREENING AND MOLECULAR DYNAMICS SIMULATION TECHNIQUES

UZMA KHANAMa, BALWANT KISHAN MALIKb, PUNITI MATHURa, BHAWNA RATHIa

aCentre for Computational Biology and Bioinformatics, Amity Institute of Biotechnology, Amity University Noida, bDepartment of Biotechnology, Sharda University, Greater Noida
Email: brathi@amity.edu        

 Received: 21 Mar 2016 Revised and Accepted: 17 May 2016

ABSTRACT

Objective: The objective of this study was to discover a therapeutic natural lead compound against mitogen-activated protein kinase kinase 4 (MKK4) employing in silico studies.

Methods: In the present study, natural compounds database was first screened for potent inhibitory activity by employing a high throughput virtual screening and molecular docking. The molecular dynamic simulation was used to analyze the stability of ligands.

Results: Top ten hit compounds obtained from virtual screening and molecular docking were analyzed for their binding poses. Molecular docking studies reveal that all ten compounds bind into the same binding pocket. Molecular dynamic simulation of ZINC06091752-MKK4 and ZINC00391412-MKK4 complexes revealed stable and potential binding mode of MKK4 to ZINC06091752 and ZINC00391412.

Conclusion: The potential binding mode of MKK4 to ZINC06091752 and ZINC00391412 was explored through molecular dynamic simulations. ZINC06091752 and ZINC00391412 have been identified as potential inhibitors of MKK4. Analysis of ligand efficiency profiles through assays would add more value to the current findings and may be beneficial in prostate cancer therapy.

Keywords: Molecular Dynamics Simulation, Virtual Screening, Molecular Docking, Prostate Cancer, ERK Kinase-1, MKK4

© 2016 The Authors. Published by Innovare Academic Sciences Pvt Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/)

INTRODUCTION

Prostate cancer is the most commonly occurring form of cancer in western men, and the second-leading cause of death due to cancer [1].The death rate is increasing due to the prostate cancer metastasis. It has been reported in the literature that mitogen-activated protein kinase kinase 4 (MKK4) is a prometastatic and involved in tumor progression. MKK4 is located on 17p11.2 chromosomal segment. MKK4 is a dual-specificity protein kinase and the member of mitogen-activated protein kinase(MAPK) family.MKK4 plays an important role in MAPK signaling pathway.MKK4 directly phosphorylates serine/threonine, and also tyrosine residues [2] and, activates two downstream pathways, c-Jun N-terminal kinase (JNK) as well as p38 MAPK [3] This activated p38 MAPK triggers the metastatic cascade steps [4], followed by the epithelial to mesenchymal transition, cellular invasion, and metastatic colonization [5].MKK4 has been found a therapeutic target for the prostate cancer metastasis. Prostate cancer therapies which are used for prostate cancer treatment have lots of side effects on human health.Natural compounds have been found very effective from ancient time. Therefore, there is a need to discover natural compounds for prostate cancer treatment.As drug discovery process is very expensive and time-consuming process, in silico studies provide a cost effective platform for discovery of effective molecules. In the present work in silico studies were performed to discover natural compounds as MKK4 inhibitor.

Three crystal structures of human MKK4 have been reported till date having different resolutions [6]. Crystal Structure of human non-phosphorylated MKK4 kinase domain complexed with phospho-amino phosphonic acid-adenylate ester having high resolution among all three crystal structures was selected for this study. Virtual screening and molecular docking were carried out to identify the new MKK4 inhibitors. The molecular dynamics (MD) simulation was carried out to study the interaction pattern and stability of the ligand-MKK4 complex. MD simulation studies not only exhibits the stable interaction of the ligand with MKK4 but also accentuated the significance of given interactions in MKK4 inhibition. Prospective experimental validations can be performed from the predictions obtained from the study.

MATERIALS AND METHODS

Protein preparation and structure refinement

X-ray crystal structures for the MKK4 (3ALN) with a resolution of 2.3 Å was downloaded from protein data bank [7]. Protein structure was prepared employing protein preparation wizard of the Schrodinger suite [8] to ensure the quality and reliability of the structure. OPLS _2005 force field was applied for the protein optimization and minimization.

Library selection and ligand preparation

ZINC biogenic compounds (Zbc), a commercially available primary and secondary metabolite database, originally includes 189466 compounds. Ligands were prepared using LigPrep module from the Schrodinger package [9]. Variants of the same ligand with different stereochemical, tautomeric, and ionization properties were generated. Approximately 276784 chemical structures were produced. Energy minimization for all the produced structures were done in order to obtain conformationally relaxed structures employing OPLS_2005 force field.

Receptor grid generation

Glide [10] was used for the grid generation. The prepared 3D structure of target protein was used to generate glide scoring grid for the successive docking calculations. The receptor grid for protein was generated with the centroid of the active site residues.

High-throughput virtual screening

The prepared library was used for virtual screening workflow (VSW). In VSW, the compounds were screened for ADME property calculation using QikProp module [11]. ADME descriptors are usually calculated at the last stage of the drug discovery process; however, in the present study, we performed the ADME property calculation in the preliminary stages of drug development to eliminate the false-positive molecules with poor pharmacokinetic properties in the initial stage of the drug discovery process. QikProp was executed for physicochemical properties of all ZINC [12] compounds. Lipifilter properties were used to assess the drug likeliness of all ZINC compounds based on lipinski rule of five [13]. Virtual screening was performed in the VSW, to identify ligands that were likely to interact with binding sites of the MKK4. The prepared ligands were docked against the receptor grid. The high throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP) module of glide was assigned. Glide score (G-score) was the bases for the identification of the top lead compounds. The best conformation for each ligand was identified and, glide energy was calculated using glide score in order to obtain the output for a specific ligand single best pose.

Molecular dynamics (MD) simulation

Before MD simulation, XP docked complexes were prepared in the same manner as prepared for virtual screening. MD simulations were performed to obtain the most stable conformation of MKK4 and ligand complex, The Desmond 3.1 MD package [14] was used for MD simulations with OPLS-AA force field parameter. An orthorhombic box with periodic boundary conditions was used to solvate the protein by adding SPC water molecules.

The protein−ligand complex was energy minimized with OPLS-AA force field parameters. The energy minimized complex was subjected to steepest descent method at which maximum force being obtained in smaller than 1000 kJ/mol·nm. Prior to thermalization, the whole system underwent for energy minimization up to 1000 steps. At a constant temperature of 300 K and pressure of 1.01325 bar isothermal, isobaric (NPT) ensemble was run after the thermalization. Nose-Hover chain thermostat [15] and Martyna-Tobias-Klein barostat [16] were used to maintain the system. SHAKE [17] was applied for 2-fs time step. Particle mesh ewald method was used to treat long-range electrostatics interactions. Energies and their coordinates for the OPLS-2005 force field simulations were recorded for 2.4 ps for total 50 ns for the systems.

RESULTS AND DISCUSSION

Virtual screening and molecular docking

The increasing interest in kinase inhibitors as novel therapeutics has created a demand for the development of the selective inhibitors for mitogen-activated protein kinases. ATP-binding pockets are the primary target sites for most of the kinase inhibitors [18, 19]. Some successful approaches targeted the ATP-binding site. PD98059 and U0126 were identified specific MKK1 MKK2 inhibitors that act by binding to the inactive conformation of the enzyme to prevent its activation by Raf [20]. SB203580 and SP600125 blocked p38 and JNK kinase activity, respectively, by interacting with the ATP-binding domain [21]. HWY336 reversibly inhibit MKK4 and MKK7 [22]. In this study virtual screening was done selectively for MKK4 inhibitory activity. Crystal structure of the MKK4 shows auto-inhibition when peptide binds to its allosteric site [6]. In the present study inhibition of MKK4 take place without the interaction of substrate to the activation loop. Three-dimensional structure of MKK4 was downloaded for the use in structure-based drug design. After the cleaning of crystallo-graphic structure binding pockets were identified. Receptor grid was generated around the allosteric site. Library of biogenic compounds were docked into the receptor grid to predict the docking pose and to calculate the score. ZINC library (276784 compounds) was compressed to 164530 compounds after the filtration with QikProp and Lipifilter. In high throughput screenings, large combinatorial libraries are pre-filtered for drug likeliness and to avoid false negatives. First, a large collection of Zbc was assessed through HTVS. Subsequently, the top-ranked fraction of the screening library was gone through SP docking. In the end, lead candidates were re-docked by XP docking method. These workflows are very common in present virtual screening protocols that typically consist of a cascade of different filter approaches. Flexible docking was carried out to calculate the glide score and glide energy. 1071 compounds were screened out through HTVS; SP docking had generated 108 compounds, which were further subjected to XP docking. Finally, 91 compounds were generated in pose viewer file after the XP docking. fig. 1 represents the 2D conformation of the ligands with the top ten highest scores obtained from virtual screening. Selective top ten hits showed a significantly high affinity towards MKK4 and all ten compounds binds to the same site (Lys131) (fig. S1-fig. S8). Obtained glide energy of top ten hits are mentioned in table 1. Two compounds namely, ZINC06091752 and ZINC00391412 showed high affinity and appeared to be key lead compounds.

These compounds had shown stable interactions with MKK4. Hydrogen bond gives the most stable and specific interaction [23,24] in the processes of biological recognition. Two hydrogen-bond interactions were noted between ZINC06091752 and the active site residues (fig. 2). One hydrogen bond was formed between the oxygen atom (O20) of the ligand and hydrogen of NH2 of Lys 231 (Lig: O2···Lys 231, 2.02 Å), and the second hydrogen bond was formed between the oxygen atom (O18) of ligand and oxygen atom of Ser 233 (Lig: O2···Ser 233, 2.05 Å). ZINC06091752 was exposed to hydrophobic interaction with Leu 236, Asp 185, Ser 184, Lys 187, Arg 110 and Lys 131. The Ionic interaction was found around Lys 131 and Lys 231, Mg++was involved in the ionic interaction between the ligand and protein. Mg++showed two ionic interactions with oxygen atoms (O23 and O18) of ligand.

Two hydrogen bond interactions were noted between ZINC00391412 and the active site residues (fig. 3). The oxygen atom (O27) of the ligand formed a hydrogen bond with the oxygen atom of Ser 233 (Lig: O2···Ser 233, 2.32 Å), and the second hydrogen bond was formed between the oxygen atom (O24) of the ligand and oxygen atom of Gly 109 (Lig: OH···Gly 109, 1.79 Å). ZINC00391412 was exposed to hydrophobic interaction with Val 116, Ser 115, Arg 110, Gly 109, Gly 114, Gly 111, Lys 231, Lys 131, Ser 233, Asn 234, Ser 184 and Leu 236. The Same pattern had been observed in ionic interaction around Lys 131 and Lys 231 and, Mg++involvement in ionic binding for both the compounds. Docking studies demonstrated that ZINC06091752 exhibited the highest glide score (-10.50 kcal/mol) and glide energy (-44.87 kcal/mol). Compound ZINC00391412 demonstrated second highest glide score (-10.44 kcal/mol) and glide energy (-43.59 kcal/mol).

Table 1: Top-10 hit compounds obtained from virtual screening against MKK4 along with docking score and glide energy

 

ZINC id

G-score (kcal/mol)

Glide energy (kcal/mol)

1

ZINC06091752

-10.501

-44.879

2

ZINC00391412

-10.444

-43.591

3

ZINC02349698

-10.408

-40.729

4

ZINC00389778

-10.388

-46.149

5

ZINC14817730

-10.328

-47.315

6

ZINC00142380

-10.188

-39.435

7

ZINC00167997

-10.153

-49.897

8

ZINC00388545

-10.085

-41.996

9

ZINC13526482

-10.029

-43.788

10

ZINC01641166

-9.94

-44.845

Fig. 1: Representing top 10 hit compounds obtained from virtual screening, 2D structures and their corresponding ZINC ids


Fig. 2: (a) Binding pattern of MKK4 and, (b) 2D-interaction map of MKK4, obtained through XP docking of ZINC06091752. It illustrates existing hydrogen bonds and ionic interactions between ligand and MKK4. Lys 231 and Ser 233 shows hydrogen bonding between ligand and receptor whereas Lys 131 and Lys 231 making ionic interaction between ligand and MKK4


Fig. 3: (a) Binding pattern of MKK4 and, (b) 2D-interaction map MKK4, obtained through XP docking of ZINC00391412. It illustrates existing hydrogen bonds and ionic interactions between ligand and MKK4. Gly 109 and Ser 233 shows hydrogen bonding between ligand and receptor whereas Lys 131 and Lys 231 making ionic interaction between ligand and MKK4

Protein-ligand complex simulation

Molecular dynamics simulations of MKK4 complexes were performed in order to validate the results of virtual screening and docking studies, as well as to analyze the interaction and pattern of system stability. The structure stability was analyzed by the root mean square deviation (RMSD) values obtained from MD trajectories. Only two complexes were selected for this study based on glide score and binding energy (ZINC06091752 and ZINC00391412). In the case of bound form with ZINC06091752, it was observed that RMSD profile (2.5 Å to 4) was increasing in between 11–23 ns, while 23-50 ns significant stability was observed in RMSD (fig. 4a).

Fig. 4b shows the RMSD of the ZINCC00391412 complex following molecular dynamics simulation for 50 ns. Analysis of trajectories revealed that an increase in the RSMD profile (1.5-4.5 Å) was observed between 0-13 ns in Cα backbone, while MKK4-ZINCC00391412 system achieved lower RMSD values in the range of 0.5-2.5 Å. Analysis of RMSD plot indicates that RMSD profiles of MKK4 backbone remained stable during 13-50 ns, where was from 13-50 ns it had shown the significant stability of backbone RMSD throughout the MD runs. Around 22-23 ns a sharp disruption had been showed by the Cα backbone.

Conformational changes of the structure were analyzed throughout MD trajectories by calculating the root mean square fluctuation (RMSF) values. In ZINC06091752 complex, higher RMSF values were observed for Ala 112 (2.56 Å), Lys 140 (3.75 Å), Gln 253 (3.57 Å), Thr 318 (10.86) and Ser 84 (3.15 Å) (fig. 5a), while residues involved in inhibitor binding namely, Asp 229, Lys 231, Asn 234 and Asp 247 were quite stable. However, In ZINCC00391412 complex, higher RMSF values were observed for ALA112 (3.51 Å), Tyr 284 (3.48 Å), Thr 318 (8.90), Ser 84 (3.15 Å), Tyr 284 (4 Å). Pro 308 (3.56 Å) and Pro 389 (2.39 Å) (fig. 5b), while residues involved in inhibitor binding namely, Val 116, Ser 233, Asn 234 and Asp 247 were quite stable. RMSF values of binding site residues (Asn 234 and Asp 247) exhibited a lower tendency in both the systems. The conformational changes were observed in the nearby residues of MKK4 binding site. These changes induced significant flexibility in the protein to accommodate the lead compound in the binding site. The stability of a system can be measured through its potential energy. Potential energy plot demonstrated that both the systems remained stable throughout 50 ns MD simulations run. ZINC-06091752 system had shown lower potential energy (-225993 kcal/mol) compared to ZINCC00391412 (-225746 kcal/mol) system (fig. 6). Interaction of MKK4 with ZINC06091752 and ZINCC-00391412 were monitored throughout the MD simulation. These interactions are categorized and summarized in the form of the plot in fig. 7. The binding characteristics of MKK4 with ZINC06091752 and ZINCC00391412 were analyzed through intermolecular hydrogen bonds and ionic interaction.

In the case of MKK4-ZINCC00391412 complex, Lys 131 involved in hydrogen bonding for 30% course of the simulation. Ser 233 had shown intermolecular hydrogen bond and this interaction was maintained throughout 79% of the simulation time. Rest of the time and Lys 131 and Ser 233 involved in water bridge formation. However, in the case of ZINC06091752 system Lys 231 had shown hydrogen bond and it was maintained throughout 80% of the simulation time and, rest of the time Lys 231 involved in the ionic interaction. Asp 247 and Asn 234 shown the same ionic interactions in both the systems throughout the MD simulation. One more ionic interaction had been observed by Asp 229 in ZINC06091752 system. Compared to MKK4-ZINCC00391412 system, ionic interaction increased in the ZINC-06091752 system. Overall, interactions remained stable throughout the 50 ns simulation run.

Overall analysis substantiated that ZINCC00391412 showed more stable binding to MKK4 compared to ZINC06091752, which shows a significant distribution of RMSD values in trajectory conformation pairs. Binding in both the systems was mediated by the metal ion, Asn 234 and Asp 247 in ZINCC00391412 and, Asp 229, Asn 234, Asp 247 in ZINC06091752. These residues had not shown any structural dynamics during the MD simulation. Lys 131 shows strong contact with ligand in first 9 ns course of trajectory after that contact between the ligand and Lys 131 forms and breaks on a regular interval (fig. S9).

ZINCC00391412 inhibited MKK4 phosphorylation within the ATP-binding site. Hydrogen bonding interactions between ZINCC-00391412 and MKK4 demonstrated that this compound bound reversibly. In general, phosphorylation at specific residues (Ser/Thr and Tyr) within functions as a molecular switch for cellular communication. Our finding suggests ZINCC00391412 specifically inhibited the kinase activity of human prostate MKK4. Further experimental studies regarding the effects of ZINCC00391412 may provide important information that could be used to overcome the specificity concerns in the development of MKK4 inhibitors. These results revels the reversible inhibition of MKK4. With these findings, the identified compounds could be validated for potential inhibition through in vivo and in vitro assays.

Fig. 4: Illustrating RMSD plot obtained from 50 ns MD simulation run. (a) RMSD plot of the MKK4-ZINC06091752 system (magenta), and MKK4 without ligand (blue). (b) RMSD plot of the MKK4-ZINC00391412 system (magenta), and MKK4 without ligand (blue). X and Y axes represent time (ns), and root means square deviation (RMSD), respectively


Fig. 5: Illustrating RMSF plot MD simulation run for 50 ns. Blue color represents the unbound MKK4 and green colors shows the ligand contact. (a) ZINC06091752, (b) ZINCC00391412. X and Y axes represent residue index and RMSF, respectively


Fig. 6: Illustrating potential energy plot of MD simulation run for 50 ns. (a) Energy plot of the MKK4-ZINC06091752 system. (b) Energy plot of MKK4-ZINCC00391412 system


Fig. 7: Protein-ligand contact plot for 50 ns MD simulation. Protein-ligand contact plot shows the binding interactions, hydrogen bond (green), ionic interaction (pink), water bridge (blue) and hydrophobic interaction (purple). (a) ZINC06091752, (b) ZINCC00391412

CONCLUSION

In the current study results, top lead compounds as potential inhibitors for MKK4. All ten identified compounds bind into the same binding pocket. Lys 131 is involved in the inhibition of MKK4 by identified compounds. Out of the top, ten significant leads obtained MD simulations was used to characterize the interactions between the ligands viz. ZINC06091752 and ZINCC00391412, and MKK4 as they had the highest docking score. The study has provided an optimal result which needs to be validated in vitro and in vivo. The in silico studies have reduces the time and the cost of experimentation by providing potential leads.

ABBREVIATION

MKK4: Mitogen-Activated Protein Kinase Kinase 4, MAPK: Mitogen-Activated Protein Kinase, JNK: C-Jun N-Terminal Kinase, MD: Molecular Dynamics, Zbc: ZINC Biogenic Compounds, VSW: Virtual Screening Workflow, ADME: Absorption Distribution Metabolism Excretion, HTVS:High Throughput Virtual Screening, SP: Standard Precision, XP: Extra Precision, RMSD: Root Mean Square Deviation, RMSF: Root Mean Square Fluctuation, G-Score: Glide Score.

CONFLICT OF INTERESTS

Declared none

REFERENCES

  1. Siegel R, Miller K, Jemal A. Cancer statistics. CA Cancer J Clin 2016;66:7-30.
  2. Cuenda A. Mitogen-activated protein kinase kinase 4 (MKK4). Int J Biochem Cell Biol 2000;32:581–7.
  3. Whitmarsh AJ, Davis RJ. Role of mitogen-activated protein kinase kinase 4 in cancer. Oncogene 2007;26:3172–84.
  4. Pavese JM, Ogden IM, Voll EA, Huang X, Xu L, Jovanovic B, et al. Mitogen-activated protein kinase kinase 4 (MAP2K4) promotes human prostate cancer metastasis. PloS One 2014;9:e102289. Doi:10.1371/journal.pone.0102289. [Article in Press]
  5. Szmulewitz RZ, Clark R, Lotan T, Otto K, Taylor J, Veneris K, et al. MKK4 suppresses metastatic colonization by multiple highly metastatic prostate cancer cell lines through a transient impairment in cell cycle progression. Int J Cancer 2012;130:509–20.
  6. Matsumoto T. Crystal structures of MKK4 kinase domain reveal that peptide substrate binds to an allosteric site and induces an auto-inhibition state. Biochem Biophys Res Commun 2010;400:369–73.
  7. Berman HM, Battistuz T, Bhat TN, Bluhm WF, Bourne PE, Burkhardt K, et al. The protein data bank. Acta Crystallogr Sect D: Biol Crystallogr 2002;58:899–907.
  8. MacroModel. Schrodinger LLC, New York: USA; 2011.
  9. LigPrep. Schrodinger LLC, New York: USA; 2011.
  10. Friesner RA, Murphy RB, Repasky MP, Frye LL, Greenwood JR, Halgren TA, et al. Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J Med Chem 2006;19:6177–96.
  11. Qikprop. Schrodinger LLC, New York, USA; 2011.
  12. Irwin JJ, Sterling T, Mysinger MM, Bolstad ES, Coleman RG. ZINC: a free tool to discover chemistry for biology. J Chem Inf Model 2012;52:1757–68.
  13. Lipinski CA, Lombardo F, Dominy BW, Feeney PJ. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Delivery Rev 2001;46:3–26.
  14. Shivakumar D, Williams J, Damm W, Shelley J, Sherman W. Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the opls force field. J Chem Theory Comput 2010;6:1509–19.
  15. Branka A. Nose-hoover chain method for nonequilibrium molecular dynamics simulation. Phys Rev 2000;61:4769–73.
  16. Martyna GJ, Tobias DJ, Klein ML. Constant pressure molecular dynamics algorithms. J Chem Phys 1994;101:4177-89.
  17. Vincent K, van G, Philippe HH. A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J Comput Chem 2001;22:501–8.
  18. Kinnings SL, Jackson RM. Binding site similarity analysis for the functional classification of the protein kinase family. J Chem Inf Model 2009;49:318–29.
  19. Liao JJ. Molecular recognition of protein kinase binding pockets for the design of potent and selective kinase inhibitors. J Med Chem 2007;50:409–24.
  20. Alessi DR, Cuenda A, Cohen P, Dudley DT, Saltiel AR. PD 098059 is a specific inhibitor of the activation of mitogen-activated protein kinase in vitro and in vivo. J Biochem 1995;270:27489–94.
  21. Gum RJ, McLaughlin MM, Kumar S, Wang Z, Bower MJ. Acquisition of sensitivity of stress-activated protein kinases to the p38 inhibitor, SB 203580, by alteration of one or more amino acids within the ATP binding pocket. J Biol Chem 1998;273:15605–10.
  22. Kim N, Park J, Gadhe CG, Cho SJ, Oh Y. A protoberberine derivative HWY336 selectively inhibits MKK4 and MKK7 in mammalian cells: the importance of activation loop on selectivity. PLoS One 2014;9:e91037. Doi:10.1371/ journal. pone.0091037. [Article in Press]
  23. Scheiner S, Kar T, Gu Y. The strength of the C H O hydrogen bond of amino acid residues. J Biol Chem 2001;276:9832–7.
  24. Bissantz C, Kuhn B, Stahl M. A medicinal chemist’s guide to molecular interactions. J Med Chem 2010;53:5061–84.