Molecular Simulation Studies on the Binding ... - ACS Publications

Mar 20, 2017 - Molecular Simulation Studies on the Binding Selectivity of Type‑I ... Institute of Bioinformatics and Medical Engineering, School of ...
0 downloads 0 Views 4MB Size
Subscriber access provided by University of Newcastle, Australia

Article

Molecular simulation studies on the binding selectivity of type-I inhibitors in the complexes with ROS1 versus ALK Yuanxin Tian, Yonghuan Yu, Yudong Shen, Hua Wan, Shan Chang, Tingting Zhang, Shanhe Wan, and Jiajie Zhang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00019 • Publication Date (Web): 20 Mar 2017 Downloaded from http://pubs.acs.org on March 21, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Molecular Simulation Studies on the Binding Selectivity of Type-I Inhibitors in the Complexes with ROS1 versus ALK Yuanxin Tian †,1, Yonghuan Yu†,1, Yudong Shen§, Hua Wan‖, Shan Chang ‡,Tingting Zhang†, Shanhe Wan †, Jiajie Zhang *† †

Guangdong Provincial Key Laboratory of New Drug Screening, School of

Pharmaceutical Sciences, Southern Medical University, Guangzhou, 510515, PR China §

College of Food Sciences, South China Agricultural University, Guangzhou, 510642,

PR China ǁ

College of Mathematics and Informatics, South China Agricultural University,

Guangzhou, 510642, PR China ‡

Institute of Bioinformatics and Medical Engineering, School of Electrical and

Information Engineering, Jiangsu University of Technology, Changzhou, 213001, PR China

1

These authors contributed equally to this work.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT: ROS1 and ALK are promising targets of anti-cancer drugs for non small cell lung cancer. Since they have 49% amide acid sequence homology in the kinases domain and 77% identity at the ATP binding area, some ALK inhibitors also showed some significant responses for ROS1 in the clinical trial, such as the type-I binding inhibitor crizotinib and PF-06463922. As a newly therapeutic target, selective ROS1 inhibitor is relatively rare. Moreover, the molecular basis for the selectivity of ROS1 versus ALK still remains unclear. In order to disclose the binding preference towards ROS1 over ALK and to aid the design of selective ROS1 inhibitors, the specific interactions and difference of conformational changes in the dual and selective ROS1/ALK inhibitors systems were investigated by molecular dynamics (MD) simulation and principle component analysis (PCA) in our work. Afterward, binding free energies (MM/GBSA) and binding free energies decomposition analysis indicated that the dominating effect of Van de Waals interaction drives the specific binding process of type-I inhibitor, and residues of P-loop and DFG-motif would play an important role to selectivity. Based on the modeling results, the new designed compound 14c was verified as a selective ROS1 inhibitor versus ALK and SMU-B was a dual ROS1/ALK inhibitor by the kinase inhibitory study. These results are expected to facilitate the discovery and rational design of novel and specific ROS1 inhibitors.

ACS Paragon Plus Environment

Page 2 of 43

Page 3 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

INTRODUCTION Oncogenic drivers have been discovered as drugable targets in patients with non-small cell lung cancer (NSCLC) ,1,2 and tyrosine kinases inhibitors (TKI) have emerged as a particularly successful treatment morality.3 ALK and ROS1 are fusion-type protein tyrosine kinases highly expressed in a variety of tumor cell line, contributing to tumorigenesis or tumor progression.4-6 They regarded as potentially therapeutic targets for cancer therapy.7-9 Since they have 49% amide acid sequence homology in the kinase domain and 77% identity at the ATP binding area, there are several ALK inhibitors that show the high activity against ROS1, such as crizotinib, PF-06463922, and so on.10-14 After crizotinib was approved by FDA for ALK-rearranged lung cancer, the clinical trial for crizotinib in ROS1-rearranged NSCLC patients also displayed significant responses.15,16 With the development of remarkable achievement about the function and signaling pathway about ROS1, more and more evidences support that the ROS1 plays an important role in different malignancies.17 Hence, discovery some small molecules as selective and potent ROS1 inhibitors are still received additional consideration in the medicinal chemistry. In the previous studies, most ROS1 inhibitors have been disclosed during research aiming to other kinases inhibitors, such as the crizotinib/PF-06463922 against ALK and foretinib as c-met inhibitor.11, 18, 19 Since ROS1 is a novel therapeutic target for different malignancies with high homologous with ALK, few inhibitors have been identified with optimal selectivity for ROS1 relative to ALK except the S. H. Lee’s team .20-22 Furthermore, the molecular basis for the selectivity of ROS1 versus ALK

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

still remains unclear, which hinders the rational design for selective ROS1 inhibitor. Recently, B.J Duker and their team23 disclosed the subtle structural differences of their kinases domain and the selective binding site for ROS1/ALK through the systems of the foretinib and cabozantinib as selective inhibitors for ROS1 versus ALK. They considered the dual ROS1/ALK inhibitors (such as ceritinib, PF-06463922 and AZD-3463) almost adopted a typical type-I binding mode (DFG-in conformation), while the ROS1-selective ones (such as foretinib and cabozantinib) exhibit the type-II binding mode (DFG-out conformation). They believed that the type-II binding mode had larger permit accommodation with the movement of Cα-helix and A-loop in an induced-fit mechanism and formed the specific site pocket. However, the ALK-selective inhibitor alectinib exhibits no ROS1 inhibitory activity although it displayed typical type-I binding mode. Moreover, the compound 7c from S. H. Lee team (structure shown in Figure1) was reported to be a highly potent inhibitor with 170 fold selectivity more than that for ROS1 relative to ALK, while its binding mode was also similar to crizotinib and different from the type-II binding mode based on their modeling experiment. What bring out the specific interactions for ROS1 versus ALK for the type-I binding mode on earth? Do the conformational changes occur when the selective compound binding to ROS1 in type-I mode? In order to probe above issues, the compounds 7c (ROS1-selective), alectinib (ALK-selective), and crizotinib (dual ROS1/ALK inhibitor) with type-I binding mode to their receptors were chosen to investigate the mechanism of ligand selectivity for type-I inhibitor to ROS1 versus ALK. Molecular dynamics (MD) simulations and

ACS Paragon Plus Environment

Page 4 of 43

Page 5 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

principle component analysis (PCA) were used to explore the subtle structure differences of binding sites and the conformation changes of type-I binding mode systems. The former was widely used to describe the dynamic behavior of the complex at atomic-level,24-27 the latter was applied to explore the functional dynamics and conformational changes of the four systems.28-32 In addition, the binding free energy calculations and free energy decomposition provided some interesting clues to find some pivotal residues and interaction relative to the selectivity of ROS1 versus ALK, which were in accord with the conformations changes. Our results demonstrated that the active site composed of P-loop, DFG-motif, and Cα-helix would be responsible for the selectivity of 7c to ROS1 versus ALK. Three key residues (Gly24 (1957) from P-loop, Lys47 (1980) from β3-sheet, and Asp169 (2102) from DFG motif) may contribute to the activity and high selectivity for ROS1 by analysis of movement motion and binding free energies. The results provide the basis for achieving selectivity and may facilitate the rational design and development of potential and selective type-I ROS1 inhibitors. Our new compounds 14c and SMU-B designed based on the results and their kinase inhibitory activities proved our results and the design strategy.

2. SYSTEMS AND METHODS 2.1 Initial Structure Preparation The three complex crystal structures were obtained from the Protein Data Bank: selective ALK inhibitor (alectinib, PDB entry:3AOX), dual ALK/ROS1 inhibitor (crizotinib with ALK:2XP2; crizotinib with ROS1:3ZBF). The missing residues were

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

modeled by sybyl7.3 program.33 The 3D structure of complex ROS1 with 7c was obtained by surflex-dock. The docking mode was consistent with the literature.21 The partial atomic charges for the ligands atoms were calculated using the restrained electrostatic potential (RESP)34 protocol after structure optimization with a level of quantum mechanical treatment by B3LYP method at the HF/6-31G* level of the Gaussian 09 suite,35 implemented in AmberTools12.36 The force field parameters for the small molecules generated using the Antechamber program were described by the general amber force field (GAFF).37 The AMBER ff99SB force field38 was used for standard and the ionization state of amino acid residues and waters according to the standard protocol. The treated structures were used as the initial structures of MD simulation. The missing hydrogens of protein were first added using LEaP program from AMBER12.39

Figure 1.The chemical structure of compounds 2.2 Molecular Dynamics Simulation: : The MD simulations were performed using AMBER12 software package.39 Each complex was immersed in a periodic truncated octahedron box of TIP3P40 water molecules with a margin distance of 10.0 Å. Then Cl- or Na+ ions were added to neutralize the systems using LeaP module in AMBER12. Particle Mesh Ewald

ACS Paragon Plus Environment

Page 6 of 43

Page 7 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

(PME)41 was employed to calculate the long-range electrostatic interactions. The cutoff distance for the long range van der waals (VDW) energy term was 12.0Å. In order to remove any steric conflicts induced during system setup, structural optimizations were first performed on the relaxed water molecules and counter ions in two steps with the harmonic constraint potential of 2.0 kcal·(mol·Å2)−1 on all heavy atoms of both protein and ligands. Afterward, the whole system was minimized without any restraint. The above steps were all executed by 2500 cycles of steepest descent minimization followed by 5000 cycles of conjugate gradient minimization. After system optimization, running MD simulations was started on the systems by gradually heating each system in the NVT ensemble from 0 to 300 K in 50 ps using a Langevin thermostat with a coupling coefficient of 1.0/ps and with a force constant 2.0 kcal·(mol·Å2)−1 on the complex. And then 500 ps of density equilibration with a force constant 2.0 kcal·(mol·Å2)−1 on the complex were performed. Subsequently the systems were again equilibrated for 500ps by releasing all the restrains. Finally, production runs for 20 ns MD simulations were performed under the constant temperature of 300 K in the NPT ensemble with periodic boundary conditions for each systems During MD procedure, the SHAKE algorithm42 was applied for the constraint of all covalent bonds involving hydrogen atoms. The time step was set to 2 fs. 2.3. Principal Component Analysis (PCA) Principal components analysis (PCA) is an effective tool to filter global slow motions from local fast motions and extract the slow and functional motions from the

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

MD trajectories by dimensional reduction method.28,32 It elucidates the essential motion governing conformation transition during the simulation. This analysis based on the calculation and diagonalization of the convariance matrix of Cα atomic fluctuation to obtain orthogonal eigenvectors (also called principal components, PCs, give the direction of the principal mode) and the corresponding eigenvalues (represent the magnitude of the motions along direction).31, 43, 44 Usually, the first few PCs with higher eigenvalues describe the slow motion modes of the systems, which contain the functional motion and display the most important conformation changes of a bio-molecular system. In our work, the PCA analysis was performed with Gromacs4.5 package29 to obtain the significant motions during the simulation in the four systems.45 2.4 Binding Free Energy Calculations and Per-Residue Free Energy Decomposition Analysis. The binding free energy was calculated using the MM/GBSA procedure implemented in AMBER 12.39 After the stability of each system were verified by RMSD for the backbone atoms, the average 1000 snapshots were extracted from the last 4ns stable MD trajectory at 4 ps intervals. For each snapshot, binding free energy was calculated as the difference between the free energy of the complex (G complex) and the total of the free energies of the receptor (G receptor) and the ligand (G ligand), shown in the following equation:

ACS Paragon Plus Environment

Page 8 of 43

Page 9 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

∆G

(binding)

∆G(binding

=G

(complex)

− (G

(receptor)

+G

(ligand)

)

) = ∆G(MM) + ∆G(sol) − TS

∆G(MM) = ∆G(ELE) + ∆G(VDW) ∆G(sol) = ∆Gele,sol ∆Gnp,sol

+ ∆Gnp,sol

= γSASA + β

where ∆G(MM), ∆G(sol), TS denote the total molecular mechanics energy(gas phase, including the non-bonded electrostatic energy ∆G(ELE)and van der Waals interaction ∆G(VDW)), the solvation free energy ∆G(sol) and the entropy term(-TS), respectively. The solvation free energy(∆G(sol)) is composed of two components: the polar contribution to solvation(∆Gele,sol) and the nonpolar contributions to solvation(∆Gnp,sol). The former was obtained by solving the generalized Born (GB) model by means of MM/GBSA methods.46 The parameters used for the GB calculation were developed by Onufriev et al (GBOBC, igb = 2). The latter was estimated as a function of the solvent accessible surface area (SASA), which calculated using LCPO (linear combination of pairwise overlaps) method with a solvent-probe radius of 1.4 Å using the amber Molsurf module., where γ is the surface tension (0.0072 kcal /mol Å2) and β is constant 0.47, 48 The analyses of entropy contributions were carried out using normal-mode analysis with the NMODE program in AMBER12. Due to the high computational demand of this approach, only 20 snapshots evenly extracted from the last 2 ns MD trajectory were extracted to calculate the entropic contribution. To obtain a detailed view of the protein-ligand binding and identify the key residues responsible for the binding, free energy decomposition to each residue was performed using the MM/GBSA method. The decomposition was performed only for molecular

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

mechanics and solvation energies but not for entropies. The binding energies were decomposed into contributions from protein-ligand interaction pairs, which consists of four terms: ∆G(ELE), ∆G(VDW), ∆Gele,sol, ∆Gnp,sol. All energy decomposition analyses were performed on the same snapshots which were used in above calculations.

3. RESULTS AND DISCUSSIONS: 3.1 RMSD and RMSF to check the systems equilibration and flexibility: The dynamic stability of complexes should be explored firstly when 20 ns MD simulations have been successfully finished for four systems in explicit aqueous solution. The RMSD is the root-mean-square deviation of the simulated positions of the backbone atoms of proteins from those in the initial X-ray crystal structures. The RMSDs for backbone atoms of proteins from four systems relative to their starting structures were analyzed for four systems, as shown in Figure 2.

Figure 2. Plots of RMSD for all the protein backbone atoms (Å) relative to their initial conformation vs simulation time (ns) for four systems.

ACS Paragon Plus Environment

Page 10 of 43

Page 11 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

As can be seen in the plots, after 7.5 ns simulation, only small RMSD fluctuations indicate the four systems are stable and equilibratory. As a whole, more fluctuations are displayed in ROS1 systems than in ALK systems, which is consistent with the resolution of their X-ray crystal structures. For each target, the selective system has comparatively little larger fluctuations than the dual system. The stability of the systems was also verified by analysis of the hydrogen bonds between ligands and receptors.

Figure 3. (A)The labeled key loops of ROS1:7c complex. Hinge region is colored by blue; DFG motif is colored by orange; P-loop is colored by red and Alpha-loop is colored by magenta. (B) The diagram of key H-bonds of ROS1:7c complex during the MD simulation. The flexibility of each residue is assessed by its root-mean square fluctuation (RMSF). The residues with higher fluctuation values indicate they are more flexible during the molecular dynamic. In order to depict the difference of the fluctuation of

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

key loops in detail, the important loops to the binding were labeled in Figure 3(A). The RMSFs of main chain versus the residue number for complexes were also calculated for each MD-simulated complex, illustrated in Figure 4.

Figure 4. RMSF of each residue of the target ALK(left) and ROS1(right) systems obtained from 20 ns MD simulation. As shown in Figure 4, the hinge regions and DFG-motifs of two targets have similar small fluctuation patterns, whereas the P-loop and alpha-loop exhibits similar larger fluctuations, indicating the P-loop and alpha-loop are more flexible than hinge regions and DFG-motif. In addition, the P-loops and alpha-loops in dual inhibitor systems (ALK:crizotinib and ROS1:crizotinib), have increasing fluctuations than those in selective inhibitor systems (ROS1:7c and ALK:alectinib), indicating the binding of selective inhibitors leads to more stable than binding to the dual inhibitors, especially in ROS1:7c complex. The decreased flexibility in the certain areas can be explained by the more ligand-receptor interactions between 7c and ROS1. Moreover, the compound 7c, which has more rotatable bonds in the binding site, is relatively larger than the dual inhibitor crizotinib with bulkier group. Therefore, it needs much more space for the binding. The additional binding site that probably created by

ACS Paragon Plus Environment

Page 12 of 43

Page 13 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

P-loop, Cα-helix and DFG-motif would be formed mainly by the movement of P-loop and alpha-loop. Then its conformations would change to match each other and achieve the optimal binding mode. Above all, it was presumed that the interactions between ligand and P-loop and alpha loop would be contributed to the selectivity of ROS1.

3.2. Hydrogen bonds analysis to find the key residues The crystallographic and former docking study indicated that the hydrogen bonds mediated the important interactions between protein and inhibitor. The hydrogen bonds from the equilibrium trajectories of the four systems were calculated by VMD 1.9 with a distance cutoff value of 3.5 Å and an angle cutoff value of 35°. The results are listed in Tables 1 with occupancy over 10%. For the dual inhibitor crizotinib, the mono- or double H-bonds with high occupancy in hinge region of two targets displayed the important interaction (for ALK:Glu105 and Met107; for ROS1:Glu94 and Met96). However, for selective inhibitor, the H-bonds from hinge region decreased. The P-loop involved in the interaction between protein and ligand. For instance, in the selective ROS1:7c systems(shown in Figure 2(B)), more hydrogen bonds with shorter distances from P-loop (Gly24 and Gly21 ) were formed with ligand, while the hydrogen bonds with Glu94 or Met96 were unstable with low occupancy 0.22% and larger distances. Furthermore, the residue Lys47 from β3 sheet and DFG-motif also displayed a long life H-bond. To quantitatively analyze the hydrogen bonds, we compared the lengths of H-bond during the simulation time (see supplementary Figure S1). It could be seen that the H-bonds made by Gly24 with the

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 43

phenolic hydroxyl group and Asp169 with CN group were shorter than those of formed by Met96 and Lys47 during the last 2ns. The P-loop and DFG motif served as a pincer to hold the ligand 7c and the hydrophobic interaction from P-loop and Leu153 aided to form favorable interaction. However, in selective ALK:alectinib system, the H-bonds formed by hinge regions played a key roles with high occupancy, and the H-bonds from Lys65 and Arg124 with low occupancy could stabilize the conformation of ligand in binding sites. The above analysis displayed the P-loop and DFG-motif would be contributed to the selective to ROS1. Therefore, structural optimization to enhance the hydrogen bond in P-loop and DFG-motif may be a reasonable solution to improve the selectivity and activity to ROS1. Table 1. The hydrogen bonds analysis of the four systems from the MD trajectories ALK:alectinib ALK:crizotinib donor

acceptor

occupancy

donor

MET114-Main-N

MOL317-Side-O1

99.61%

MOL310-Side-N5

GLU105-Main-O

99.40%

ARG124-Side-NH2

MOL317-Side-O2

13.93%

MET107-Main-N

MOL310-Side-N4

98.74%

LYS65-Side-NZ

MOL317-Side-N2

11.92%

ASP178-Main-N

MOL310-Side-F1

42.18%

GLY31-Main-N

MOL310-Side-N3

12.57%

ROS1:7c

acceptor

occupancy

ROS1:crizotinib

MOL298-Side-O1

GLY24-Main-O

74.20%

MOL298-Side-N2

GLU94-Main-O

99.48%

LYS47-Side-NZ

MOL298-Side-N3

72.80%

MET96-Main-N

MOL298-Side-N1

97.40%

PHE170-Main-N

MOL298-Side-N2

68.20%

ASP169-Main-N

MOL298-Side-F1

54.09%

MOL298-Side-O2

ASP169-Side-OD1

31.20%

MOL298-Side-N5

ASP100-Side-OD1

20.21%

MOL298-Side-O2

ASP169-Side-OD2

29.40%

MOL298-Side-N5

LEU18-Main-O

15.02%

MOL298-Side-O1

GLY21-Main-O

13.00%

ACS Paragon Plus Environment

Page 15 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

3.3. Structure analysis of dynamic results to understand the selectivity mechanism The X-ray cocrystal/docking and MD results of four complexes revealed that three inhibitors shared a similar typically DFG-in binding mode. Figure 5 displayed the structure and conformation comparison of representative snapshots from MD of four systems. As seen in Figure 5, the hinge regions of four systems were almost the same. Obviously, it was not the specific loop responsibility for the selectivity. For ALK target (Figure 5A), the selective inhibitor alectinib and the dual inhibitor crizotinib shared good similar conformation binding to the receptor, except that the former have larger group in solvent region near the hinge regions. It is noteworthy that the P-loop of selective system(ALK:alectinib, yellow) extended down to the active site, which was distinct from that of dual inhibitor system(ALK:crz, cyan). The alpha-loop seems have no much difference. For ROS1 target (Figure 5B), the P-loop has similar conformation, while the alpha-loop displays obvious difference. The alpha-loop of selective system (ROS1:7c, blue) appears behind that of dual system(ROS1:crz, orange), with the results of amplified the active site. For dual inhibitor systems (Figure 5C), the conformations of P-loop and alpha-loop have minute difference although the ligand is same, except that the β2-sheet linked by P-loop of ROS1 extends up to active site. For ALK/ROS1 selective inhibitor systems (Figure 5D), the distinct difference of P-loop and alpha-loop could be seen apparently. The P-loop of selective ROS1 system lies upper and behind than that of the selective ALK system (7.1Å between the Ala41:ALK and Ala22:ROS1). Additional, the alpha-loop of selective ROS1 system lies behind that of selective ALK system (3.5Å between

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Cys71 and Ser53). The divergence leads to the different conformation, which probably responsible for the selectivity of ROS1 versus ALK. As a result, the two loops should get our attention.

Figure 5. Structure comparison between different system snapshot from MD of: (A) selective(yellow) and dual(cyan) ALK complexes, (B)selective(blue) and dual(orange) ROS1 complexes, (C)dual ALK(cyan) and ROS1(pink) system, (D)selective ALK(yellow) and ROS1(blue) complexes. The middle small pictures are the conformation comparison of ligands. In order to show the main differences of selective ALK/ROS1 systems, the residues except the key loops were colored by gray.

3.4. PCA analysis to disclose the P-loop motion mode In order to further investigate the functional motions of the four systems, the PCA analysis was performed for Cα atoms of protein based on the equilibrium trajectories.

ACS Paragon Plus Environment

Page 16 of 43

Page 17 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 6 compares the first and second slowest motion modes and motion correlation of the selective inhibitor systems. The orientation of the cone describes motive direction, and the length of the cone is positively correlated with motive magnitude. The heat map presented the correlated fluctuations of selective system graphically. Positive correlations were mapped in the upper left triangle, while negative correlations were mapped in the lower right triangle. Deeper color indicated more correlated (or anti-correlated). The P-loop and alpha-loop domain was highlighted by black cycle. As shown in Figure 6, it was obvious that the first and second slow modes of the P-loop and alpha-loop in two systems were quite different. For ALK-selective system, the P-loop moved to the left of ATP-binding site, and displayed the small fluctuation. The alpha loop moved to the front of ATP-binding sites, and alpha helix also had small fluctuation. It is in agreement with the RMSF analysis. Comparing the dual inhibitor of ALK system, the P-loop of ALK has twisted movements to keep it relatively stable position (shown in supporting information Figure S2). For ROS1-selective system, the P-loop move as up-down mode, shown as arrow. When the P-loop moved down to the ATP binding sites, the alpha-loop moved to the front. Both of them showed the small motive magnitude. The second motion mode showed that when the P-loop move upper with the alpha-loop move to the rear-left by large magnitude. With the movement of alpha-loop, the β3-sheet would rise and then appear another active site composed with the alpha-helix which can interact with ligands. At the same time, the alpha loop and alpha-helix combining P-loop can

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

maintain the new binding sites. In the dual-ROS1 system, P-loop and alpha-loop have larger movement magnitude and similar movement direction comparing to that of in selective system (shown in Figure S2 dual-ROS1 systems, PC2). Moreover, in the dual ALK/ROS1 systems, we can also see the distinct movements’ differences of P-loop and alpha-loop.

Figure 6. The first and second slowest motion modes(PC1 and PC2) of the ALK-selective and ROS1-selective complexes and heat map. The length of cone is

ACS Paragon Plus Environment

Page 18 of 43

Page 19 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

correlated with motive magnitude positively, and the orientation of cone indicates motive direction. The P-loop was colored by green and the alpha-loop was colored by cyan. The heat map illustrates the motion correlation. Deeper color indicated more correlated (or anti-correlated). The P-loop and alpha-loop domain are highlighted by black cycle. The heat map of motion correlation showed that the global dynamics of the systems were very similar. The ROS1-selective system has comparatively more positive correlation than ALK-selective systems since it has relatively bigger deeper color area in the upper-left triangle. Additionally, the P-loop and alpha-loop domain (highlighted by black cycle) in the ROS1-selective systems displayed more deeper color and bigger area, indicating the motion of the domains in ROS1-selective has more positive correlation. Taken together the analysis of motion and motion correlation above, it was found that the motion directions of P-loop and the alpha loop in ROS1 and ALK were quite different. In selective ROS1 inhibitor systems, the P-loop and alpha-loop played a key role to selectivity. However, in selective ALK inhibitor systems, it seemed not so important. In dual ROS1/ALK inhibitor system, the P-loop and alpha loop also displayed quite different, indicating that the interactions between targets and ligands may be dissimilar. Therefore, we speculated that the P-loop and alpha-loop had more flexible to interact with the ligands. When they interacted with the suitable volume ligands, P-loop moved to the upper and alpha-loop moved to the front of the active binding site, and then informed the other selective active site combined with the

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 43

alpha-helix and DFG-motif. The new active site may be responsible for the selectivity of type-I ROS1 inhibitors.

3.5. MM/GBSA binding energy and decomposition analysis Binding free energy was considered to reflect the binding affinity between the ligand and receptor. To evaluate the binding energy of four complexes, MM/GBSA methodology was applied to calculate binding free energies in terms of gas phase energy, solvation energy, and entropic contributions. Since the binding free energy obtained by Amber MM/GBSA is not the absolute total energy because of the calculation error of the entropy contribution, the computed values do not reproduce the absolute experimental values accurately. However, they have still been shown to correlate with the experiment values well to the same target.26, 49 Table 2 Calculated binding free energies of four systems using MM/GBSA method (kcal/mol). Energy 3aox 2xp2 3zbf ROS1:7c

ΔGVDW ΔGELE ∆Ggas ∆Gele,sol(GB) ∆Gnonpol,sol Esurf ∆Gsol ∆Ggas+∆Gsol T∆S ∆G pred(GB) K(i)(nmol/L)

-51.23(±0.10) -8.42(±0.12) -59.65(±0.15) 26.57(±0.11) -5.66(±0.01) 20.91(±0.10) -38.73(±0.09) -15.28(±2.50) -13.45 0.83

-44.04(±0.12) -13.76(±0.13) -57.80(±0.18) 27.59(±2.36) -4.97(±0.01) 22.62(±0.10) -35.18(±0.12) -17.49(±1.95) -17.69 0.74-8.2

-46.89(±0.12) -22.80(±0.15) -69.69(±0.18) 37.77(±0.12) -5.48(±0.01) 32.29(±0.11) -37.40(±0.12) -17.64(±1.80) -19.76 1.7*

-53.82(±0.15) -39.92(±0.44) -93.74(±0.40) 62.87(±0.33) -7.19(±0.01) 55.68(±0.33) -38.06(±0.15) -23.24(±1.80) -14.84 24*

* the value is IC50 Table 2 summarized computed binding free energies along with their respective enthalpic and entropic contributions for four systems and experimental Ki/IC50 values.

ACS Paragon Plus Environment

Page 21 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Overall, our results showed that the general trends of computed binding free energy values were in accordance with their experimental activities for each target. Especially for ROS1 target, the ∆G(bind) values of the selective inhibitor 7c complex and the dual inhibitor crizotinib complex are -14.84kcal/mol and -19.76 kcal/mol, responding to their IC50 are 24 and 1.7nmol/L, individually. Therefore, the MM/GBSA binding free energies are expected to reflect their relative in vivo activities. One of the advantages of the MM/GBSA method is that it enables a decomposition of the free energy into identifiable contributions. The free energy components responsible for the binding were further explored separately to get insights into driving forces for selective bindings of kinases. As seen in Table 2, the intermolecular van der Waals(∆GVDW) and the electrostatic (∆GELE) interaction provided the driving force for binding for the four system. However, the total solvation energy (∆Gsol) consisting of polar (∆Gele,sol(GB)) and nonpolar(∆Gnonpol,sol) terms were unfavorable for all the complexes. In contrast to the nonpolar solvation, the polar solvation interactions were found to make unfavorable contribution to the binding energy with the positive ∆Gele,sol(GB) values. Although there were some hydrogen bonds between the ligands and receptors, their contributions could not compensate the large desolvation penalties during the interaction, thereby always leading to an unfavorable contribution. This compensation phenomenon had been previously observed in several studies of protein-ligand interactions in solution.26, 50 Apparently, the van der Waals interactions dominated the inhibitor binding and determine the binding

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

specificity of the inhibitors as well. Comparing the two systems of dual inhibitor, they had similar van der Waals energies and different electrostatic energies. Considering the different desolvation, they had the similar binding energies and showed the similar activities. For target ALK, the selective-ALK complex had larger VDW contribution and smaller ELE contribution than the dual-ALK complex, indicating that the selective-ALK complex had larger hydrophobic surface near the hinge region, which was consistent with the nonpolar solvation energy proportional to the solvent accessible surface areas. Similarly, the VDW of selective-ROS1 complex was also more negative than that of the dual-ROS1 complex. It was surprising that the both the ∆GELE and ∆GVDW also showed the more significant negative in selective-ROS1 versus dual-ROS1 systems while the selective ROS1 inhibitor 7c displayed the lower activity. However, decomposition of the solvation energy contribution into polar and nonploar solvation components showed that the favorable contribution from hydrophobic interaction were far less than the unfavorable polar solvation contribution. Furthermore, it was worth noting that the difference of entropy since the conformational entropy determined the conformational change upon drug binding. Comparing the dual-ROS1 system, the entropy of the selective-ROS1 system has more negative, leading that the total binding free energy of selective-ROS1 less than that of the dual-ROS1 systems. Analysis of free energy components indicated the VDW interaction and larger solvent accessible surface area may contribute to the selectivity between 7c and ROS1. Combining the structure analysis and hydrogen bond analysis, we concluded

ACS Paragon Plus Environment

Page 22 of 43

Page 23 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

that the N-carbonitrile-3-(3-hydroxy-5-methyl)phenyl-1H-pyrazole group may play a critical role in the selectivity and activity of 7c, which was agreement with the estimation from So Ha Lee’ team. To further elaborate the key residues related to the binding selectivity in-depth, the calculated binding free energies were decomposed into the contribution of single-residue using the MM/GBSA method24, 47, 48, 51 with the default parameters embedded in Amber 12. Figure 7 displayed the most important residues contributing to the binding free energies for the two selective complexes. Generally, the residue would be important for the binding behavior when the decomposed energy of a residue was more negative than -1.0 kcal/mol. Undoubtedly, Asp169 from DFG motif, Val26 from β2 and Leu47 from β3 of selective ROS1 system would be the predominant contribution to the binding energy. The decreased interactions with hinge region weakened the activity to ROS1. However, in the selective ALK system, these corresponding interactions are attenuated due to downward behavior of the P-loop regions, and the contribution of Leu37 (corresponding to Leu18 in ROS1) from β1 increased. The interaction from β1-sheet and hinge region may responsible for the ALK-selective system, while the β2-sheets linked by P-loop and β3 sheet linked by alpha-loop combining the DFG motif may responsible for the ROS1-selective systems. The movement of P-loop and alpha-loop lead to the conformational changes of β2 and β3-sheets, then resulted in the different contribution of β-sheets in selective systems. Their flexibilities were in accord with the analysis of RMSF and movement motion analysis. The P-loop was considered

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

closely related to the interaction specificity reported in the recently study.47 Therefore, the P-loop, alpha-loop and the DFG motif are considered to be responsible to the selectivity of type-I selective ROS1 system.

Figure 7. Per-residue binding free energy decomposition of selective systems to find the key residues. ROS1:7c complex (top) and ALK: Alectinib complex (bottom). The total binding energy plots are shown on the left and the corresponding contributions of val der Waals and Electrostatics are shown on the right.

3.6. Experimental Validation: In order to validate our conclusion, some compounds were designed and synthesized based on the results. Among them, the compound 14c and SMU-B (structure data were displayed in supporting information) were chosen as examples to illustrate the design details. Their structures and inhibitory activities against ALK and ROS1 enzymes were displayed in Figure 8. In Figure 8, the blue portions could form

ACS Paragon Plus Environment

Page 24 of 43

Page 25 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

mono- or double H-bonds with the hinge regions, the red portions laid below the P-loop and the green portions point at the solvent near the hinge regions. According our modeling results, the active site composed by P-loop and DFG-motif would benefit to the selectivity of ROS1 versus ALK, while the modification of solvent regions would maintain the dual-inhibitory. In order to improve the selectivity to ROS1 versus ALK, the 2,6-dichloro-3-fluorobenzyloxyl (red group in crizotinib) was replaced by N-(3,4-dimethoxybenzyl)-3-methyl-benzamide(red group in 14c) and we obtained new compound 14c. It showed a relatively selectivity to ROS1 with the IC50 0.365µM/4.45µM against ROS/ALK. The compound SMU-B was obtained by modification of the solvent regions. It showed good activity to both ROS1 and ALK at nM level as we expected (kinases inhibitory activity experimental details were described in supporting information).

Figure 8. The new compounds designed based on our strategy Although the activity of 14c was not good enough, the selectivity still had been improved. Docking and MD simulation of compound 14c showed that the binding

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 43

mode of 14c was similar to that of crizotinib except the selective active site (shown in Figure 9A). The red group interacted with P-loop by hydrophobic and formed H-bond with Asp169 from DFG-motif. The first reason of unsatisfactory activity of ROS1 probably was in the steric hindrance from the ethoxyl group. During the MD simulation, the conformation of methyl turned almost 180 degree from inner to outer of the active site (shown in Figure 9B, colored by magenta). It directly led to the decrease of interaction between the hinges region/the DFG-motif and ligand. The Lys47 has not involved in the interaction and decrease the electrostatic interaction may be the other reason. The further modification for the compound 14c is underway. Comparing the 14c, SMU-B displayed potent dual ROS1/ALK activity by modification of solvent regions (green). Docking results also showed that they adopted the

similar DFG-in

conformation

to the

crizotinib.

The

group

2,6-dichloro-3-fluorobenzyloxyl did not take up the selective pocket composed by P-loop and DFG-motif and displayed no selectivity(shown in Figure 9C).

Figure 9. The binding modes of compound 14c and SMU-B. (A) The binding mode and H-bonds of 14c at 4ns of MD simulation. (B) The conformation change of 14c during the MD simulation. The ethoxyl group (magenta) turned from the inner (yellow) to the outer(green) of the active sites. (C) The conformation comparison of 14c(yellow) and SMU-B(cyan).

ACS Paragon Plus Environment

Page 27 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

4. CONCLUSION In order to investigate the binding selectivity of type-I inhibitors to target ROS1/ALK, a computational strategy combining MD simulation, movement motion analysis, MM/GBSA free energy calculation were applied to clarify the binding selectivity by comparing the selective and dual ROS1/ALK inhibitor systems. Our results unambiguously suggested that the binding selectivity was controlled by the differences of P-loop and alpha-loop, which led to the conformational changes of β-sheets. The interaction from DFG motif was beneficial to form selective binding pocket of ROS1. Besides, the energy decomposition analysis also demonstrated that the P-loop and DFG-motif were responsible for the binding specificity for ROS1, while the interactions hinge-region and solvent regions nearby were important for the dual inhibitor systems. The new designed compound 14c and SMU-B including their kinases inhibitory activities verified our modeling results and design strategy. These observations provided a better structural understanding of mechanism of type-I ROS1/ALK inhibitors. They were also expected to help guide the discovery and rational design of highly selective and potent ROS1 inhibitors. ASSOCIATED CONTENT Support information Figure S1. The distances of key H-bonds plots of ROS1:7c system during last 2ns of molecular dynamics.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure S2. The first and second slowest motion modes (PC1 and PC2) of the dual-ALK and dual-ROS1 systems. The structure data of compound 14c and SMU-B. The kinases inhibitory activities experiment protocol. Figure S3. The IC50 curve of compound 14c against ROS1 kinse. AUTHOR INFORMATION Corresponding Author * Prof. Jiajie Zhang, Tel: 86-20-61648548. E-mail:[email protected]

Present Address Current address for Prof. Zhang: School of Pharmaceutical Sciences, Southern Medical University, Guangzhou, 510515, P.R China. Notes The authors declare no competing financial interest. ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No. 81573263 and 31600591), Natural Science Foundation of Guangdong Province, China (No.2015A030313285) and Science and Technology Planning Project of Guangdong Province, China (No.2014A020210012 and 2016A020210087). We also acknowledge Professor David Case for the kind gift of AMBER 12 software. The calculations of the ligands were performed in the High Performance Cluster of Prof. Zhenhai Zhang’s group from Nanfang Hospital.

ACS Paragon Plus Environment

Page 28 of 43

Page 29 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

REFERENCES (1). Minuti, G.; D'Incecco, A.; Cappuzzo, F., Targeted Therapy for NSCLC with Driver Mutations. Expert. Opin. Biol. Ther. 2013, 13, 1401-1412. (2). Choi, C. M., Overview of ALK and ROS1 Rearranged Lung Cancer. Tuberc. Respir. Dis. (Seoul). 2013, 75, 236-237. (3). Wheeler, D. L.; Yarden, Y., Receptor Tyrosine Kinases : Family and Subfamilies. Humana Press: Cham, 2015. (4). Gainor, J. F.; Shaw, A. T., Novel Targets in Non-Small Cell Lung Cancer: ROS1 and RET fusions. Oncologist. 2013, 18, 865-875. (5). Ou, S. H.; Tan, J.; Yen, Y.; Soo, R. A., ROS1 as a 'Druggable' Receptor Tyrosine Kinase: Lessons Learned from Inhibiting the ALK Pathway. Expert. Rev. Anticancer. Ther. 2012, 12, 447-456. (6). Bergethon, K.; Shaw, A. T.; Ou, S. H.; Katayama, R.; Lovly, C. M.; McDonald, N. T.; Massion, P. P.; Siwak-Tapp, C.; Gonzalez, A.; Fang, R.; Mark, E. J.; Batten, J. M.; Chen, H.; Wilner, K. D.; Kwak, E. L.; Clark, J. W.; Carbone, D. P.; Ji, H.; Engelman, J. A.; Mino-Kenudson, M.; Pao, W.; Iafrate, A. J., ROS1 Rearrangements Define a Unique Molecular Class of Lung Cancers. J. Clin. Oncol. 2012, 30, 863-870. (7).de la Bellacasa R P; Karachaliou N; Estrada-Tejedor R, ALK and ROS1 as a Joint Target for the Treatment of Lung Cancer: a review. Transl. Lung. Cancer. Res. 2013, 2, 72-86. (8). Aisner, D. L.; Nguyen, T. T.; Paskulin, D. D.; Le, A. T.; Haney, J.; Schulte, N.; Chionh, F.; Hardingham, J.; Mariadason, J.; Tebbutt, N.; Doebele, R. C.; Weickhardt, A. J.; Varella-Garcia, M., ROS1 and ALK Fusions in Colorectal Cancer, with Evidence of Intratumoral Heterogeneity for Molecular Drivers. Mol. Cancer. Res. : MCR 2014, 12, 111-118. (9). Bos, M.; Gardizi, M.; Schildhaus, H. U.; Buettner, R.; Wolf, J., Activated RET and ROS: Two New Driver Mutations in Lung Adenocarcinoma. Transl. Lung. Cancer. Res. 2013, 2, 112-121. (10). Song, A.; Kim, T. M.; Kim, D. W.; Kim, S.; Keam, B.; Lee, S. H.; Heo, D. S., Molecular Changes Associated with Acquired Resistance to Crizotinib in ROS1-Rearranged Non-Small Cell Lung Cancer. Clin. Cancer. Res. 2015, 21, 2379-2387. (11). Zou, H. Y.; Li, Q.; Engstrom, L. D.; West, M.; Appleman, V.; Wong, K. A.; McTigue, M.; Deng, Y. L.; Liu, W.; Brooun, A.; Timofeevski, S.; McDonnell, S. R.; Jiang, P.; Falk, M. D.; Lappin, P. B.; Affolter, T.; Nichols, T.; Hu, W.; Lam, J.; Johnson, T. W.; Smeal, T.; Charest, A.; Fantin, V. R., PF-06463922 is a Potent and Selective Next-generation ROS1/ALK Inhibitor Capable of Blocking Crizotinib-resistant ROS1 Mutations. Proc. Natl . Acad. Sci . USA 2015, 112, 3493-3498. (12). Johnson, T. W.; Richardson, P. F.; Bailey, S.; Brooun, A.; Burke, B. J.; Collins, M. R.; Cui, J. J.; Deal, J. G.; Deng, Y. L.; Dinh, D.; Engstrom, L. D.; He, M.; Hoffman, J.; Hoffman, R. L.; Huang, Q.; Kania, R. S.; Kath, J. C.; Lam, H.; Lam, J. L.; Le, P. T.; Lingardo, L.; Liu, W.; McTigue, M.; Palmer, C. L.; Sach, N. W.; Smeal, T.; Smith, G. L.; Stewart, A. E.; Timofeevski, S.; Zhu, H.; Zhu, J.; Zou, H. Y.; Edwards, M. P., Discovery

of

(10R)-7-amino-12-fluoro-2,10,16-trimethyl-15-oxo-10,15,16,17-tetrahydro-2H-8,4-(m

etheno)pyrazolo[4,3-h][2,5,11]-benzoxadiazacyclotetradecine-3-carbonitrile

(PF-06463922),

a

Macrocyclic Inhibitor of Anaplastic Lymphoma Kinase (ALK) and c-ros Oncogene 1 (ROS1) with Preclinical Brain Exposure and Broad-spectrum Potency against ALK-resistant Mutations. J. Med. Chem. 2014, 57, 4720-4744. (13). Heigener, D. F.; Reck, M., Crizotinib. Recent Results in Cancer Research. Fortschritte der Krebsforschung. Progres dans les recherches sur le cancer 2014, 201, 197-205.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14). Frampton, J. E., Crizotinib: a Review of Its Use in the Treatment of Anaplastic Lymphoma Kinase-positive, Advanced Non-Small Cell Lung Cancer. Drugs 2013, 73, 2031-2051. (15). Shaw, A. T.; Solomon, B. J., Crizotinib in ROS1-rearranged Non-Small-Cell Lung Cancer. New Engl. J. Med. 2015, 372, 683-684. (16). Shaw, A. T.; Ou, S. H.; Bang, Y. J.; Camidge, D. R.; Solomon, B. J.; Salgia, R.; Riely, G. J.; Varella-Garcia, M.; Shapiro, G. I.; Costa, D. B.; Doebele, R. C.; Le, L. P.; Zheng, Z.; Tan, W.; Stephenson, P.; Shreeve, S. M.; Tye, L. M.; Christensen, J. G.; Wilner, K. D.; Clark, J. W.; Iafrate, A. J., Crizotinib in ROS1-rearranged Non-Small-Cell Lung Cancer. New. Engl. Med. 2014, 371, 1963-1971. (17).Lee, J.; Lee, S. E.; Kang, S. Y.; Do, I. G.; Lee, S.; Ha, S. Y.; Cho, J.; Kang, W. K.; Jang, J.; Ou, S. H.; Kim, K. M., Identification of ROS1 Rrearrangement in Gastric Adenocarcinoma. Cancer 2013, 119, 1627-35. (18). Davare, M. A.; Saborowski, A.; Eide, C. A.; Tognon, C.; Smith, R. L.; Elferich, J.; Agarwal, A.; Tyner, J. W.; Shinde, U. P.; Lowe, S. W.; Druker, B. J., Foretinib is a Potent Inhibitor of Oncogenic ROS1 Fusion Proteins. Proc. Natl . Acad. Sci . USA 2013, 110, 19519-19524. (19). Chin, L. P.; Soo, R. A.; Soong, R.; Ou, S. H., Targeting ROS1 with Anaplastic Lymphoma Kinase Inhibitors: a Promising Therapeutic Strategy for a Newly Defined Molecular Subset of Non-Small-Cell Lung Cancer. J. Thorac. Oncol. 2012, 7, 1625-1630. (20). M.M. Al-Sanea; A.Z. Abdelazem; B.S. Park; K.H. Yoo; T. Sim; Kwon, Y. J.; Lee, S. H., ROS1 Kinase Inhibitors for Molecular-Targeted Therapies, Curr. Med. Chem. 2016, 23, 142-160. (21).Abdelazem, A. Z.; Al-Sanea, M. M.; Park, B. S.; Park, H. M.; Yoo, K. H.; Sim, T.; Park, J. B.; Lee, S. H.; Lee, S. H., Synthesis and Biological Evaluation of New Pyrazol-4-ylpyrimidine Derivatives as Potential ROS1 Kinase Inhibitors. Eur. J. Med. Chem. 2015, 90, 195-208. (22). Park, B. S.; Al-Sanea, M. M.; Abdelazem, A. Z.; Park, H. M.; Roh, E. J.; Park, H. M.; Yoo, K. H.; Sim, T.; Tae, J. S.; Lee, S. H., Structure-Based Optimization and Biological Evaluation of Trisubstituted Pyrazole as a Core Structure of Potent ROS1 Kinase Inhibitors. Bioorg Med Chem 2014, 22, 3871-3878. (23). Davare, M. A.; Vellore, N. A.; Wagner, J. P.; Eide, C. A.; Goodman, J. R.; Drilon, A.; Deininger, M. W.; O'Hare, T.; Druker, B. J., Structural Insight into Selectivity and Resistance Profiles of ROS1 Tyrosine Kinase Inhibitors. Proc. Natl . Acad. Sci . USA 2015, 112, E5381-E5390. (24). Ravichandran, S.; Luke, B. T.; Collins, J. R., Can Structural Features of Kinase Receptors Provide Clues on Selectivity and Inhibition? A Molecular Modeling Study. J Mol Graph Model 2015, 57, 36-48. (25). Chen, S. F.; Cao, Y.; Chen, J. J.; Chen, J. Z., Binding Selectivity Studies of PKBα using Molecular Dynamics Simulation and Free Energy Calculations. J Mol Model 2013, 19, 5097-5112. (26). Liu H. X., Yao X. J. , Molecular Basis of the Interaction for an Essential Subunit PA-PB1 in Influenza Virus RNA Polymerase: Insights from Molecular Dynamics Simulation and Free Energy Calculation. Mol. Pharm 2010, 1, 75-85. (27). Zelenko, U.; Hodoscek, M.; Rozman, D.; Golic Grdadolnik, S., Structural Insight into the Unique Binding Properties of Pyridylethanol(phenylethyl)amine Inhibitor in Human CYP51. J. Chem. Inf. Model. 2014, 54, 3384-3395. (28). Tripathi, S.; Srivastava, G.; Sharma, A., Molecular Dynamics Simulation and Free Energy Landscape Methods in Probing L215H, L217R and L225M BetaI-tubulin Mutations Causing Paclitaxel Resistance in Cancer Cells. Biochem. Bioph. Res. Co 2016, 476, 273-279. (29). Wan, H.; Chang, S.; Hu, J. P.; Tian, Y. X.; Tian, X. H., Molecular Dynamics Simulations of Ternary Complexes: Comparisons of LEAFY Protein Binding to Different DNA Motifs. J. Chem. Inf. Model. 2015, 55, 784-794.

ACS Paragon Plus Environment

Page 30 of 43

Page 31 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

(30). Wan, H.; Hu, J. P.; Tian, X. H.; Chang, S., Molecular Dynamics Simulations of Wild Type and Mutants of Human Complement Receptor 2 Complexed with C3d. Phys. Chem. Chem. Phys. 2013, 15, 1241-1251. (31). Wang, F.; Wan, H.; Hu, J. P.; Chang, S., Molecular Dynamics Simulations of Wild Type and Mutants of Botulinum Neurotoxin A Complexed with Synaptic Vesicle Protein 2C. Mol. BioSys. 2015, 11, 223-231. (32). Zhou, Y.; Zhang, N.; Chen, W.; Zhao, L.; Zhong, R., Underlying Mechanisms of Cyclic Peptide Inhibitors Interrupting the Interaction of CK2α/CK2β: Comparative Molecular Dynamics Simulation Studies. Phys. Chem. Chem. Phys 2016, 18, 9202-9210. (33). Sybyl 7.3, T.I., 1699 South Hanley Rd., St. Louis, MI 63144, USA.2006 (34). Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A., A Well-behaved Electrostatic Potential Based Method using Charge Restraints for Deriving Atomic Charges: the RESP Model. J. Phys. Chem. 1993, 97, 10269-10280. (35). Gaussian 09, revision E. 01, Gaussian, Inc.,Wallingford, CT,. 2009. (36). Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A., Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Chem. Inf. Model. 2006, 25, 247-260. (37). Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A., Development and Testing of a General Amber Force Field. J. Comp. Chem. 2004, 25, 1157-1174. (38). Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C., Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct Funct Bioinforma. 2006, 65, 712-725. (39). D. Case, T. D., T. E. Cheatham III, C. Simmerling,; J.Wang, R. D., R. Luo, R.Walker, W. Zhang and K. Merz, U. o. C., San Francisco, AMBER 12. 2012. (40). Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. (41). Darden, T.; York, D.; Pedersen, L., Particle Mesh Ewald: An N. log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98,10089-10093 (42). Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C., Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-alkanes. J. Comput. Phys. 1977, 23, 327-341. (43). Kong, R.; Chang, S.; Xia, W.; Wong, S. T., Molecular Dynamics Simulation Study Reveals Potential Substrate Entry Path into Gamma-Secretase/Presenilin-1. J. Struct. Biol. 2015, 191, 120-129. (44). Amadei, A.; Linssen, A. B. M.; de Groot, B. L.; van Aalten, D. M. F.; Berendsen, H. J. C., An Efficient Method for Sampling the Essential Subspace of Proteins. J. Struct. Biol.Dyn. 1996, 13, 615-625. (45). Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J., GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701-1718. (46). Onufriev, A.; Bashford, D.; Case, D. A., Exploring Protein Native States and Large-scale Conformational Changes with a Modified Generalized Born Model. Proteins: Struct. Funct. Bioinform. 2004, 55, 383-394. (47). Kong, X.; Sun, H.; Pan, P.; Tian, S.; Li, D.; Li, Y.; Hou, T., Molecular Principle of the Cyclin-Dependent Kinase Selectivity of 4-(thiazol-5-yl)-2-(phenylamino) Pyrimidine-5-carbonitrile Derivatives Revealed by Molecular Modeling Studies. Phys. Chem. Chem. Phys. 2016, 18, 2034-2046.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(48). Chohan, T. A.; Qian, H. Y.; Pan, Y. L.; Chen, J. Z., Molecular Simulation Studies on the Binding Selectivity of 2-anilino-4-(thiazol-5-yl)-pyrimidines in Complexes with CDK2 and CDK7. Mol. BioSys. 2016, 12, 145-161. (49). Kuhn, B.; Kollman, P. A., Binding of a Diverse Set of Ligands to Avidin and Streptavidin: an Accurate Quantitative Prediction of their Relative Affinities by a Combination of Molecular Mechanics and Continum Solvent Models. J. Med. Chem. 2000, 43, 3786-3791. (50). Yang, Y.; Shen, Y.; Liu, H.; Yao, X., Molecular Dynamics Simulation and Free Energy Calculation Studies of the Binding Mechanism of Allosteric Inhibitors with P38α MAP Kinase. J. Chem. Inf. Model. 2011, 51, 3235-3246. (51). Huang, Y. Y.; Li, Z.; Cai, Y. H.; Feng, L. J.; Wu, Y.; Li, X.; Luo, H. B., The Molecular Basis for the Selectivity of Tadalafil toward Phosphodiesterase 5 and 6: a Modeling Study. J. Chem. Inf. Model. 2013, 53, 3044-3053.

ACS Paragon Plus Environment

Page 32 of 43

Page 33 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

TOC Graphic

Molecular Simulation Studies on the Binding Selectivity of Type-I Inhibitors in the Complexes with ROS1 versus ALK Yuanxin Tian †,1, Yonghuan Yu†,1,Yudong Shen§, Hua Wan‖, Shan Chang ‡,Tingting Zhang†, Shanhe Wan †, Jiajie Zhang *†

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1.The chemical structure of compounds Figure 1 172x50mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 34 of 43

Page 35 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 2. Plots of RMSD for all the protein backbone atoms (Å) relative to their initial conformation vs simulation time (ns) for four systems. Figure 2 237x137mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. (A)The labeled key loops of ROS1:7c complex. Hinge region is colored by blue; DFG motif is colored by orange; P-loop is colored by red and Alpha-loop is colored by magenta. (B) The diagram of key H-bonds of ROS1:7c complex during the MD simulation. Figure 3 269x182mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 36 of 43

Page 37 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 4. RMSF of each residue of the target ALK(left) and ROS1(right) systems obtained from 20 ns MD simulation. Figure 4 118x47mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. Structure comparison between different system snapshot from MD of: (A) selective(yellow) and dual(cyan) ALK complexes, (B)selective(blue) and dual(orange) ROS1 complexes, (C)dual ALK(cyan) and ROS1(pink) system, (D)selective ALK(yellow) and ROS1(blue) complexes Figure 5 292x223mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 38 of 43

Page 39 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 6. The first and second slowest motion modes(PC1 and PC2) of the ALK-selective and ROS1-selective complexes and heat map Figure 6 211x298mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7. Per-residue binding free energy decomposition of selective systems to find the key residues. Figure 7 257x191mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 40 of 43

Page 41 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 8. The new compounds designed based on our strategy Figure 8 167x101mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The binding modes of compound 14c and SMU-B. Figure 9 289x88mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 42 of 43

Page 43 of 43

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

TOC graphic. Comparison of energy decomposition for key residues selective ALK/ROS1 systems and binding active sites. TOC graphic 723x348mm (300 x 300 DPI)

ACS Paragon Plus Environment