Molecular docking, 3D QSAR and dynamics simulation studies of imidazo-pyrrolopyridines as janus kinase 1 (JAK 1) inhibitors

Accepted Manuscript Title: MOLECULAR DOCKING, 3D QSAR AND DYNAMICS SIMULATION STUDIES OF IMIDAZO-PYRROLOPYRIDINES AS JANUS KINASE 1 (JAK 1) INHIBITORS Author: Ramesh itteboina Srilata Ballu Sree Kanth Sivan Vijjulatha Manga PII: DOI: Reference:

S1476-9271(15)30263-2 http://dx.doi.org/doi:10.1016/j.compbiolchem.2016.04.009 CBAC 6543

To appear in:

Computational Biology and Chemistry

Received date: Revised date: Accepted date:

28-11-2015 16-4-2016 26-4-2016

Please cite this article as: itteboina, Ramesh, Ballu, Srilata, Sivan, Sree Kanth, Manga, Vijjulatha, MOLECULAR DOCKING, 3D QSAR AND DYNAMICS SIMULATION STUDIES OF IMIDAZO-PYRROLOPYRIDINES AS JANUS KINASE 1 (JAK 1) INHIBITORS.Computational Biology and Chemistry http://dx.doi.org/10.1016/j.compbiolchem.2016.04.009 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

HIGHLIGHTS ♣ CoMFA, CoMSIA models were developed for a series of 43 Imidazo- pyrrolopyridine derivatives. ♣ Molecular docking method was used to analyze possible interactions between receptors and the compounds. ♣ Results of the docking studies and molecular dynamics simulations complement each other. ♣ Seven derivatives as potential candidates of JAK1 inhibitors with good predicted activities were designed.

GRAPHICAL ABSTRACT

MOLECULAR DOCKING, 3D QSAR AND DYNAMICS SIMULATION STUDIES OF IMIDAZO-PYRROLOPYRIDINES AS JANUS KINASE 1 (JAK 1) INHIBITORS Ramesh itteboina, Srilata Ballu, Sree Kanth Sivan, Vijjulatha Manga,*

Janus kinase 1 (JAK 1) plays a critical role in initiating responses to cytokines by the JAK signal transducer and activator of transcription (JAK-STAT). This controls survival, proliferation and differentiation of a variety of cells. Docking, 3D quantitative structure activity relationship (3D-QSAR) and molecular dynamics (MD) studies were performed on a series of Imidazopyrrolopyridine derivatives reported as JAK 1 inhibitors.

MOLECULAR DOCKING, 3D QSAR AND DYNAMICS SIMULATION STUDIES OF IMIDAZO-PYRROLOPYRIDINES AS JANUS KINASE 1 (JAK 1) INHIBITORS Ramesh itteboina, Srilata Ballu, Sree Kanth Sivan, Vijjulatha Manga,*

Department of Chemistry, University College of Science, Osmania University, Hyderabad – 07, INDIA. [email protected] Molecular Modeling and Medicinal Chemistry Group, Department of Chemistry, University College of Science, Osmania University, Hyderabad – 07. Phone: 09866845408.

1

ABSTRACT Janus kinase 1 (JAK 1) plays a critical role in initiating responses to cytokines by the JAK - signal transducer and activator of transcription (JAK-STAT). This controls survival, proliferation and differentiation of a variety of cells. Docking, 3D quantitative structure activity relationship (3D-QSAR) and molecular dynamics (MD) studies were performed on a series of Imidazo-pyrrolopyridine derivatives reported as JAK 1 inhibitors. QSAR model was generated using 30 molecules in the training set; developed model showed good statistical reliability, which is evident from r2ncv and r2loo values. The predictive ability of this model was determined using a test set of 13 molecules that gave acceptable predictive correlation (r2Pred) values. Finally, molecular dynamics simulation was performed to validate docking results and MM/GBSA calculations. This facilitated us to compare binding free energies of cocrystal ligand and newly designed molecule R1. The good concordance between the docking results and CoMFA/ CoMSIA contour maps afforded obliging clues for the rational modification of molecules to design more potent JAK 1 inhibitors.

KEY WORDS CoMFA CoMSIA PLS JAK SP MD

2

1. INTRODUCTION Intracellular signaling of numerous cytokines are mediated by Janus protein tyrosine kinases (JAK 1, JAK 2, JAK 3, and TYK 2); these play critical roles in a variety of biological processes, including hematopoiesis and the regulation of immune and inflammatory responses (Song and Schindler, 2010; Shuai and Liu, 2003; Ghoreschi et al., 2009) . The intracellular domain of cytokine receptors is activated upon ligand induced receptor homo or heterodimerization; that results in the instant phosphorylation of tyrosine residues of cytokine receptor in the cytoplasmic domain were JAKs are associated with the cytoplasmic motifs. These phosphotyrosines then serve as docking sites for the cytoplasmic signal transducer and activator of transcription (STAT) proteins. Upon JAK-mediated tyrosine phosphorylation, the STATs dimerize and translocates into the nucleus where they regulate transcription of specific target genes (Levy and Darnell, 2002; Kisseleva et al., 2002). Biochemical and genetic studies have shown that JAK 1 is the most broadly used JAK, as it is involved in the signaling of the gamma common (γc), beta common (βc), gp130, type I and type II interferon, IL-6, and IL-10 subfamilies of cytokines (Song and Schindler, 2010; Schindler et al., 2007). JAK/STAT pathway transmits signals from the cell membrane to the nucleus in response to cytokines and extracellular growth factors. Activation of this pathway has been found in certain types of human tumors like gliomas whose most common site is brain. Gliomas is a kind of tumor that starts in the brain or spine. It has been proved that JAK/STAT pathway play a role in the progression of human gliomas. Correlation between the JAK/STAT pathway in human gliomas and patients’ prognosis has been investigated by Yanyang et.al. using western blotting analysis and immunohistochemical staining(Yanyang et al., 2011). JAK1, JAK2, and STAT-3 phosphorylation in both human and murine glioma cells is effectively blocked by AZD1480 as a potential therapeutic for the treatment of human glioblastmo (Braden et. al 2011). Curcumins have been reported to inhibit JAK1 2/STAT3 tyrosine-phosphorylation in a dose-dependent fashion in murine glioma cell lines (Weissenberger et al., 2010) asserting JAK1 as an important target for cancer therapy. Computational drug design approaches are vastly employed in the development and optimization of inhibitors. A detailed study of molecular interactions of JAK 1inhibitors with the receptor is utilized to design novel molecules with enhanced anticancer activity. 3

Our primary objective was to study the structural requirements for JAK 1 inhibitors and to design novel molecules. In the present article, we report receptor based 3D-QSAR (Caflisch,1997; Cho et al., 1996; Sippl et al., 2001; Sippl , 2000) studies using CoMFA (Cramer 3rd et al., 1988; Cramer 3rd et al., 1988) and CoMSIA (Klebe et al.,1994) methodologies on Imidazo-pyrrolopyridine derivatives. Partial least square (PLS) (Wold et al., 1993) based statistical analysis was carried out on 43 molecules to identify the correlation. Besides, a combined computational approach, including the docking analysis and molecular dynamics (MD) simulations were also employed to elucidate the probable binding modes of these antagonists at the binding site. The good concordance between the 3D contour maps and the docking results provided insight to identify several key features of the binding mechanism for these Imidazo-pyrrolopyridine derivatives.

2. METHODOLOGY The crystal structure of Janus kinase1 (PDB id: 3EYG) (Williams et al., 2009) was downloaded from the protein databank. GLIDE 5.6 (Schrödinger LLC (2010) Glide, Version 5.6. New York, NY) was used for molecular docking studies. The protein was prepared using protein preparation wizard applying the default parameters, that includes adding hydrogen’s, assigning correct bond orders, filling missing atoms and residues using PRIME, hydrogen bond optimization and minimization. Grid was generated around the active site of the JAK1 with receptor van der Waals scaling for the nonpolar atoms as 0.9 (Friesner et al., 2004). A set of 43 known JAK 1 inhibitors (Table 1) was selected from the literature (Janusz et al., 2012; Mark et al.,2012). These were built using Maestro build panel and prepared by Ligprep application in the Schrödinger 2010 suite. The molecular docking of the 43 molecules into the generated grid was performed by using the standard precision (SP) docking mode (Friesner et al., 2004). The cocrystallized ligand was redocked and its root mean square deviation (RMSD) was calculated to validate the docking process. The dock poses of all the molecules were analyzed for their hydrogen bond interaction with the active site residue. The best dock pose for each molecule was chosen for CoMFA and CoMSIA analysis without further alignment, docked based alignment is shown in Fig.1. The molecules were imported into SYBYL-X 2.1 (Sybyl (2013) Tripos Inc, St. Louis (MO)) molecular modeling program package, Gasteiger4

Huckel charges (Gasteiger and Marsili,1980) were assigned. A standard protocol was followed for calculating steric, electrostatic, hydrophobic, acceptor and donor parameters required for CoMFA and CoMSIA analyses using Tripos force fields (Fathima et al., 2012). A partial least squares (PLS) regression was used to generate a linear relationship that correlates changes in the computed fields with changes in the corresponding experimental values of biological activity (pKi) for the data set of ligands. The data set was divided into a training set randomly consisting of 30 molecules covering the complete activity range and a test set of 13 molecules incorporating at least 70% of activity range. Biological activity values of ligands were used as dependent variables in a PLS statistical analysis. The column filtering value (s) was set to 2.0 kcal mol-1 to improve the signal-to-noise ratio by omitting those lattice points whose energy variation were below this threshold. Crossvalidation was performed by the leave-one-out (LOO) procedure to determine the optimum number of components (ONC) and the coefficient r2loo. Table 1: Structures of the JAK 1 Inhibitors

R N N

Expt. pKi 9.3010

N H Pred. pKi ( CoMFA) 9.382

Pred. pKi (CoMSIA) 9.366

Glide Score (kcal/mol) -10.044

9.3979

9.372

9.413

-9.592

7.6197

7.613

7.569

-8.882

7.0555

7.105

7.771

-9.063

7.8239

7.808

7.716

-8.689

N

Molecule 1

R H3C N

N O

2 N

N O

3 N

4t N

5 N

5

6

7.2076

7.273

7.305

-8.892

8.6197

8.656

8.495

-8.950

7.7212

7.781

7.645

-9.342

7.5528

8.224

7.917

-9.068

8.0177

8.048

7.668

-9.117

7.3010

7.400

7.412

-9.544

7.4948

7.464

7.479

-9.939

7.3665

7.272

7.376

-9.041

7.6575

7.633

7.719

-9.779

8.0809

8.074

7.975

-9.874

8.7695

8.727

8.784

-9.147

7.4814

7.371

7.431

-8.837

7.3010

7.472

7.397

-9.939

N

7

N

8t

O N

9t

CH3

O

CH3

10t

O

O S

11 CH3 O

12

CH3 N CH3

O

13

CH3 N

H3C

O

14

CH3 N

15

N N

16

O

O S CH3

F

17

F F

18

N

6

R

N N

19t

R

Expt. pKi 8.2146

Pred. pKi ( CoMFA) 8.252

Pred. pKi (CoMSIA) 8.082

Glide Score (kcal/mol) -9.241

8.2365

8.034

8.254

-9.972

8.0087

7.898

8.108

-9.891

O

6.5528

6.555

6.450

-9.579

CH3

7.9586

8.100

7.823

-8.960

8.0315

8.015

8.100

-9.641

8.3467

8.393

8.305

-9.957

8.0222

8.037

7.970

-9.235

7.8860

7.935

8.020

-9.163

Pred. pKi (CoMSIA)

Glide Score (kcal/mol)

CH3 NH

O

CH3 CH3

O

H

20

N

H

N

H

H

21

N H

N

Table 1(continued) Molecule

22 N

23t NH

CH3 O

24 NH

25

NH

26

NH

N

O

CH3

O CH3

27

N

O

CH3

O H3C

R N

N

N

Table 1(continued) Molecule

R

Expt. pKi

N H

Pred. pKi ( CoMFA)

7

28

8.000

8.301

8.164

-9.130

8.8860

8.896

8.927

-8.710

30t

8.8860

8.293

8.091

-8.949

31t

8.5528

9.038

8.736

-8.635

8.9208

8.583

8.486

-9.244

8.0315

8.275

8.121

-9.214

8.6989

8.301

8.490

-9.301

8.8239

8.226

8.449

-7.665

8.2146

8.291

8.415

-9.294

8.5850

8.245

8.321

-9.306

8.2365

8.384

8.178

-9.556

NH

29

N

O

32

F

33

F

34t

NH

NH

F F NH

35t

36

NH

N

NH

CH3

O

37t

NH

O CH3

O

38

OH

8

39t

OH

40

8.1739

8.382

8.176

-9.552

8.7447

8.322

8.400

-9.351

8.1366

8.278

8.393

-9.376

8.3187

8.234

8.219

-9.468

8.2676

8.230

8.216

-9.456

OH

41 OH

42t

43

OH

OH

t=Test set molecule

Fig. 1 Docked base alignment of imidazo - Pyrrolopyridines derivatives

An optimum number of components obtained are then used to derive the final QSAR model using all of the training set molecules with non-cross validation and to get the conventional regression coefficient (r2). To validate the CoMFA and CoMSIA derived models, the predictive ability of the test set molecules (expressed as r2pred) was determined (Fatima et al., 2012). Since the statistical parameters were found to be best for the model from the LOO method, it was employed for further prediction of activity of designed molecules. The designed molecules were also constructed, minimized and docked into the 9

protein active site. ADME (absorption, distribution, metabolism and excretion) properties of designed molecules were evaluated computationally using Qikprop module of Schrodinger (Qikprop 3.4, Schrӧdinger, LLC, New York, NY, 2010). The physically significant descriptors and pharmaceutically relevant properties of organic molecules with compliance to Lipinski’s rule of five are calculated by Qikprop. Molecular dynamics simulations Desmond 3.8 (D. E. Shaw Research, DESMOND (Version 3.8), New York, 2014) was used for performing MD simulations to understand the conformational changes on the binding interaction and stability of docked molecule. This was compared with the selective binding modes of the crystal structure ligand (PDB ID: 3EYG) and newly designed molecule R1 (Table 2) utilizing their 3D structures obtained from docking results. The coordinates for MD simulations have been achieved from the docking results. The two complexes were soaked using simple point charge (SPC) water molecules (Berendsen et al., 1987) in an orthorhombic box of 15 Å ×15 Å ×15 Å size. The OPLS 2005 force field (Kaminski and Friesner, 2001) was employed, one chlorine ion was added randomly to balance, the net charges of the system and salt concentration was maintained at 0.15 molL-1 during the entire course of the simulation. The systems were then minimized twice, with restraints on solute initially and without any restraints later using the steepest descent integrator for maximum 1000 iterations, and the limited memory Broyden-FletcherGoldfarb-Shanno (LBFGS) algorithms until the value converged to 5 kcal mol−1 Å−1. The systems were then relaxed before extensive MD simulations using a relaxation protocol composed of small (12 and 24 ps) MD Simulations. During MD simulations, the physiological temperature was maintained at 300K using NoseHoover thermostat approach with a relaxation time of 1ps and the pressure was maintained at 1.01325 bars by employing the Martyana-Tobaisklein barostat approach using coupling style isotropic with a relaxation time of 2 ps. The cutoff distance for computing non-bonded interactions was truncated at 9 Å and long range electrostatic interactions were calculated using particle-mesh Ewald (PME) method (Strahan et al., 1998) with Ewald tolerance of 1e-09. The SHAKE algorithm (Strahan 10

et al., 1983) applied to all bonds involving H-bonds. High quality trajectories were recorded at 4.8 ps during MD simulations and the recording interval energy was set at 1.2 ps. Finally, 5ns MD simulations were performed with a two fs time step at NPT canonical ensemble within the periodic boundary conditions. The root mean square deviations (RMSD) concerning simulation time, interaction during simulation were examined and binding free energy calculations were performed for investigating the energetic contribution of protein-ligand binding affinities. Table 2: Glide Score and Predicted Activity of Newly Designed Molecules R

R

2

R

N

1

3

O

N N

N

N H

Molecule

R1

R2

R3

R1

CH3

OH

C6H5

9.719

9.847

-9.812

R2

C2H5

OH

C6H5

9.597

9.862

-9.597

R3

C2H5

OH

C5H9

9.550

9.648

-9.348

R4

CH3

OH

C6H11

9.515

9.541

-9.294

R5

CH3

OH

C5H9

9.498

9.662

-9.141

R6

CH3

OH

C4H7

9.429

9.613

-9.433

R7

C2H5

OH

C4H7

9.392

9.615

-9.436

Pred. pKi ( CoMFA)

Pred. pKi (CoMSIA)

Glide Score (kcal/mol)

3. RESULTS AND DISCUSSION Evaluation of the accuracy of docking process is determined by measuring how closely the lowest energy pose (binding conformation) predicted by the object scoring function (Glide score), resembles an experimental binding mode as determined by X-ray crystallography. In the present study Glide (SP) docking was validated by redocking the crystal structure ligand into the binding site of JAK 1. An exquisite agreement between the localization of the inhibitor upon docking and the position of the ligand in the crystal structure was observed. The root mean square deviation between the docked

11

conformation and the observed X-ray crystallographic conformation of molecule 3{(3R,4R)-4-methyl-3-[methyl(7H-pyrrolo[2,3-d] pyrimidin-4-yl) amino] piperid-1-yl}-3oxo propane nitrile equaled 0.30Å (Fig.2), suggesting the reliability of Glide docking in reproducing the experimentally observed binding mode for Janus kinase 1 inhibitors.

Fig. 2 Superimposition of crystal structure pose (cyano) on dock pose (red) of co-crystallized ligand. The rms deviation is 0.3069 A0

The analysis of 3D-QSAR CoMFA and CoMSIA was carried out using imidazopyrrolopyridine derivatives reported as inhibitors of Janus kinase 1. Molecules having inhibitory activity against Janus kinase 1 enzyme with exact Ki values were selected. A total of 43 molecules were used for derivation of the model, these were divided into a training set of 30 and test set of 13 molecules. The CoMFA and CoMSIA statistical analysis is summarized in (Table 3), Statistical data show q2loo of 0.504 for CoMFA, 0.518 for the CoMSIA models, r2ncv of 0.948 for CoMFA and 0.951 for CoMSIA indicating a good internal predictive ability of the models. Predictive ability of the models was tested by taking 13 molecules that were excluded from the model derivation as the test set. The predictive correlation coefficient r2pred of 0.52 for CoMFA and 0.53 for the CoMSIA models indicate a good external predictive ability of the models. The experimental and predicted activity from CoMFA and CoMSIA model is given in Table 1. A scatter plot of experimental vs. predicted pKi values is represented in (Fig.3). The basic intension of performing 3D QSAR analysis is to obtain a model that predicts the activity of molecules accurately and also help in developing novel molecules by structural modification, CoMFA and CoMSIA contour maps help in determining structural requirements for receptor binding and biological effect. The CoMFA/CoMSIA 12

results

were

graphically

interpreted

by

field

contribution

maps

using

the

“STDEV*COEFF” field type. Table 3: PLS Result Summary Statistical parameters CoMFA q2looa 0.504 Number of molecules in training set 30 Number of molecules in test set 13 ONCb 5 SEEc 0.162 r2 d 0.948 Fratioe 87.050 r2predf 0.52 Fraction of fields contributions Steric 71.2% Electrostatic 28.8 Hydrophobic -Acceptor -Donor a – Cross-validation correlation coefficient by leave one out method, b – Optimum number of components c – Standard error of estimate d – Conventional correlation coefficient e – Fisher test value f – Cross-validation correlation coefficient on the test set.

11.7% 14.8% 29.8% 17.5% 17.5%

Scatter plot of Expt.pki vs Pred pki - CoMSIA

Scatter plot of Expt.pki vs Pred pki-CoMFA

10

Predicted - pki

10

Predicted-pki

CoMSIA 0.518 30 13 6 0.160 0.951 74.333 0.53

9 8 y = 0.9478x + 0.4211

7

9

8

y = 0.951x + 0.3954

7

R2 = 0.951

R2 = 0.9477

6

6

6

7

8

Experimental-pki

9

10

6

7

8

9

10

Experimental - pki

Fig.3 Scatter plot of predicted vs. experimental pKi values (The test set is represented as squares and training set represented as triangles)

Fig.4 and 5 shows the contour maps derived from the CoMFA and CoMSIA models respectively. The most potent analogue, molecule 2 was embedded in the map (a, 13

c), while least active molecule 22 was embedded in the map (b, d) to visualize steric and electrostatic regions around the molecules effecting activity. The areas of green indicate a favourable contribution to potency while yellow region the converse. The blue and red contours indicate positive electrostatic charge and negative charge regions respectively. Substitutions on piperidine attached to imidazole ring orients towards favoring green region indicate a substitution at this position with a bulky group leads to increase in the activity, as in the case of most active molecules 1 and 2. In the case of lowest active molecule 22, substitutions on cyclohexane attached to imidazole ring orients towards disfavored region explaining the lower activity of the molecules. Electrostatic contour are covering the substitution attached to core indicating charged groups at this region will affect the activity.

(a)

(b)

(c) (d) Fig. 4 CoMFA strict standard deviation (S. D.* coefficient) contour maps illustrating steric and electrostatic features in combination with molecule 2 and 22. (a, b) Green contours show favorable bulky group substitution at that point while yellow regions show disfavorable bulky group for activity. (c, d) Red contours indicate negative charge favoring activity. Whereas blue contours indicate positive charge favoring activity.

14

(a)

(c) Fig. 5 CoMSIA S.D.

*

(b)

(d) coefficient contour maps illustrating steric and electrostatic features in

combination with molecule 2 and 22. (a,b) The green contour indicates a sterically favored region; yellow maps calls for a reduction of this potential to improve activity. (c, d) Blue indicates a positive charge preferred region to improve activity

Fig.6, shows the hydrophobic contour maps derived from the CoMSIA for molecule 2 and 22. Hydrophobic yellow region is seen near the imidazole ring and a hydrophilic white region near the piperidine ring. Substitution at the imidazole ring with a hydrophobic group and a substitution on the piperidine ring with a hydrophilic group should enhance the activity. Fig.7 (a, b) represents hydrogen bond donor contours, were cyan and purple contour represents favorable and disfavorable hydrogen bond donor groups respectively. Fig.7 (c, d) represents hydrogen bond acceptor contours, were magenta and red contours represent favorable and unfavorable H-bond acceptor groups respectively. Magenta contour near the carbonyl group of the piperidine ring suggests that the carbonyl group is an important substituent for the molecule. In Fig.7a, the CH2 group of the piperidine ring is directed towards the cyan contour expressing a donor 15

favored group at this position will enhance the activity. Fig.7c shows a red contour (an acceptor disfavored) at the same position where the blue patch appeared in the donor contour suggesting consistency in the analysis.

(a)

(b)

Fig. 6 CoMSIA S.D.*coefficient contour maps illustrating hydrophobic features in combination with molecule 2 and 22.(a, b) The yellow contour for hydrophobic favored region, white indicates the hydrophilic favored region

(a)

(c)

(b)

(d)

Fig. 7 CoMSIA S.D.* coefficient contour maps illustrating acceptor and donor features in combination with molecule 2 and 22. (a, b) The cyan contour for H-bond donor group increase activity, purple indicates the disfavored region. (c,d) The magenta contour for H-bond acceptor group increase activity, red indicates the disfavored region

16

Structural modification for improving the inhibitory activity of JAK 1 was performed on reference molecule 2 based on structural requirements depicted in figure 8. Based on the hydrophobic contours the second position of the imidazole ring was substituted with hydrophobic groups like methyl and ethyl to improve the potency of the molecules. A hydroxy group has been substituted at the 3rd position in the piperidine ring, to satisfy the structural requirements based on the hydrogen bond acceptor and donor contour. Acetonitrile side chain has been substituted with bulky hydrophobic groups (phenyl, cyclopentyl and cyclobutyl) as analyzed by the contour maps. The newly designed molecules were docked into the protein active site and the dock poses were used to predict the activity by applying the 3D-QSAR model. Figure 9 shows the molecular interaction between inhibitors and receptor obtained from docking studies. The hydrophobic region of the receptor is also depicted in the figure 9, to clarify further the role of substitution modification during lead optimization. Bulky hydrophobic group substituted at the place of acetonitrile side in molecule R1 penetrates into the hydrophobic pocket increasing the affinity of the molecule compared to molecule 2. It also maintained the hydrogen bond interaction as the reference molecule 2. The designed molecules showed a better dock score and predicted activity compared to the existing

Hyd reg roph ion obic

Fav or e

d

ones. d Steric,Electrostatic,H-Bon Donor Favored region Ster i Acc c,Electr epto o r Di static F sfav ored avored a regi on nd N N i ob N ph red o r O yd avo ed N c,H Disf avor i t F sta r tro ono ptor c e e ,El d D Acc N N ric H Ste -Bon ond B c,H Hd n a ion reg

Fig. 8 Structure requirements for binding and inhibitory activity of Indole hydrazone

17

(a)

(b) Fig. 9 Docked pose of (a) molecule 2 and (b) newly designed molecule R1 in the protein active site, showing hydrogen bond interaction (yellow lines) with Glu 957 and Leu 959 residues. Pink surface represents hydrophobic cavity in the receptor.

Molecular dynamics simulations Molecular Dynamics (MD) simulation is a useful implicit mechanism to understand molecular interactions possible in a protein-ligand complex that is responsible for its stability. To gauge and compare the stability and molecular interactions of newly designed 3D QSAR based molecule, R1 and cocrystal ligand, MD simulations were performed for 5ns. Fig.10 shows the root mean square deviations (RMSD) of protein backbone with simulation time for each complex from 18

the initial structure. Initially, RMSD reached a value of 1.782 Å from 0.731 Å during 0-500ps and then retained between 1.387-2.13 Å throughout the simulation up to 5ns for cocrystal ligand complex whereas RMSD of molecule R1 complex reached a value of 1.828 Å initially and then retained between 1.261-2.381 Å during the entire 5ns simulation time. The averaged RMSD of the crystallized complex was found to be around 1.545 Å, whereas for R1- protein complex was around 1.601 Å indicating similar stability for both the complexes. Change in binding energies during the simulation period using MD simulation trajectory was calculated.

Fig.10 Root mean square deviations (RMSD) of protein back bone as function of simulation time of each complex with respect to initial structure. Blue colour shows RMSD of protein back bone of crystal structure ligand where as Red colour shows RMSD of protein back bone with molecule R1.

Fig.11 shows changes in binding energies (dG) with simulation time for both the protein-ligand complexes. From Fig.11 it is observed that both the complexes showed stable dG during simulation time, dG bind values ranged from -74.692 to 92.249kcal/mol for crystal structure ligand and -61.630 to -84.361kcal/mol for molecule R1. Average dG bind for the crystal structure ligand was found to be 81.309kcal/mol and for molecule R1 it was around -77.813kcal/mol. A similar trend of binding energies over a period of simulation for both the complexes shows their relative stability, strengthening the assumption that molecule R1 can act as a potent inhibitor against JAK 1. Analysis of protein-ligand interactions is the most 19

significant part in MD simulations; it illustrates the changes in the mode of binding during simulations. Protein-ligand interactions or contacts are mainly of Leu 881, Val 889, Ala 906, Lys 908, Glu 957, Phe 958, Leu 959, Leu 1010 residues. Glu 957, Leu 959 were considered as the critical residue for ligand binding. Fig.12 shows protein-ligand contact interaction over trajectory. The Y-axis is normalized over the course of the trajectory. From Fig.12(b) it is evident that the newly designed molecule R1 has not only retained all the interactions as that of cocrystal ligand but also showed some new interactions with residues like Arg 879, Val 989, Met 956 and Asn 1008. Moreover molecule R1 maintained the hydrogen bond interactions with same Glu 957 for 99% and Leu 959 for 89% of the simulation time. The cocrystal ligand maintained hydrogen bond interactions with Glu 957 for about 100% and Leu 959 for about 86% of the simulation time (Fig.13(b)). Further molecule R1 showed enhanced interactions with Asn 1008 than crystal structure ligand. Docking and MD simulations showed similar interactions in case of the crystal structure ligand, but whereas in molecule R1 Glu 957 and Leu 959 interactions were retained during the course of simulation, but interactions with Lys 908 was lost and a new hydrogen bond interaction was envisaged with Asn 1008 for more than 47% of the simulation time. The MD simulation study of both the complexes sufficiently explains the binding modes of inhibitors and enhances the acceptability of 3D QSAR model for the design of new molecules.

Fig.11 Change in Binding energy (dG) with respect to simulation time. Red colour represents protein complexed with crystal structure ligand and blue colour represents protein complexed with molecule R1.

20

Fig.12 Protein ligand contact interaction over trajectory with respect to (a) Co-crystal ligand and (b) Molecule R1.

Fig.13 Ligand interaction diagram of (a) Co-crystal ligand and (b) newly designed molecules R1 with protein active site during the course of dynamic simulation showing average interaction percentage with active site residues.

Predicted ADME properties Newly designed molecules were analyzed for drug-likeness by assessing their physiochemical properties and by applying Lipinski's rule of five (Table 4). The Lipinski's rule for the drug-like molecules, states that the molecule should have molecular weight <650 Daltons, H-bond donors <5, H-bond acceptors <10, and a log P of <5. For the designed molecules, the partition coefficient (QPlogPo/w) and water solubility 21

(QPlogS) is critical for estimating the absorption and distribution of drugs within the body, which ranged between 1.734-2.631 and 3.544-4.425, respectively. Crossing the blood-brain barrier (BBB), which is a prerequisite for the entry of drugs into CNS, was found to be in the acceptable range (1.15- to 1.292), indicating that the molecules can be considered for further development in treating gliomas. Caco-2 cell permeability (QPPCaco), a model governing gut-blood barrier, ranged from 131.998 to 197.216. MDCK cell permeability (QPPMDCK), a model that mimics the blood-brain barrier, ranges from 106.751 to 155.667. Further, the predicted percentage human oral absorption, ranges from 76.064 to 82.707. All these pharmacokinetic parameters were found to be within the acceptable range Table 4 : ADME Properties of Newly Designed Molecule Molecule a

QPlogPo/wb

QPlogSc

QPPCacod

QPlogBBe

QPPMDCKf

R1

2.425

-4.334

131.998

-1.292

106.751

%Human oral absorptiong 79.099

R2

2.631

-4.202

141.000

-1.281

114.207

80.815

R3

2.508

-4.425

197.216

-1.188

155.667

82.707

R4

2.212

-4.135

144.382

-1.200

113.869

78.546

R5

2.067

-4.117

166.292

-1.178

125.500

78.799

R6

1.734

-3.544

150.306

-1.150

118.025

76.064

R7

2.249

-3.990

161.553

-1.217

127.328

79.639

a. Newly designed molecules. b. Predicted octanol/water partition coefficient log P (Acceptable range-2.0 to 6.5) c. Predicted aqueous solubility in mol/L (Acceptable range -6.5 to 0.5) d. Predicted Caco cell permeability in nm/s (Acceptable range :< 25 is poor and >500 is great) e. Predicted Blood Brain Barrier permeability (Acceptable range-3 to 1.2) f. Predicted apparent MDCK cell permeability in nm/s (Acceptable range in nm/s (Acceptable range :< 25 is poor and >500 is great) g. Percentage of human oral absorption (Acceptable range: <25 is poor and >80% is high.

CONCLUSIONS Molecular docking based 3D-QSAR and MD simulation studies are widely used tools for understanding binding modes of the molecule to the enzymatic receptors and rationalize the structural requirements for the inhibitory activity of the molecule. CoMFA and CoMSIA methodologies were used to build models for JAK 1 inhibitory activity of imidazo-pyrrolopyridines and C-2 methyl imidazo-pyrrolopyridine derivatives. The generated models have the statistical reliability that is evident from high r2ncv, q2loo and 22

r2pred values for all the models. Based on detailed contour map analysis, improvement in JAK 1 binding affinity can be achieved by substitutional modification at the imidazole and piperidine ring. The designed molecules based on these parameters showed better predicted activity than reference molecules, this indicates QSAR models generated have a good predictive ability to design potent inhibitors. The critical role of Glu 957, Leu 959 in the binding of inhibitors to JAK 1 revealed from docking studies was validated through MD simulations. Further MD simulations of protein-ligand complex also revealed that JAK 1 inhibitors may undergo a conformational change to form a hydrogen bond with Asn 1008 retaining hydrogen bond interactions with Glu 957 and Leu 959. These Newly designed molecules act as leads with required pharmacokinetics that are used to obtain a greater number of molecules for clinical studies.

ACKNOWLEDGMENTS We greatly acknowledge Tripos Inc, USA and Schrödinger LLC, New York for providing us the software. This research was made possible through grants from DST-SERB (SB/EMEQ-004/2013), CSIR (01/ (2436)/10/EMR-II), UGC (42233/2013(SR)), and UGC/UPE/FAR/OU/2014 New Delhi, India The authors RI and SB would like to acknowledge financial support from UGC for research fellowships. We wish to express our gratitude to Department of Chemistry, Osmania University for providing facilities to carry out the research work.

REFERENCES Andersen, H.C., 1983. Rattle: A “velocity” version of the shake algorithm for molecular dynamics

calculations.

J.

Comput.

phys.

52,

24-34.doi:http://dx.doi.org/doi:

10.1016/0021-9991(83)90014-1. Berendsen, H.J.C., Grigera, J.R., Straatsma, T.P., 1987. The missing term in effective pair potentials. J. Phys. Chem. 91, 6269–6271. doi:http://dx.doi.org/doi: 10.1021/j100308a038. Braden, C.M., Ma J.Y., Catherine, P.L., Gillespie, G.Y., Yu, H., Zheng, Y., Susan, E.N., Dennis, H., Etty N. B., 2011. Therapeutic Potential of AZD1480 for the Treatment of Human Glioblastoma. Mol Cancer Ther; 10(12), 2384-2393. DOI: 10.1158/1535-7163.MCT-110480 23

Caflisch, A., Ehrhardt, C., 1997. Structure-based combinatorial ligand design, in: P. Veerapandian (Ed), Structure based drug design, Marcel Dekker, New york, pp. 541- 558. Cho, S.J., Garsia, M.L.S., Bier, J., Tropsha, A., 1996. Structure- Based Alignment and Comparative Molecular Field Analysis of Acetylcholinesterase Inhibitors. J. Med. Chem. 39, 5064–5071. doi:http://dx.doi.org/ doi: 10.1021/jm950771r. Cramer, R.D.3rd., Patterson, D.E., Bunce, J.D., 1988.

Comparative molecular field

analysis(CoMFA). 1. Effect of shape on binding of steroids to carrier proteins. J.Am. Chem. Soc.110, 5959-5967. doi:http://dx.doi.org/ doi: 10.1021/ja00226a005. Cramer, R.D.3rd., Patterson, D.E., Bunce, J.D., 1988. Cross validation, bootstrapping, and partial least squares compared with multiple regression in conventional QSAR studies. Quant. Struct. Act. Relat. 7, 18-25. doi:http://dx.doi.org/doi: 10.1002/qsar.19880070105. Fatima, S., Raju B., Sree Kanth S., Vijjulatha M., 2012. Molecular docking and 3D-QSAR studies on inhibitors of DNA damage signaling enzyme human PARP-1. Journal of receptors and signal transduction. 32(4), 214224.doi:http://dx.doi.org/doi:10.3109/10799893.2012.693087. Friesner, R.A., Banks, J.L., Murphy, R.B., Halgren, T.A., Klicic, J.J., Mainz, D.T., Repasky, M.P., Knoll, E.H., Shelley, M., Perry, J.K., Shaw, E.D., Francis, P., Shenkin, P.S., 2004. Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J. Med. Chem.47, 1739–1749. doi:http://dx.doi.org/doi: 10.1021/jm0306430. Gasteiger, J., Marsili, M., 1980. Iterative partial equalization of orbital electro negativity a rapid access to atomic charges. Tetrahedron.36, 3219-3228. doi:http://dx.doi.org/oi:10.1016/0040-4020 (80)80168-2. Ghoreschi, K., Laurence, A., O’Shea, J.J., 2009. Janus Kinases in Immune Cell Signalling. Immunol.Rev.228, 273−287. doi:http://dx.doi.org/doi: 10.1111/j.1600-065X.2008.00754.x. Janusz, J.K., Wade, B., Richard, J.B., Christine, C., Gauri, D., Hazel, D., Charles, E., Nico, G., Paul, G., Trevor, K.H., Peter, R.H., Marya, L., Christopher, A.H., Adam, J., Tony, J., Jane, R.K., Pawan, B.K., Robert, J.M., Rohan, M., Kyle, M., Jeremy, M., Raman, N., Steven, S., Micah, S., Savita, U., Mark, U., Anne, V.A., Stuart, I.W., Bohdan, W., Mark, Z., 2012. Identification of Imidazo-Pyrrolopyridines as Novel and Potent JAK1 Inhibitors. J. Med. Chem. 55, 5901-5921. doi:http://dx.doi.org/doi: 10.1021/jm300438j. 24

Kaminski, G.A., Friesner, R.A., 2001. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. chem.105, 6474–6487. doi:http://dx.doi.org/doi: 10.1021/jp003919d. Kisseleva, T., Bhattacharya, S., Braunstein, J., Schindler, C.W., 2002. Signalling Through the JAK/STAT Pathway Recent Advantages and Future Challenges. Gene. 285, 1−24. doi:http://dx.doi.org/doi:10.1016/S0378-1119(02)00398-0. Klebe, G., Abraham, U., Mietzer, T., 1994. Molecular similarity indices in a comparative analysis (CoMSIA) of drug molecules to correlate and predict their biological activity. J. Med. Chem. 37, 4130-4146. doi:http://dx.doi.org/doi: 10.1021/jm00050a010. Levy, D.E., Darnell, J.E., 2002. Jr. STATS: Transcriptional Control and Biological Impact. Nat. Rev. Mol. Cell. Biol. 3, 651−662. doi:http://dx.doi.org/doi:10.1038/nrm909. Mark, Z., Rohan, M., Mercedesz, B., Kathy, B., Philippe, B., Wade, S.B., Christine, C., Gauri, D., Jason, D., Peter, S.D., Charles, E., Nico, G., Paul, G., Stefen, G., Chris, H., Emily, J.H., Eric, H., Peter, R.H., Christopher, A.H., Tian, J., Adam, J., Tony, J., Jane, R.K., Michael, F.T.K., Pawan, B.K., Janusz, J.K., Sharada, L., Jiangpeng, L., Marya, L., Zhonghua, L., Patrick, J.L., Robert, J.M., Jeremy, M.M., Rebecca, P., Madelein, R., Scott, S., Steven, S., Micah, S., Savita, U., Mark, U., Anne, V.A., Stuart, I.W., X. Ling, X., Yisong, X., 2012. Discovery and Optimization of C-2 Methyl Imidazo-pyrrolopyridines as Potent and Orally Bioavailable JAK1 Inhibitors with Selectivity over JAK2. J. Med. Chem. 55, 6176-6193. doi:http://dx.doi.org/doi: 10.1021/jm300628c. Schindler, C., Levy, D.E., Decker, T., 2007. JAK-STAT Signaling: FromInterferons to Cytokines. J. Biol. Chem. 282, 20059−20063. doi:http://dx.doi.org/doi: 10.1074/jbc.R700016200. Schrödinger LLC (2010) Glide, Version 5.6. New York, NY. Schrӧdinger LLC (2010) Qikprop, version 3.4, New York,NY. Shaw, D.E., 2014. Research, DESMOND (Version 3.8), New York. Shuai, K., Liu, B., 2003. Regulation of JAK-STAT Signalling in the Immune System. Nat. Rev. Immunol. 3, 900−911. doi:http://dx.doi.org/oi:10.1038/nri1226. Sippl, W., 2000. Receptor-based 3D QSAR analysis of estrogen receptor ligands merging the accuracy of receptor-based alignments with the computational efficiency of ligand-based methods.J.Comput.Aided.Mol.Des.14, 559–572. doi:http://dx.doi.org/doi: 10.1023/a:100811591378. 25

Sippl,W., Contreras, J.M., Parrot, I., Rival, Y.M., Wermuth, C.G., 2001. Structure-based 3D QSAR and design of novel acetylcholinesterase inhibitors. J. Comput. Aided. Mol. Des. 15, 395–410. doi:http://dx.doi.org/doi: 10.1023/a:1011150215288. Song, L., Schindler, C., 2010. JAK-STAT Signalling, Second nd ed., Bradshaw, R.A.,

Dennis

E.A., (eds.). The Handbook of Cell Signalling, Elsevier Inc, San Diego C A. 3 pp. 2041−2048. doi:http://dx.doi.org/doi: 10.1023/a:1011150215288. Strahan, G.D., Keniry, M.A., Shafer, R.H., 1998. NMR structure refinement and dynamics of the K+-[d(G3T4G3)]2 quadruplex via particle mesh Ewald molecular dynamics simulations. Biophys. J. 75, 968–981. doi:http://dx.doi.org/doi: 10.1016/S0006-3495(98)77585-X. Sybyl-X-2.1 version (2013) Tripos Inc., St. Louis (MO) Weissenberger, J., Priester, M., Bernreuther, C., Rakel, S., Glatzel, M., Seifert, V., Kögel, D., 2010. Dietary curcumin attenuates glioma growth in a syngeneic mouse model by inhibition of the JAK1,2/STAT3 signaling pathway. Clin. Cancer. Res. 16, 5781-95. doi: 10.1158/1078-0432.CCR-10-0446. Williams, N.K., Bamert, R.S., Patel, O., Wang, C., Walden, P.M., Wilks, A.F., Fantino, E., Rossjohn, J., Lucet, I.S., 2009. Dissecting specificity in the Janus kinases: the structures of JAK-specific inhibitors complexed to the JAK1 and JAK2 protein tyrosine kinase domains. J. Mol. Biol. 387, 219-232. doi:http://dx.doi.org/doi: 10.1016/j.jmb.2009.01.041. Wold, S., Johansson, A., Cochi, M., 1993. PLS-Partial least squares projection to latent structures, in: H. Kubinyl (ed), 3D-QSAR in drug design: theory, methods and application, ESCOM, Liedin, pp. 523-550. Yanyang, T., Yuexia, Z., Jianfang, F., Yizhan, C., Guoqiang, F., Xiaoxi, T., Boliang W., 2011. Activation of JAK/STAT signal pathway predicts poor prognosis of patients with gliomas. Med.Oncol. 28,15-23. doi: 10.1007/s12032-010-9435-1. Epub 2010 Feb 5.

26