Theoretical Analysis of Activity Cliffs among Benzofuranone Class

Subscriber access provided by READING UNIV

Article

Theoretical Analysis of Activity Cliffs among Benzofuranone Class Pim1 Inhibitors Using the Fragment Molecular Orbital Method with Molecular Mechanics Poisson-Boltzmann Surface Area (FMO+MM-PBSA) Approach Chiduru Watanabe, Hirofumi Watanabe, Kaori Fukuzawa, Lorien J. Parker, Yoshio Okiyama, Hitomi Yuki, Shigeyuki Yokoyama, Hirofumi Nakano, Shigenori Tanaka, and Teruki Honma J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00110 • Publication Date (Web): 07 Nov 2017 Downloaded from http://pubs.acs.org on November 27, 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 59

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

Theoretical Analysis of Activity Cliffs among Benzofuranone Class Pim1 Inhibitors Using the Fragment Molecular Orbital Method with Molecular Mechanics Poisson-Boltzmann Surface Area (FMO+MM-PBSA) Approach

Chiduru Watanabe,†,‡ Hirofumi Watanabe,† Kaori Fukuzawa,*,‡,§ Lorien J. Parker, ¶, # Yoshio Okiyama,† Hitomi Yuki,†Shigeyuki Yokoyama,¶ Hirofumi Nakano,┴ Shigenori Tanaka,∟ and Teruki Honma*,†

†

RIKEN Center for Life Science Technologies, 1-7-22 Suehiro-cho, Tsurumi-ku, Yokohama 230-0045, Japan

‡

Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan

§

Department of Physical Chemistry, School of Pharmacy and Pharmaceutical Sciences, Hoshi University, 2-4-41 Ebara, Shinagawa, Tokyo 142-8501, Japan

¶

RIKEN Structural Biology Laboratory, 1-7-22 Suehiro-cho, Tsurumi-ku, Yokohama 230-0045, Japan 1

ACS Paragon Plus Environment


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 2 of 59

Department of Structural Biology, St. Vincent's Institute, 9 Princes St, Fitzroy, Victoria 3065, Australia

┴

Drug Discovery Initiative, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,

Tokyo 113-0033, Japan

∟

Graduate School of System Informatics, Kobe University, 1-1 Rokkodai, Nada-ku, Kobe 657-8501, Japan

Abbreviations MM: molecular mechanics FMO: fragment molecular orbital QM: quantum mechanics QM/MM: quantum mechanics/molecular mechanics HF: Hartree-Fock MP2: second-order MØller-Plesset perturbation PBSA: Poisson-Boltzmann surface area GBSA: Generalized Born surface area IFIE: inter-fragment interaction energy PIEDA: pair interaction energy decomposition analysis Cpd 1: compound 1 ES: electrostatic EX: exchange-repulsion DI: dispersion interaction CT+mix: charge transfer with higher order mixed terms MM-opt: MM-optimized QM/MM-opt: QM/MM-optimized SAR: structure-activity relationship

2


Page 3 of 59

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: Significant activity changes due to small structural changes (i.e., activity cliffs) of serine/threonine kinase Pim1 inhibitors were studied theoretically using fragment molecular orbital method with molecular mechanics Poisson-Boltzmann surface area (FMO+MM-PBSA) approach. This methodology enabled quantum-chemical calculations for large biomolecules with solvation. In the course of drug discovery targeted Pim1, six benzofuranone-class inhibitors were found to differ only in the position of the indole-ring nitrogen atom. By comparing the various qualities of complex structures based on X-ray, classical molecular mechanics (MM)-optimized, and quantum/molecular mechanics (QM/MM)-optimized structures, we found that the QM/MM-optimized structures provided the best correlation (R2 = 0.85) between pIC50 and the calculated binding energy of FMO+MM-PBSA. Combining the classical solvation energy with the QM binding energy was important for increasing the correlation. In addition, decomposition of the interaction energy into various physico-chemical components by pair interaction energy decomposition analysis suggested that CH-π and electrostatic interactions mainly caused activity differences.

3



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 4 of 59

1. INTRODUCTION Pim1 is a member of the Ser/Thr kinase and plays an important role in cell proliferation.

1,2

Overexpression of Pim1 has been frequently observed in hematologic

malignancies and prostate cancer.

3,4

Potent Pim1 inhibitors with nanomolar-level IC50

values have been reported to induce G1/S cell-cycle arrest and growth inhibition in human prostate cancer and cultured leukemic cells.

1,2,5

To identify more potent and

selective Pim1 inhibitors, we have performed in silico screening based on information about Pim1 crystal structures and known inhibitors. We have discovered several new 6

Pim1 inhibitors, including benzofuranone. Moreover, the X-ray structure of Pim1 and the benzofuranone class inhibitor (IC50: 410 nM, Figure S1) complex was solved to 7

experimentally determine the binding mode (PDB-ID: 3UMX). On the basis of the structural information, an additional basic amine on the left aliphatic ring of the benzofuranone class inhibitor was introduced to interact with Asp128. The IC50 of the 8

modified compound 1 (Cpd 1) was 2 nM.

In the process of exploring

structure-activity relationships (SARs) on the right indole ring of Cpd 1 to improve ADME profile such as solubility and metabolic stability, significant activity changes (more than 100-fold) were observed when only one atom was replaced, as shown in 8

Figure 1. Recently, such an activity change has been called as “activity cliff” 9-14.

4


Page 5 of 59

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


Such activity cliffs have attracted attention from medicinal chemists, because their intentional use creates wider opportunities for hit-to-lead and lead optimization phases.

12, 14

However, when a structural change does not provide any obvious new

interaction, such as a hydrogen bond or steric repulsion, an explanation or accurate prediction of such an activity cliff is very difficult using conventional force-field-based scoring functions. In this study, we tried to quantitatively reproduce the rank order of IC50 of the activity cliff of benzofuranone Pim1 inhibitors using quantum mechanics (QM) together with a solvent continuum model. Molecular mechanics (MM) force fields are widely employed for evaluation of protein-ligand interaction energies because of reduced calculation costs compared with QM. However, most current MM force fields are based on electrostatic interactions with fixed atomic charges and van der Waals interactions with a spherically symmetric Lennard-Jones potential. In contrast, a QM calculation can incorporate the effects of donating and withdrawing electrons and can appropriately deal with cation-π, CH-π, and π-π interactions. The fragment molecular orbital (FMO) method

15–20

can deal with

the electronic states of a protein-ligand complex with hundreds of amino acid residues. Furthermore, the contribution(s) of the interaction(s) from each fragment pair can be more efficiently analyzed by using the inter-fragment interaction energy (IFIE).

19

5



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 6 of 59

Moreover, we have been able to deal with the many-body corrected IFIE (e.g., FMO3-IFIE) and the pair interaction energy decomposition analysis (PIEDA) to 20,21

evaluate the interactions by functional group unit

and the contribution of energies

such as electrostatic and dispersion forces, separately, to facilitate drug design.

22,23

The FMO method has already been applied to analysis of the interactions within various protein-ligand systems such as estrogen receptors, influenza neuraminidases,

30,31

21,24–26

and FK506-binding proteins.

HIV proteases,

27–29

32,33

, where these system

calculated under the gas phase. As reported by Fukuzawa et al.,

24

high correlations

between experimental affinities and calculated binding energies for estrogen receptors and their ligands have been observed, even though the FMO calculation did not incorporate solvation effects. However, because proteins and ligands are surrounded by solvent molecules under physiological conditions, many cases require that solvation effects be considered when evaluating binding affinity.

34–37

For example, when a polar

group of ligands binds to the inner part of a ligand-binding pocket, the prediction of binding affinity should consider the effect of the solvation and de-solvation energy on the polar group. Since the late 2000s, the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) and the molecular mechanics Generalized Born surface area (MM-GBSA) have been proposed to estimate the binding enerigies38-41. These

6


Page 7 of 59

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


MM-based methods together with solvation models significantly improved the correlation between calculated binding energies and experimental IC50 values. However, the MM-based methods can not sufficiently take QM-based interaction energies such as dispersion force and charge transfer into consideration. In a recent study, solvation effects have been incorporated into the FMO method via explicit consideration of water molecules 44–47

polarized continuum model

42,43

or use of implicit solvation models such as the

and Poisson-Boltzmann equation.

48

More convenient

approaches, which combine the interaction energy obtained from a QM calculation with the solvation energy obtained from a MM calculation, have also been undertaken to predict binding affinity.

32, 49–52

In addition to the methods used to calculate

interactions and solvent effects, attention should be given to structure preparation by appropriate QM-based or other methods to perform more accurate calculations. The use of MM optimization often fails to reproduce the appropriate hydrogen bonding distance, and the coordinates of such an optimized structure are inadequate, for example, in the QM framework. High-resolution, three-dimensional structure coordinates, which are based on optimizations with quantum/molecular mechanics (QM/MM)

53

54

or FMO,

are thus critically important, because in the QM calculation,

7



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 8 of 59

0.1–0.5 Å differences of atomic coordinates cause significant differences in the interaction energy.

55

In this paper, we used the FMO method to reproduce the SAR between Pim1 and its six benzofuranone-class inhibitors that show activity cliff, which would be caused by differences of electronic structures accompanied by the replacement of one or two 8

atom(s). The six inhibitors differed only in the position of the nitrogen atom in the indole ring (Cpd 1) or azaindole ring (Cpds 2–6), as shown in Figure 1A. The SAR of these inhibitors provides a typical example of the activity cliff, because the tiny differences in the structures of these compounds resulted in more than a 200-fold range of activity. Figures 1B and 1C indicate structures around the ligand-binding pocket of Pim1 and the Cpd 1 complex. Near the indole ring, there were important hydrogen bonds with the hinge region (Glu121) and CH-π interactions

56,57

with hydrophobic

residues (Leu44, Ala65, and Leu174) in β sheets. Although the indole ring of ligands, which is associated with replacement of the nitrogen atom, binds to the inner part of the ligand-binding pocket, part of the indole ring and other substructures of a benzofuranone and piperazine partially contact solvent molecules (Figure S3). We tried to quantitatively reproduce the SAR with the FMO method by incorporating solvation effects with the molecular mechanics Poisson-Boltzmann surface area

8


Page 9 of 59

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


(FMO+MM-PBSA) method

49,50

. This method was dividing up the calculation

according to thermodynamic cycle shown in Figure 2. In addition, we compared the binding affinities estimated with three types of structure preparations (X-ray structures as is [X-ray], the MM-optimized structures [MM-opt], and the structures optimized by the combination of QM and MM calculations [QM/MM-opt]). We then discuss which calculation method (FMO and MM with/without PBSA) and structure preparation method improved the correlation between the experimental IC50 values and the estimated binding energies. We also identified key residues and key energy components that contributed to the affinity difference by using the many-body corrected IFIE and an analysis of its decomposition.

2. COMPUTATIONAL METHODS Based on the flow chart in Scheme 1, we evaluated protein-ligand binding energies by the FMO method and solvation free energy by the MM-PBSA method. This scheme gave the best correlation between calculated binding energies and IC50 values in our study. Starting from given X-ray crystal structures, three types of geometry were prepared for Pim1-inhibitor complex structures to examine the influence of geometry optimization on the correlation of interaction energy with IC50 values: (1) individual

9



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 10 of 59

X-ray complex structures were used for each ligand with minimal molecular modeling (X-ray); (2) complex structures were modeled based on one-template complex structures with MM optimization (MM-opt); and (3) complex structures were further optimized with a QM/MM calculation (QM/MM-opt). With these structures, we used FMO calculations at MP2/6-31G* levels in vacuo to evaluate protein-ligand interaction energies, followed by a classical MM-PBSA calculation to evaluate the solvation/desolvation energy and the internal energy by using the PBSA term. To evaluate the deformation energy of the free ligand, solvated ligand structures were prepared with MM-based conformational searches and QM-based optimization and energy calculations in implicit water. Details of the calculation procedure of FMO and FMO+MM-PBSA methods are discussed in sections 2.4 and 2.5. 2.1. Dataset. In this work, Pim1 kinase and its six benzofuranone-class of inhibitors, 8

with IC50 values of 2–447 nM, were selected for the calculations (Figure 1A). Because the differences of these inhibitors are very small with respect to the position of the nitrogen atoms in the azaindol ring, a comparison of their activities is a typical example of an activity cliff with a 200-fold difference in activity. The X-ray structures 8

of Pim1 were experimentally obtained for Cpds 1, 3, 5, and 6. Data collection and refinement statistics for crystal complexes of Pim1 kinase domain in complex with

10


Page 11 of 59

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


compounds 1, 5 and 6 are provided in Table S1. Statistics for the compound 3 complex has been published previously8.

The PDB accession codes for the structures of the

compound 1, 3, 5 and 6 complexes are 5VUC, 3UMW, 5VUA and 5VUB, with resolutions of 2.00, 2.08, 2.20, and 2.00 Å, respectively (Table S1). Superimposed X-ray structures complexed with these four compounds around the ligand binding pocket indicated that the binding modes of these compounds were almost identical (Figure S3). The implication was therefore that the differences of the electron distributions in the ligands primarily contributed to the activity difference. Because the complex structures have not yet been solved for Cpds 2 and 4, their complex model structures were constructed as described in the next section. 2.2. X-ray crystal analysis of complexes between Pim1 and Cpds 1, 5, and 6. Pim1 expression and purification for crystallization. The kinase domain of human Pim1 was synthesized by the Escherichia coli cell-free system using the dialysis 58-60

method

. Details of the purification and expression methods have been reported

8

previously . Pim1 crystallization. Human Pim1 was crystallized under the conditions previously 61

reported by Jacobs et al. (2005) . Soaking involved the addition of the solid compound in a well solution containing 3% DMSO for 24 h. The use of the soaking

11



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 12 of 59

method to determine the crystal structures of Pim1 with this range of inhibitors has 7

been discussed in detail in Parker et al., 2012 . Crystals were cryoprotected with 30% glycerol and flashcooled in liquid nitrogen. Data were on the BL26B1 beamline at SPring-8, Harima, Japan. 62

Pim1 structure solution. MOSFLM was employed for processing and integration 63

and SCALA was used for scaling . The structures were solved using an in-house 8

64

model of Pim1 based on 3UMW and the program PHENIX . The space group is P65 and there is one molecule in the asymmetric unit. Refinement was performed using 65

PHENIX. Model building was accomplished with Coot . Accession numbers. The coordinates and structure factors of the complexes have been deposited in the PDB with accession codes 5VUC, 5VUA and 5VUB for Cpds 1, 5 and 6, respectively. X-ray structure solution. The crystals of the complexes of Pim1 with Cpds 1, 5 and 6 diffracted to 2.0, 2.2 and 2.0 Å resolutions, respectively. The electron densities for the compounds were immediately evident in the Fo–Fc and 2Fo–Fc maps. Compounds 1, 5 and 6 were fitted with occupancies of 100% and had average B-factors of 28.1, 32.0 and 31.0 Å2, respectively. The average protein B-factors were 37.6, 34.4 and 31.3 Å2 for all non-hydrogen atoms included in refinement.

12


Page 13 of 59

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


2.3. Molecular Modeling. X-ray structures. Complex structures were prepared based on the X-ray crystal structures. The amino acid residues of #37-305 were employed because coordinates of these residues are commonly given by X-ray crystallography. All crystallization waters were removed. Because the complex X-ray structures of Cpds 1, 3, and 5 have missing residues and atoms, we replaced the missing atomic coordinates with those of the complex structure of Cpd 6. These complemented atoms and added hydrogen atoms were optimized concurrently using 66

the Amber99 force field

67

Computing Group Inc.).

with the MOE graphical software package (Chemical

Because of the similarity of the ligand-binding mode of

these complexes, the locations and directions of side chains around the ligands were also very similar (Figure S3). On the basis of these structural observations, (1) the two hydrogen bonds with the main chain C=O of Glu121 and the side chain NH3+ of Lys67 (Figure 1B) and (2) the several CH-π interactions with the side chains of Leu44, Val52, Ala65, Leu174, and Ile185 (Figure 1C) were thought to be important interactions with compounds. We refer to such prepared structures as “X-ray” structures in this paper. Modeling of complex structures based on the template X-ray structure of Cpd 1. In the case of “X-ray”, we used X-ray complex structures of Pim1 with Cpds 1, 3, 5, and

13



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 59

6 as is. On the other hand, for the “MM-opt” and “QM/MM-opt” cases, all the structures were modeled using one template structure. As the template structure, we selected the X-ray structure of Pim1-Cpd 1 complex because it showed the highest resolution (2.00 Å) and Cpd 1 is the most potent compound (IC50: 2 nM) among published complex structures. Starting from the X-ray structure of complex between Pim1 and Cpd 1, we replaced the appropriate carbon atom of Cpd 1 with a nitrogen atom to form Cpds 2, 3, 4, 5, and 6. The structure optimizations of initial structure were carried out under the partially restraint conditions (Figure 3) in order not to largely change from template structure based on the X-ray crystal structure. We then performed the MM optimization with the Amber99 force field. During the optimization, the atoms of indole/azaindole ring as belonged to compound colored red in Figure 3A were unconstrained, and the other atoms were fixed during the optimization. This modeling was performed with the MOE. We refer to the structures built with these procedures as “MM-opt” structures. Next, we performed the QM/MM geometry optimization starting with the MM-opt structures. We employed the ONIOM-QM/MM method at the HF/6-31G*:UFF level of theory using the Gaussian09 program package.

68

The QM region was defined as all

atoms of the ligand and partial atoms of main-chain on Glu121 and Arg122. Hydrogen

14


Page 15 of 59

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


bonds joined these residues with the indole/azaindole ring of the ligand. The MM region was defined by the other atoms of the receptor. In addition to the atoms that were not constrained in the MM-opt calculation, the red C, H, O, and N atoms based on indole/azaindole ring of compound and main chain between Glu121 and Arg122 in Figure 3B were not constrained, and the other atoms were fixed in the QM/MM geometry optimization. We refer to this QM/MM optimized structure as “QM/MM-opt” structures. The hydrogen bond lengths between ligand and amino acid residues of Pim1, r1–r4 (Figure 3), are summarized in Table S2. The QM/MM-optimized hydrogen bond lengths of r1 were calculated to be 1.76–1.85 Å, whereas those of the MM optimization had less variation, 1.69–1.73 Å. The other hydrogen bond lengths, r2, r3, and r4, were fixed at the MM-opt structure of Cpd 1; 1.75, 1.59, and 2.32 Å, respectively. Ligand structures in water. To realize Pim1-ligand binding in solution, it was necessary to consider the ligand structure in water and its deformation to a complex structure. For the calculation of solvation/desolvation energies of the ligand, we performed an MM-level conformational search to select representative and energetically stable conformations. This step was followed by QM geometry optimization in an implicit water model. First, we carried out a systematic

15



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

conformational

search

using

the

MMFF94x

force

field

Page 16 of 59

with

MOE;

the

root-mean-square deviation limit was 0.15 Å, and solvation effects were taken into account with the generalized Born model. The most stable conformations for Cpds 1– 3, 5, and 6 were structures that were flipped at the bond between the benzofuranone ring and the azaindole ring, in contrast to the conformations in X-ray complex structures (Figure S4). However, in the case of Cpd 4, the most stable conformation at the MM level was not the flipped structure. The second most stable structure of Cpd 4 in the MM method was the flipped structure. The stable conformations obtained with the MM method were then optimized at the HF/6-31G* level with polarizable 68,69

continuum model (PCM) calculations using Gaussian 09 software.

The results

indicated that the flipped conformations of Cpd 4 were more stable in an aqueous solution: the energies of flipped conformations were 1.3 kcal/mol lower on average. Taken together, the ligand conformation in water was determined by the following procedure. (1) Ligand conformation searches were performed by MM calculation. (2) For two stable conformations (no-flip and flip structures) geometry optimization at the PCM-HF/6-31G* level was performed. Among them, we employed the most stable conformation (flipped structure) as the ligand structure in solution for calculation of ligand deformation energy (Figure 2).

16


Page 17 of 59

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


2.4. The FMO method and FMO-based interaction energy analysis. In the FMO method, a protein-ligand structure is divided into small fragments, and the total energy of the system is constructed by summing the partial energies obtained with molecular orbital (MO) calculations of each fragment and its combinations, taking into account the environmental electrostatic potential (ESP) from other fragments. Whereas the two-body approximation (abbreviated to FMO2) is generally used for a protein-ligand interaction analysis at the compound/amino acid level of resolution, the three-body approximation (abbreviated to FMO3) has been employed for such analyses at a 20,21

functional group level of resolution;

total energies with these approximations are

written for fragment indices I, J, and K in the following equations:

~ FMO2 Etotal = ∑ EI′ + ∑∆EIJ I

(1)

I >J

~

FMO3 FMO2 Etotal = Etotal +

∑ ∆E

IJK

(2)

I > J >K

~ where E I′ is the monomer energy without the environmental ESP contribution; ∆E IJ

( = ∆ E IJFMO2 ) is the two-body correction energy, the so-called inter-fragment interaction ~ energy (IFIE or FMO2-IFIE); and ∆E IJK is the three-body correction energy. By

using these correction terms, the IFIE with three-body correction (FMO3-IFIE) can be defined as follows:

17



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

∆EIJFMO3 = ∆EIJFMO2 + In addition, PIEDA

Page 18 of 59

1 ~ ∆EIJK ∑ 3K

22,23

(3)

was used to analyze the energy components of FMO2-IFIE,

∆E IJFMO2 : electrostatic (ES), exchange-repulsion (EX), dispersion interaction (DI), and

charge transfer with higher order mixed terms energies (CT+mix), where

~ ~ ~ ~ ~ ∆EIJ = ∆EIJES + ∆EIJEX + ∆EIJCT+mix + ∆EIJDI..

(4)

In the FMO3-IFIE calculations at the MP2 level, only the DI term, which is synonymous to electron correlation energy, was decomposed as follows:

~ ~ ~ ~ ∆EIJMP2 = ∆EIJHF + ∆EIJCORR = ∆EIJHF + ∆EIJDI..

(5)

Therefore, in accord with the interaction energy analysis in section 3.2, the FMO2-IFIE was divided into four PIEDA components and the FMO3-IFIE into two components to reveal the difference of the electronic mechanism that invokes an activity cliff. To obtain ligand-binding energies, we summed the IFIE of all ligand-residue pairs. The FMO contribution of the ligand-binding energy in vacuo was defined using FMO3-IFIE interaction energies as follows: int ∆EFMO = ∑ ∆EIJFMO3 .

(6)

I =ligand J ≠I

In this work, all FMO calculations were performed at the FMO3-MP2/6-31G* level 70

with the Cholesky decomposition approximation (called CDAM)

by using the 18


Page 19 of 59

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


ABINIT-MP program.

18,71,72

The Mulliken atomic orbital approximation (called

ESP-AOC) was applied to the calculations for consideration of the environmental ESP.

19

Each protein-ligand system was divided into 504 fragments. The main chain

and side chain of each residue were treated as different fragments, and the ligand was divided into two fragments as shown in Figure 4. In the IFIE analysis of each amino acid fragment, the residue name subscripts “M” and “S” designate the main and side chain fragments, respectively (e.g., Arg122M and Lys67S). IFIE and PIEDA analyses 71,72

and their visualizations were performed with BioStation Viewer.

2.5. Binding Free Energy in Solution: The FMO+MM-PBSA Method. Although the FMO calculation gave us accurate information about protein-ligand interactions based on enthalpic effects, estimation of solvation free energy is necessary to explain 73

SARs and activity cliffs. In contrast, the MM-PBSA method

is often used to evaluate

protein-ligand binding free energies with incorporation of solvation effects. However, treatment of weak interactions, such as hydrogen bonds, charge transfers, and CH-π and π-π interactions are not accurate enough in the MM energy calculations. Complementary use of the MM and QM methods can provide a solution if the MM part of the MM-PBSA energy function is replaced by FMO-QM counterpart. Here we employed the solvation energy obtained with the MM-PBSA calculation to describe

19



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 59

desolvation during complex formation. In addition, we employed the value obtained with the MM-PBSA method for the internal energy of the FMO+MM-PBSA. Whereas 32

in the MM-PBSA calculation the entropy term is added when needed , we here neglected the entropy term because the entropic contributions of the six compounds in free energy showed approximately equal in the complexation (Table S3); in the case of the Pim1 dataset, the difference of structure is due to only the single substitution of nitrogen for carbon. We

defined

the

protein-ligand

binding

free

energy

obtained

with

the

FMO+MM-PBSA method (Figure 2) by using the FMO-based interaction energy in def int , the MM-based deformation energy, ∆E MM , that arises from vacuo (eq 6), ∆ E FMO

solv conformational changes, and desolvation free energy ∆G MM − PBSA as follows:

bind int def solv ∆ G FMO = ∆ E FMO + ∆ E MM + ∆ G MM − PBSA .

(7)

Here, effect of electron relaxation on ligand binding could not be incorporated into int def summation of IFIE, ∆ E FMO . Thus, the ligand deformation energy ∆E MM was treated

by MM method to remove electron relaxation from evaluation of ligand binding energies. The deformation energy is defined by the difference in MM-based internal internal internal energies of a ligand between its complex form, Eligand(C) , and solvated one, E ligand :

def internal internal . ∆E MM = Eligand(C) − Eligand

(8)

20


Page 21 of 59

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 MM-based desolvation free energy is given as follows:

(

solv solv solv solv ∆GMM − PBSA = Gcomplex − G protein + Gligand

)

(9)

solv solv solv , Gprotein , and Gligand are the solvation energies of the complex, protein, where G complex

and ligand, respectively, and all are derived from the MM-PBSA calculation. In this supermolecule calculation, geometry of free protein was fixed at the complex structure, but regarding the geometry of free ligand, the structure obtained by MM conformation search and QM optimization in water were used. These equations correspond to the following equations in the general MM-PBSA scheme that are used to obtain the MM-based binding free energy: bind bind solv ∆ G MM = ∆ E MM + ∆G MM − PBSA .

(10)

The enthalpic term is written as follows:

(

)

bind MM MM MM ele vdW def ∆E MM = Ecomplex − E protein + Eligand = ∆EMM + ∆E MM + ∆EMM

(11)

MM MM MM where E complex , Eprotein , and E ligand are the MM-based total energies for the

ele complex, protein, and ligand, respectively, and ∆EMM

vdW and ∆ E MM

are the

non-bonding electrostatic and van der Waals interaction energies between the protein and a ligand, respectively. Finally, we defined the relative binding free energy for each compound as the difference in binding free energy from Cpd 1:

21



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 22 of 59

bind bind bind ∆∆ GCpd X = ∆GCpd X − ∆GCpd 1 ( X = 2 − 6 )

(12)

We carried out this calculation by using both the FMO+MM-PBSA (eq 7) and the MM-PBSA (eq 10) methods to estimate the affinities of the six compounds to the Pim1 protein. All of the MM and related PBSA energies were calculated with the mm_pbsa.pl program in the AMBER 10 package

74

75

using the Amber ff99SB force field 76

proteins and the general Amber force field (GAFF)

77,78

with the AM1-BCC

for

charge

for ligands. Dielectric constants of solute molecules and solvent water with a 1.4-Å probe radius were set to 1.0 and 80.0, respectively, to evaluate solvation free energy.

3. RESULTS AND DISCUSSION 3.1. Binding energies and their correlations with IC50 values. Based on the derived equations for the binding energy calculations in section 2, three types of bind bind int relative binding energies, − ∆∆GFMO , − ∆∆ E FMO , and − ∆∆GMM , were evaluated

with the experimental pIC50 values of Cpds 1–6. The three types of optimized geometries prepared in section 2.3, X-ray, MM-opt, and QM/MM-opt, were used in each binding energy calculation. All correlation diagrams between pIC50 values and

22


Page 23 of 59

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


bind int relative binding energies are plotted in Figure 5. Numerical values of ∆GFMO , ∆EFMO

bind , and ∆GMM are listed in Table 1 and Table S4.

bind The relative binding free energies, − ∆∆GFMO , determined by the FMO+MM-PBSA

method are shown in Figures 5A–C. In the case of the X-ray structures, the correlation with pIC50 was relatively poor (R2 = 0.49). In particular, although Cpds 1 and 3 have similar IC50 values, the FMO+MM-PBSA calculation yielded significantly different energies (ca. –18 kcal/mol) that were reversed in relation between pIC50 and binding energy. In the case of the MM-opt structures, the correlation with pIC50 values was worse (R2 = 0.24). It was difficult to discriminate between the higher pIC50 group (Cpds 1–3) and the lower pIC50 group (Cpds 4–6). The implication was that the difference of activities could not be appropriately evaluated with the binding energies obtained from the QM calculations with the MM optimized structures. In contrast, when the QM/MM-opt structures were used, there was a good correlation with pIC50 values (R2 = 0.85). The reason for the good correlation could be improvement of the hydrogen bond length, r1 (Figure 3 and Table S2), between the N-H of the indole/azaindole ring of the ligand and the C=O of Glu121 by use of the QM/MM optimization. The r1 values of Cpds 1–6 varied from 1.76 to 1.85 Å in the QM/MM-opt structures, whereas those of the MM-opt structures were in the range of 1.69–1.73 Å

23



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 24 of 59

(Table S2). The relatively flexible hydrogen bond lengths of the QM/MM-optimized structures may be attributed to the better estimation of binding energies. In particular, bind the − ∆∆GFMO outliers for Cpds 2 and 6 in the MM-opt structures were corrected in

the QM/MM-opt structures (see Figures 5B and 5C). The application of the QM/MM-opt structures consequently improved the classification of compounds into groups with higher (Cpds 1, 2, and 3) and lower (Cpds 4, 5 and 6) experimental activities. At this level of correlation, we could quantitatively reproduce the IC50 value of each compound. To study the effects of the solvation energy on the binding free energy, we compared the results without the solvation term of the PBSA. Ligand-binding energies int based on pure FMO calculations, − ∆∆E FMO , which include only the contribution of

enthalpy, are shown in Figures 5D–F. The ligand binding energies determined by only the FMO resulted in a lower correlation with pIC50 values, even for the QM/MM-opt structure (R2 = 0.67). Such results suggest that the incorporation of solvation free energy was very important, even for a dataset with subtle differences of compound structures, where cancellation of solvation energy was expected to occur in the relative binding energy calculation. We compared the FMO theoretical levels between Hartree-Fock (HF) and second-order MØller-Plesset perturbation theory (MP2) (Table

24


Page 25 of 59

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


1) and examined the contribution of ligand fragments (Figure 4B). There was a moderately positive correlation between both the HF and MP2 levels and pIC50 values. Although the R2 values for both the total ligand, Fragment (1)+(2), and Fragment (1) at the MP2 level were improved from the MM-opt to the QM/MM-opt results, the R2 value for Fragment (2) was unchanged. These results may have been caused by the fact that the geometry optimization was only partially performed around Fragment (1). Next, we compared the accuracy of the protein-ligand binding energy between the FMO and MM methods. Figures 5G–I show the correlation between the pIC50 values bind and the estimated − ∆∆GMM obtained from the MM-PBSA calculations. The

MM-PBSA calculations resulted in negative correlations for the X-ray and MM-opt structures; the correlation was poor (R2 = 0.34) for the QM/MM-opt structure. It was hard to distinguish between the lower and higher pIC50 groups.

The energy

ele vdW def solv components of the MM-PBSA, ∆E MM , ∆EMM , ∆EMM , and ∆GMM - PBSA , are listed

ele in Table S4. Only the electrostatic component, − ∆EMM , was positively correlated

with pIC50 values for all types of structures. These results suggest that there would be some difficulty in evaluating van der Waals interactions such as CH-π and π-π interactions by using the MM-PBSA method.

25



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 59

Estimation of binding energy by using the FMO+MM-PBSA method with QM/MM-optimized geometries therefore gave us results that most accurately reproduced experimental activities (Figure 5C, R2 = 0.85). The only outlier in Figure 5C was Cpd 1. Because only Cpd 1 does not have an additional nitrogen atom with a lone pair on the indole ring, there are several reasons can be considered for its outlier position; solvation energy was not accurate enough because PBSA term was calculated with MM-based method, not with QM-based method; the rearrangement of conserved water molecules would occur. Further consideration should be investigated in a future study including combination with molecular dynamics simulation. 3.2. Analyses of interactions between Pim1 and ligands. We performed interaction energy analyses of fragment units for Pim1-ligand complexes using the IFIE and PIEDA energies. In this paper, the main chain and side chain of each amino acid residue were assigned to different fragments, and the ligand was divided into two fragments. Such a fragmentation scheme would provide a detailed analysis in terms of the approximate contribution of each functional group. By employing the PIEDA 22,23

method

, we decomposed interaction energies into four components, ES, EX,

CT+mix, and DI, between the ligand and the amino acid residues around the ligand binding pocket. Electrostatic potentials of ligands were also examined.

26


Page 27 of 59

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


Interaction energies determined by IFIE and PIEDA. Figure 6 and Table 2 show the results of the FMO3-IFIE and PIEDA calculations. Whereas the HF method evaluated mainly electrostatic interactions, including hydrogen bonds, the MP2 method incorporated the effect of dispersion interactions in addition to the other interactions (eq 5). In the FMO3-IFIE analysis at the MP2 level (Figure 6A), strong attractive interactions between the ligand fragment (1) and surrounding charged amino acid residues were observed for Lys67S (Lys67 side chain), –34.9 kcal/mol, and for Arg122M (main chain carbonyl oxygen of Glu121), –17.3 kcal/mol. In contrast, the interactions between the ligand fragment (1) and Glu89S and Asp186S were repulsive (17.3 and 7.0 kcal/mol, respectively). In addition, there were attractive interactions involving hydrophobic residues (Leu44S, Phe49S, Val52S, Leu174S, and Ile185s), –1.7, –2.1, –4.3, –2.2, and –2.7 kcal/mol, respectively, whereas the FMO3-IFIE analysis at the HF level showed slightly repulsive interactions for these fragment pairs (data not shown). The implication is that consideration of the electron correlation energy, which is regarded as dispersion interaction energy, is important to the description of interactions with hydrophobic residues. These findings from IFIE analyses are consistent with observations from X-ray complex structures (Figures 1B–C); there may be two hydrogen bonds between the ligand fragment (1) and Lys67S fragment and

27



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 28 of 59

Arg122M fragment (Figure 1B); many CH-π interactions have also been suggested between the ligand and the hydrophobic residues (Figure 1C). With respect to the int correlation between the interaction energy, − ∆EFMO , and pIC50, the MP2 energy and

electron correlation energy were shown to have similar R2 values, the suggestion being that the dispersion term is important for characterization of the SAR. To elucidate the important energy component, we performed an energy decomposition analysis with the FMO2-PIEDA. Although the FMO3-MP2 energies are more quantitatively accurate than those of FMO2-MP2, we found that the implications of the results with these two levels of theory were qualitatively similar in Table 2. Figure 6B is a color-coded representation of the most significant energy components for ligand interactions with each main/side chain fragment. This analysis confirmed that the main components of the interaction energy were the ES term with hydrophilic residues and the DI term with hydrophobic residues.

Table 2 shows the

correlation coefficients between each energy component and pIC50. For fragment (1), both the ES and DI terms were better correlated with pIC50 (R2 = 0.65 and 0.62, respectively) than was the total interaction energy (R2 = 0.53); in the case of the correlation for the whole ligand, fragment (1)+(2), the correlations for both the ES and DI terms were even better (R2 = 0.74 and 0.74, respectively). Furthermore, even better

28


Page 29 of 59

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


correlations were obtained for the sum of the ES and DI terms, i.e., R2 = 0.75 for fragment (1) and R2 = 0.80 for fragment (1)+(2). Tables S5 and S6 show results of a detailed analysis at the residue level of the ES and DI components, respectively, between ligand fragment (1) and the surrounding fragments of amino acid residues (Figure 6B). In these tables, we selected the residues with large deviations among the six ligands. The values of amino acid residue fragments within a distance of 4.5 Å from a ligand are listed in Table S5 and S6. In Table S5 and Figure 6B, the ES interactions forming hydrogen bonds between Lys67S and Arg122M with the ligand were the strongest (–36 to –40 kcal/mol and –17 to –24 kcal/mol, respectively). However, the correlations with pIC50 were negative for these two fragments, although the ligand-binding pose may be determined by the strongly attractive interaction of hydrogen bonds. Positive correlations with pIC50 were obtained for Arg122S, Asp186S, Pro123M, Leu120M, Ala65M, and Ile104M; their ES energies indicated interactions with the fragment (1) of the ligand that ranged from moderately attractive to repulsive, –2 to +11 kcal/mol; no hydrogen bonds were observed. These results suggest that the positive/negative charges of these side chain fragments and the polarization of the main-chain fragments contributed to the ES interactions and to the characteristics of the ligand-binding affinities. For the DI

29



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 30 of 59

interactions shown in Table S6 and Figure 6B, several hydrophobic residues, Ile104S, Leu120S, and Leu44S, provided positive correlations with pIC50 for fragment (1), whereas the correlations for the hydrophobic residues Val52s, Leu174S, and Ile185S were relatively weakly negative. These interactions with hydrophobic side chains suggested the existence of attractive CH-π interactions (see Figure 1C), and the difference in strength of the CH-π interactions contributed to the differences of binding affinities. In particular, Cpd 1 showed the strongest DI interactions with the Leu44S fragment (Table S6) and the largest total DI energies (Table 2), the suggestion being that the CH-π interactions were one of the causes of the very small IC50 value (2 nM) of Cpd 1. The charge distributions and electrostatic potentials of the indole/azaindole ring were consistent with the significance of the ES interactions. Natural population analysis (NPA) and electrostatic potential mapping were performed at the HF/6-31G* level. While the total NPA charges of the indole/azaindole ring in Cpds 1–6 remained around zero (qindole = 0.05e–0.09e), the NPA charge on individual atoms at the introduced nitrogen atom and neighboring carbon atoms varied greatly (see Figure S5); the atomic charge decreased at the introduced nitrogen atom (∆q ~ –0.3e) and increased at the neighboring carbon atoms (∆q ~ 0.3e). The introduction of a nitrogen atom also

30


Page 31 of 59

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


affected the electrostatic potentials of ligands. Totally positive electrostatic potentials around the ligand were derived from the positive charge of the ligand, and locally negative electrostatic potentials were observed around the carbonyl oxygen atom and the introduced nitrogen atom (Figure 7). Such differences of electrostatics for individual compounds should affect the ES interaction between the ligand and surrounding amino acid residues. We performed a detailed analysis of the PIEDA component between the ligand and the nearest fragments of amino acid residues (Figure 8). The ES interaction of Arg122S was found to be basically repulsive, but slightly attractive for Cpds 2 and 5. In these two compounds, a nitrogen atom with a negative charge was introduced at the closest position to Arg122 and therefore caused an electrostatic attraction with the positive charge of the Arg122 guanidyl group. For Ile104M (the carbonyl group of the peptide bond between Val103 and Ile104), Cpd 3, with a nitrogen atom at the position closest to Ile104M, was shown to have an attractive ES interaction compared to the other compounds. The distances between the introduced nitrogen atom of the ligand and the positively charged atom of Arg122S/Ile104M ranged from 6 to 8 Å, in a range intermediate between the first and second layers of the protein around the ligand. These results suggest that weak

31



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 32 of 59

electrostatic interactions of a fragment that did not directly connect to the ligand contributed to the IC50 difference. In summary, the PIEDA analysis showed that multiple influences of both DI and ES interaction energies between indole/azaindole rings and surrounding residues affected the SAR between Pim1 and the six compounds.

4. CONCLUSIONS In this paper, we successfully reproduced the SAR for the complexes between Pim1 and its six inhibitors (Cpds 1–6), a target system that includes activity cliffs, by a combination of QM calculations of the whole system and semi-empirical solvation-free energies, the FMO+MM-PBSA method. A good correlation (R2 = 0.85) was obtained between experimental pIC50 values and calculated binding free energies, whereas the classical MM-PBSA calculation failed to quantitatively reproduce the SAR (R2 = 0.34). Consideration of solvation effects (PBSA terms) in the former was important for obtaining higher correlations with pIC50 values. Concerning the quality of the structures used for the FMO+MM-PBSA calculations, the calculated binding free energies that use the bare X-ray structures and MM-opt structures showed poorer correlations with pIC50 values (R2 = 0.49 and 0.24, respectively). Use of these structures made it impossible to discriminate between higher and lower IC50 groups. In 32


Page 33 of 59

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


contrast, the QM/MM-opt structures provided a much higher correlation (R2 = 0.85). As suggested in this case, the X-ray structures did not always provide sufficiently accurate bond lengths for structure-based interaction energy analysis by the FMO method. Through FMO-based interaction energy analyses, both ES and DI interactions were suggested to be important factors for the SAR. Although two strong hydrogen bonds were observed for ligands with Lys67 and Glu121 (carbonyl oxygen of backbone), the strengths of these hydrogen bonds were not correlated with IC50 values. Not only specific hydrogen bonding residues but also surrounding residues with ES interactions contributed to the characteristics of the ligand-binding affinities. It became apparent that electrostatic potential modulation caused by differences of the nitrogen atom position on the indole/azaindole ring and CH-π interactions between indole/azaindole rings and surrounding amino acid residues affected the SAR. Interaction energy analyses based on the FMO method and decomposition analyses by PIEDA were quite useful for clarifying characteristics of ligand binding and the importance of each amino acid residue in a ligand. In the present study, we employed ligand deformation and solvation free energy terms based on MM force field, in spite of including theoretical inconsistency of

33



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 34 of 59

mixing MM-based energy into FMO energy. We already examined to replace MM-based ligand deformation energies by QM-based ones; however, the correlation between pIC50 and FMO+MM-PBSA energy values was not improved in this dataset (Supporting Information section F). Thus, binding free energy will be more reliably estimated by directly being coupled with QM-based solvation and deformation energy calculations including QM-polarizable medium79, which is under development as full FMO-PBSA. We also neglected entropic terms, because the structural degrees of freedom for compounds used here were almost identical. However, improved treatment with incorporation of an entropic term that is adaptable to the QM method is needed to extend the applicability of the method to more diverse compounds. These quantum mechanical approaches should provide a powerful tool for in silico rational drug design in the context of structure-based drug design.

ASSOCIATED CONTENT Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI:XXXX. Section A shows initial hit of the benzofuranone inhibitor (Figure S1). Section B concerns the CH-π (CHPI) analysis method (Figure S2). Section C provides 34


Page 35 of 59

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


information on X-ray crystal structures of Pim1 and inhibitor complexes (Figure S3 and Table S1), bond lengths of the hydrogen bond network between ligand and amino acid residues of Pim1 (Table S2), and stable structures of ligands in complex and with solvation (Figure S4). The entropy term was calculated by the stable structure of each ligand in water (Table S3). Section D includes MM-PBSA results (Table S4). Section E is a NPA charge analysis (Figure S5). Section E describes the PIEDA (Tables S5 and S6). Section F compares FMO+MM-PBSA binding energy including MM-based ligand deformation energy with that including QM-based ligand deformation energy (Table S7).

AUTHOR INFORMATION Corresponding authors

*K. F.: E-mail: [email protected].

*T. H.: E-mail: [email protected]

ACKNOWLEDGEMENTS The authors would like to thank Prof. Yuji Mochizuki and Dr. Toru Sengoku for technical support and help with preparation of the manuscript, and Dr. Noriaki Hashimoto, Dr. Daisuke Takaya, and Dr. Tomohiro Sato for fruitful discussions. This

35



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 36 of 59

research was a part of the “Research and Development of Innovative Simulation Software” project supported by the Research and Development for Next-generation Information Technology of the Ministry of Education, Culture, Sports, Science and Technology (MEXT). This work was supported by JSPS KAKENHI Grant Number JP15K05397, JP15K21632 and 26460035. This research is supported by the Platform Project for Supporting in Drug Discovery and Life Science Research (Platform for Drug Discovery, Informatics, and Structural Life Science) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), and Japan Agency for Medical Research and Development (AMED). This research in part made use of the computational resources of the RICC (RIKEN Integrated Cluster of Clusters) and the π-VizStudio at the Kobe university. References 1. Brault, L.; Gasser, C.; Bracher, F.; Huber, K.; Knapp, S.; Schwaller, J. PIM Serine/Threonine Kinases in the Pathogenesis and Therapy of Hematologic 2.

Malignancies and Solid Cancers. Haematologica 2010, 95, 1004–1015. Magnuson, N. S.; Wang, Z.; Ding, G.; Reeves, R. Why Target PIM1 for Cancer

3.

Diagnosis and Treatment? Future Oncol. 2010, 6, 1461–1478. Amson, R.; Sigaux, F.; Przedborski, S.; Flandrin, G.; Givol, D.; Telerman, A. The Human Protooncogene Product P33pim Is Expressed During Fetal Hematopoiesis

4.

and in Diverse Leukemias. Proc. Natl. Acad. Sci. U.S.A. 1989, 86, 8857–8861. Dhanasekaran, S. M.; Barrette, T. R.; Ghosh, D.; Shah, R.; Varambally, S.; Kurachi, K.; Pienta, K. J.; Rubin, M. A.; Chinnaiyan, A. M. Delineation of

5.

Prognostic Biomarkers in Prostate Cancer. Nature 2001, 412, 822–826 Pogacic, V.; Bullock, A. N.; Fedorov, O.; Filippakopoulos, P.; Gasser, C.; Biondi, A.; Meyer-Monard, S.; Knapp, S.; Schwaller, J. Structural Analysis Identifies Imidazo[1,2-b]Pyridazines as PIM Kinase Inhibitors with In vitro Antileukemic

6.

Activity. Cancer Res. 2007, 67, 6916–6924. Tsuganezawa, K.; Watanabe, H.; Parker, L.; Yuki, H.; Taruya, S.; Nakagawa, Y.; Kamei, D.; Mori, M.; Ogawa, N.; Tomabechi, Y.; Handa, N.; Honma, T.; Yokoyama, S.; Kojima, H.; Okabe, T.; Nagano, T.; Tanaka, A. A Novel Pim-1 Kinase Inhibitor Targeting Residues That Bind the Substrate Peptide. J. Mol. Biol.

7.

2012, 417, 240–252. Parker, L.; Watanabe, H.; Tsuganezawa, K.; Tomabechi, Y.; Handa, N.; Shirouzu, M.; Yuki, H.; Honma, T.; Ogawa, N.; Nagano, T.; Yokoyama, S.; Tanaka, A. Flexibility of the P-loop of Pim-1 Kinase: Observation of a Novel Conformation 36


Page 37 of 59

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


Induced by Interaction with an Inhibitor. Acta Crystallogr Sect F Struct Biol Cryst 8.

Commun. 2012, 68, 860–866. Nakano, H.; Saito, N.; Parker, L. J.; Tada, Y.; Abe, M.; Tsuganezawa, K., Yokoyama, S.; Tanaka, A.; Kojima, H.; Okabe, T.; Nagano, T. Rational Evolution of a Novel Type of Potent and Selective Proviral Integration Site in Moloney Murine Leukemia Virus Kinase 1 (PIM1) Inhibitor from a Screening-Hit

9.

Compound. J. Med. Chem. 2012, 55, 5151–5164. Stumpfe, D.; Bajorath, J. Exploring Activity Cliffs in Medicinal Chemistry. J. Med.

Chem. 2012, 55, 2932–2942. 10. Hu, Y.; Stumpfe, D.; Bajorath, J. Advancing the Activity Cliff Concept. F1000Research 2013, 2:199. 11. Stumpfe, D.; de la Vega de León, A.; Dimova, D.; Bajorath, J. Follow Up: Advancing the Activity Cliff Concept, Part II. F1000Research 2014, 3:75. 12. Stumpfe, D.; Hu, Y.; Dimova, D.; Bajorath, J. Recent Progress in Understanding Activity Cliffs and Their Utility in Medicinal Chemistry. J. Med. Chem. 2014, 57, 18–28. 13. Pérez-Villanueva, J.; Méndez-Lucio, O.; Soria-Arteche, O.; Medina-Franco, J. L. Activity Cliffs and Activity Cliff Generators Based on Chemotype-related Activity Landscapes. Mol Divers 2015, 19, 1021–1035. 14. Pennington, L. D.; Moustakas, D. T. The Necessary Nitrogen Atom: A Versatile High-Impact Design Element for Multiparameter Optimization. J. Med. Chem. 2017, 60, 3552–3579. 15. Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: an Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701–706. 16. Fedorov, D. G.; Kitaura, K., The Fragment Molecular Orbital Method: Practical Application to Large Molecular System; CRC Press: Boca Raton, FL, 2009. 17. Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electron-Correlated Fragment-Molecular-Orbital Calculations for Biomolecular and Nano Systems. Phys. Chem. Chem. Phys. 2014, 16, 10310–10344. 18. Nakano; T., Kaminuma; T., Sato; T., Akiyama; Y., Uebayasi; M., Kitaura; K. Fragment Molecular Orbital Method: Application to Polypeptides. Chem. Phys. Lett. 2000, 318, 614–618. 19. Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Use of Approximate Electrostatic Potential. Chem. Phys. Lett. 2002, 351, 475–480. 37



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 38 of 59

20. Nakano, T.; Mochizuki, Y.; Yamashita, K.; Watanabe, C.; Fukuzawa, K.; Segawa, K.; Okiyama, Y.; Tsukamoto, T.; Tanaka, S. Development of the Four-body Corrected Fragment Molecular Orbital (FMO4) Method. Chem. Phys. Lett. 2012, 523, 128–133. 21. Watanabe, C.; Fukuzawa, F.; Okiyama, Y.; Tsukamoto, T.; Kato, A.; Tanaka, S.; Mochizuki, Y.; Nakano, T. Three- and Four-body Corrected Fragment Molecular Orbital Calculations with a Novel Subdividing Fragmentation Method Applicable to Structure-based Drug Design. J. Mol. Graphics Modell. 2013, 41, 31–42. 22. Fedorov, D. G.; Kitaura, K. Pair Interaction Energy Decomposition Analysis. J. Comp. Chem. 2007, 28, 222-237. 23. Tsukamoto, T.; Kato, K.; Kato, A.; Nakano, T.; Mochizuki, Y.; Fukuzawa, K. Implementation of Pair Interaction Energy Decomposition Analysis and Its Applications to Protein-Ligand Systems. J. Comput. Chem. Jpn. 2015, 14, 1-9. 24. Fukuzawa, K.; Kitaura, K.; Uebayasi, M.; Nakata, K.; Kaminuma, T.; Nakano, T. Ab Initio Quantum Mechanical Study of the Binding Energies of Human Estrogen Receptor α with Its Ligands: An Application of Fragment Molecular Orbital Method. J. Comput. Chem. 2005, 26, 1-10. 25. Fukuzawa, K.; Mochizuki, Y.; Tanaka, S.; Kitaura, K.; Nakano, T. Molecular Interactions between Estrogen Receptor and Its Ligand Studied by the ab Initio Fragment Molecular Orbital Method. J. Phys. Chem. B 2006, 110, 16102-16110. 26. Watanabe, C.; Fukuzawa, K.; Tanaka, S.; Aida-Hyugaji, S. Charge Clamps of Lysines and Hydrogen Bonds Play Key Roles in the Mechanism to Fix Helix 12 in the Agonist and Antagonist Positions of Estrogen Receptor α: Intramolecular Interactions Studied by the Ab Initio Fragment Molecular Orbital Method. J. Phys. Chem. B 2014, 118, 4993-5008. 27. Ishikawa, T.; Mochizuki, Y.; Amari, S.; Nakano, T. Tokiwa, H.; Tanaka, S.; Tanaka, K. Fragment Interaction Analysis Based on Local MP2. Theor. Chem. Acc. 2007, 118, 937-945. 28. Yoshida, T.; Yamagishi, K.; Chuman, H. QSAR Study of Cyclic Urea Type HIV-1 PR Inhibitors Using Ab Initio MO Calculation of Their Complex Structures with HIV-1. QSAR Comb. Sci. 2008, 27, 694-703. 29. Yoshida, T.; Fujita, T.; Chuman, H. Novel Quantitative Structure-Activity Studies of HIV-1 Protease Inhibitors of the Cyclic Urea Type Using Descriptors Derived from Molecular Dynamics and Molecular Orbital Calculations. Curr. Comp.-Aided Drug Des. 2009, 5, 38-55.

38


Page 39 of 59

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


30. Hitaoka, S.; Harada, M.; Yoshida, T.; Chuman, H. Correlation Analyses on Binding Affinity of Sialic Acid Analogues with Influenza Virus Neuraminidase-1 Using ab Initio MO Calculations on Their Complex Structures. J. Chem. Inf. Model. 2010, 50, 1796-1805. 31. Hitaoka, S.; Matoba, H.; Harada, M.; Yoshida, T.; Tsuji, D.; Hirokawa, T.; Itoh, K.; Chuman, H. Correlation Analyses on Binding Affinity of Sialic Acid Analogues and Anti-Influenza Drugs with Human Neuraminidase Using ab Initio MO Calculations on Their Complex Structures-LERE-QSAR Analysis (IV). J. Chem. Inf. Model. 2011, 51, 2706-2716. 32. Nakanishi, I.; Fedorov, D. G.; Kitaura, K. Molecular Recognition Mechanism of FK506 Binding Protein: An All-electron Fragment Molecular Orbital Study. Proteins: Struct., Funct., Bioinf. 2007, 68, 145-158. 33. Watanabe, H.; Tanaka, S.; Okimoto, N.; Hasegawa, A.; Taiji, M.; Tanida, Y.; Mitsui, T.; Katsuyama, M.; Fujitani H. Comparison of Binding Affinity Evaluations for FKBP Ligands with State-of-the-art Computational Methods: FMO, QM/MM, MM-PB/SA and MP-CAFEE Approaches. Chem-Bio Informatics Journal 2010, 10, 32-45. 34. Wissner, A.; Berger, D. M.; Boschelli, D. H.; Floyd, M. B. Jr.; Greenberger, L. M.; Gruber, B. C.; Johnson, B. D.; Mamuya, N.; Nilakantan, R.; Reich, M. F.; Shen, R.; Tsou, H.-R.; Upeslacis, E.; Wang, Y. F.; Wu, B.; Ye, F.; Zhang, N. 4-Anilino-6,7-dialkoxyquinoline-3-carbonitrile Inhibitors of Epidermal Growth Factor Receptor Kinase and Their Bioisosteric Relationship to the 4-Anilino-6,7-dialkoxyquinazoline Inhibitors. J. Med. Chem. 2000, 43, 3244-3256. 35. Chen, J. M.; Xu, S. L.; Wawrzak, Z.; Basarab, G. S.; Jordan, D. B. Structure-Based Design of Potent Inhibitors of Scytalone Dehydratase: Displacement of a Water Molecule from the Active Site. Biochemistry 1998, 37, 17735-17744. 36. Manta, S.; Xipnitou, A.; Kiritsis, C.; Kantsadi, A. L.; Hayes, J. M.; Skamnaki, V. T.; Lamprakis, C.; Kontou, M.; Zoumpoulakis, P.; Zographos, S. E.; Leonidas, D. D.; Komiotis, D. 3’-Axial CH2OH Substitution on Glucopyranose does not Increase Glycogen Phosphorylase Inhibitory Potency. QM/MM-PBSA Calculations Suggest Why. Chem Biol Drug Des 2012, 79, 663–673. 37. Biela, A.; Khayat, M.; Tan, H.; Kong, J.; Heine, A.; Hangauer, D.; Klebe, G. Impact of Ligand and Protein Desolvation on Ligand Binding to the S1 Pocket of Thrombin. J. Mol. Biol. 2012, 418, 350–366.

39



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 40 of 59

38. Xu, L.; Huiyong, S.; Li, Wang, J.; Hou, T. Assessing the Performance of MM/PBSA and MM/GBSA Methods. 3. The Impact of Force Fields and Ligand Charge Models. J. Phys. Chem. B 2013, 117, 8408–8421. 39. Sun, H.; Li, Y.; Tian, S.; Xu, L.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys. 2014, 16, 16719-16729. 40. Sun, H.; Li, Y.; Shen, M.; Tian, S.; Xu, L.; Pan, P.; Guan, Y.; Hou,T. Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 2014, 16, 22035-22045. 41. Hou, T.; Li, N.; Li, Y.; Wang, W. Characterization of Domain–Peptide Interaction Interface: Prediction of SH3 Domain-Mediated Protein–Protein Interaction Network in Yeast by Generic Structure-Based Models. J. Proteome Res. 2012, 11, 2982−2995. 42. Komeiji, Y.; Ishida, T.; Fedorov, D. G.; Kitaura, K. Change in a Protein’s Electronic Structure Induced by an Explicit Solvent: An Ab Initio Fragment Molecular Orbital Study of Ubiquitin. J. Comput. Chem. 2007, 28, 1750-1762. 43. Fukuzawa, K.; Kurisaki, I.; Watanabe, C.; Okiyama, Y.; Mochizuki, Y.; Tanaka, S.; Komeiji, Y. Explicit Solvation Modulates Intra- and Inter-Molecular Interactions within DNA: Electronic Aspects Revealed by the Ab Initio Fragment Molecular Orbital (FMO) Method. Comp. Theor. Chem. 2015, 1054, 29-37. 44. Fedorov, D. G.; Kitaura, K.; Li, H.; Jensen, J. H.; Gordon, M. S. The Polarizable Continuum Model (PCM) Interfaced with the Fragment Molecular Orbital Method (FMO). J. Comp. Chem. 2006, 27, 976-985. 45. Fedorov, D.G.; Kitaura, K. Energy Decomposition Analysis in Solution Based on the Fragment Molecular Orbital Method. J. Phys. Chem. A 2012, 116, 704–719. 46. Hou, Z.; Nakanishi, I.; Kinoshita, T.; Takei, Y.; Yasue, M.; Misu, R.; Suzuki, Y.; Nakamura, S.; Kure, T.; Ohno, H.; Murata, K.; Kitaura, K.; Hirasawa, A.; Tsujimoto, G.; Oishi, S.; Fujii, N. Structure-Based Design of Novel Potent Protein Kinase CK2 (CK2) Inhibitors with Phenyl-azole Scaffolds. J. Med. Chem. 2012, 55, 2899−2903. 47. Koyama, Y.; Ueno-Noto, K.; Takano, K. Affinity of HIV-1 Antibody 2G12 with Monosaccharides: A Theoretical Study Based on Explicit and Implicit Water Models. Comp. Biol. Chem. 2014, 49, 36-44.

40


Page 41 of 59

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. Watanabe, H.; Okiyama, Y.; Nakano, T.; Tanaka, S. Incorporation of Solvation Effects into the Fragment Molecular Orbital Calculations with the Poisson-Boltzmann Equation. Chem. Phys. Lett. 2010, 500, 116-119. 49. Mazanetz, M. P.; Ichihara, O.; Law, R. J.; Whittaker, M. Prediction of Cyclin-dependent Kinase 2 Inhibitor Potency Using the Fragment Molecular Orbital Method. J. Cheminformatics. 2011, 3, 2. 50. Shigemitsu, Y., Quantum Chemical Study on Molecular-Level Affinity of DJ-1 Binding Compounds, Int. J. Quant. Chem. 2013, 113, 574-579. 51. Wichapong, K.; Rohe, A.; Platzer, C.; Slynko, I.; Erdmann, F.; Schmidt, M.; Sippl, W. Application of Docking and QM/MM-GBSA Rescoring to Screen for Novel Myt1 Kinase Inhibitors. J. Chem. Inf. Model. 2014, 54, 881−893. 52. Lu, H.; Huang, X.; AbdulHameed, M. D. M.; Zhan, C.-G. Binding Free Energies for Nicotine Analogs Inhibiting Cytochrome P450 2A6 by a Combined Use of Molecular Dynamics Simulations and QM/MM-PBSA Calculations. Bioorganic & Medicinal Chemistry 2014, 22, 2149-2156. 53. Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H.-B.; Ding, L.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678-5796. 54. Tsukamoto, T.; Mochizuki, Y.; Watanabe, N.; Fukuzawa, K.; Nakano, T. Partial Geometry Optimization With FMO-MP2 Gradient: Application to TrpCage. Chem. Phys. Lett. 2012, 535, 157-162. 55. Tokuda, K.; Watanabe, C.; Okiyama, Y.; Mochizuki, Y.; Fukuzawa, K.; Komeiji, Y. Hydration of Ligands of Influenza Virus Neuraminidase Studied by the Fragment Molecular Orbital Method. J. Mol. Graphics Modell. 2016, 69, 144-153. 56. Umezawa, Y.; Nishio, M. CH/π Interactions as Demonstrated in the Crystal structure of Guanine-nucleotide Binding Proteins, Src homology-2 Domains and Human Growth Hormone in Complex with their Specific Ligands. Bioorg. Med. Chem. 1998, 6, 493-504. 57. Nishio, M. The CH/π Hydrogen Bond in Chemistry. Conformation, Supramolecules, Optical Resolution and Interactions Involving Carbohydrates. Phys. Chem. Chem. Phys. 2011, 13, 13873-13900. 58. Kigawa, T.; Yabuki, T.; Matsuda, N.; Matsuda, T.; Nakajima, R.; Tanaka, A.; Yokoyama, S. Preparation of Escherichia Coli Cell Extract for Highly Productive Cell-Free Protein Expression. J. Struct. Funct. Genomics 2004, 5, 63–68.

41



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 42 of 59

59. Kigawa, T.; Matsuda, T.; Yabuki, T.; Yokoyama, S. Bacterial Cell-Free System for Highly Efficient Protein Synthesis. In Cell-Free Protein Synthesis Methods and Protocols.; Spirin, A. S.; Swartz, J. R.; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim , Germany, 2008; pp. 83–97. 60. Matsuda, T.; Koshiba, S.; Tochio, N.; Seki, E.; Iwasaki, N.; Yabuki, T.; Inoue, M.; Yokoyama, S.; Kigawa, T. Improving Cell-Free Protein Synthesis for Stable-Isotope Labeling. J. Biomol. NMR 2007, 37, 225–229. 61. Jacobs, M. D.; Black, J.; Futer, O.; Swenson, L.; Hare, B.; Fleming, M.; Saxena, K. Pim-1 Ligand-bound Structures Reveal the Mechanism of Serine/Threonine Kinase Inhibition by LY294002. J. Biol. Chem. 2005, 280, 13728–13734. 62. Leslie, A.G.W. Recent Changes to the MOSFLM Package for Processing Film and Image Plate Data. Jnt CCP4/ESF–EACBM Newsl. Protein Crystallogr. 1992, 26. 63. Winn, M. D.; Ballard, C. C.; Cowtan, K. D.; Dodson, E. J.; Emsley, P.; Evans, P. R.; Keegan, R. M.; Krissinel, E. B.; Leslie, A. G. W.; McCoy, A.; McNicholas, S. J.; Murshudov, G. N.; Pannu, N. S.; Potterton, E. A.; Powell, H. R.; Read, R. J.; Vagin, A.; Wilson, K. S. Overview of the CCP4 Suite and Current Developments. Acta Cryst. 2011, D67, 235–242. 64. Adams, P. D.; Afonine, P. V.; Bunkóczi, G.; Chen, V. B.; Davis, I. W.; Echols, N.; Headd, J. J.; Hung, L.-W.; Kapral, G. J.; Grosse-Kunstleve, R. W.; McCoy, A. J.; Moriarty, N. W.; Oeffner, R.; Read, R. J.; Richardson, D. C.; Richardson, J. S.; Terwilliger, T. C.; Zwart, P. H. PHENIX: A Comprehensive Python-based System for Macromolecular Structure Solution. Acta Cryst. 2010, D66, 213–221. 65. Emsley, P.; Cowtan, K. Coot: Model-building Tools for Molecular Graphics. Acta Cryst. 2004, D60, 2126–2132 66. Wang, J., Cieplak, P., Kollman, P.A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049–1074. 67. Molecular Operating Environment (MOE), 2013.08; Chemical Computing Group Inc., 1010 Sherbooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2013. 68. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; 42


Page 43 of 59

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


Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A.; Peralta, Jr., J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford CT, 2009 69. Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999-3093. 70. Okiyama, Y.; Nakano, T., Yamashita, K.; Mochizuki, Y.; Taguchi, N.; Tanaka, S. Acceleration of Fragment Molecular Orbital Calculations with Cholesky Decomposition Approach. Chem. Phys. Lett. 2010, 490, 84-89. 71. BioStation: ABINIT-MP and BioStation Viewer. The program package is available at: http://www.ciss.iis.u-tokyo.ac.jp/english/dl/ (accessed July 25, 2017) 72. MIZUHO/BioStation Viewer, version 3.0; Mizuho information and research institute Inc: Tokyo, 2013. 73. Srinivasan, J.; Cheatham, T. E.; Cieplak, P.; Kollman, P. A.; Case, D. A. Continuum Solvent Studies of the Stability of DNA, RNA and Phosphoramidate-DNA Helices. J. Am. Chem. Soc. 1998, 120, 9401−9409. 74. Case, D.A.; Darden, T.A.; Cheatham III, T.E.; Simmerling, C.L.; Wang, J.; Duke, R.E.; Luo, R.; Crowley, M.; Walker, R.C.; Zhang, W.; Merz, K.M.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossváry, I.; Wong, K.F.; Paesani, F.; Vanicek, J.; Wu, X.; Brozell, S.R.; Steinbrecher, T.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Mathews, D.H.; Seetin, M.G.; Sagui, C.; Babin, V.; Kollman, P.A. AMBER 10; University of California: San Francisco, 2008. 75. 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 2006, 65, 712-725. 76. Wang, J; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157-1174. 77. Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21, 132-146. 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

Page 44 of 59

78. Jakalian, A; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-quality Atomic Charges. AM1-BCC Model: II. Parameterization and validation. J. Comput. Chem. 2002, 23, 1623-1641. 79. Gogonea, V.; Merz, K. M. Fully Quantum Mechanical Description of Proteins in Solution. Combining Linear Scaling Quantum Mechanical Methodologies with the Poisson–Boltzmann Equation. J. Phys. Chem. A 1999, 103, 5171–5188.

44


Page 45 of 59

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


Scheme caption Scheme 1. Flow of calculations used to determine ligand binding energies with the FMO+MM-PBSA method.

Figures caption Figure 1. (A) Compound structures and experimental IC50 values. Substructures of benzofuranone, indole, and piperazine are shown in broken line circles (i), (ii), and (iii), respectively. (B) Pim1 and hydrogen bond network. (C) CH-π interactions between ligand and Cpd 1. The typical and weak CH-π interactions are indicated by purple and orange lines, respectively, which were analyzed using the CHPI program, as shown in SI Figure S2.

Figure 2. Schematic diagram of the FMO+MM-PBSA method.

Figure 3. Optimization schemes of Pim1 and inhibitor complex with a MM calculation (A) and with a QM/MM calculation (B). Orange-shaded area in (B) was treated as a QM region.

45



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 46 of 59

Figure 4. Fragmentation of amino acid residues (A) and compounds (B). Cpd 1 is shown as an example. Similar fragmentation scheme was also employed for other compounds.

Figure 5. Comparison of correlation between pIC50 and binding energies among three structures: X-ray, MM-opt, and QM/MM-opt. The binding free energies obtained by bind FMO+MM-PBSA, − ∆∆GFMO , are shown in (A), (B), and (C). The binding energies

int , are shown in (D), (E), and (F). The calculated with the FMO method, − ∆∆EFMO

bind binding free energies based on the MM-PBSA, − ∆∆GMM , are shown in (G), (H), and

(I). Root mean squared error (RMSE) of the binding energies for each method is depicted in each graph. The vertical axis is a binding energy in kcal/mol.

Figure 6. Visualization of interaction energies between Pim1 and Fragment (1) of Cpd 1. (A) Interaction energy (IFIE) analysis by the FMO3-MP2 method. For Fragment (1) of Cpd 1 (shown in yellow), main/side chain fragments of Pim1 with attractive and repulsive interactions are represented by red and blue, respectively. (B) Visualization of interaction energies by the PIEDA method. For Fragment (1) of Cpd 1 (shown in

46


Page 47 of 59

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


yellow), the main components of the stabilizing interaction of the main/side chain fragments are represented by color scheme (ES: red and blue; DI: green and white).

Figure 7. Electrostatic potential of Cpd 1 (A) and Cpd 2 (B). Negative electrostatic potentials of ligands in the complex between Pim1 and Cpds 1–6 (C–H).

Figure 8. PIEDA between fragment(1) of ligands and selected amino acid residues of Pim1.

47



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 48 of 59

Tables bind int , and binding free energies, ∆GFMO , Table 1. FMO-based binding energies, ∆E FMO

between Pim1 and ligands for three types of complexes (X-ray, MM-opt, and QM/MM-opt). int (kcal/mol) ∆EFMO

Fragment (1)

Fragment (2)

bind (kcal/mol) ∆GFMO

Fragment (1) + (2)

Fragment (1) + (2)

Structure

Cpd

IC50 (nM)

HF

MP2

HF

MP2

HF

MP2

HF

MP2

X-ray

1

2

–22.69

–77.53

–371.39

–388.69

–394.07

–466.22

5.12

–67.03

3

3

–41.45

–97.88

–375.93

–392.44

–417.37

–490.32

–12.15

–85.10

5

92

–25.89

–72.99

–376.48

–391.31

–402.37

–464.30

2.75

–59.18

6

447

–24.36

–78.58

–367.24

–383.73

–391.60

–462.32

9.89

–60.83

R with –pIC50

0.41

0.51

0.42

0.61

0.46

0.61

0.61

0.70

R2 with –pIC50

0.17

0.26

0.17

0.37

0.21

0.37

0.37

0.49

Fragment (1)

Fragment (2)

Fragment (1) + (2)

Fragment (1) + (2)

Structure

Cpd

IC50 (nM)

HF

MP2

HF

MP2

HF

MP2

MM-opt

1

2

–23.58

–78.07

–370.13

–387.37

–393.71

–465.44

4.78

–66.95

2

2

–28.70

–81.65

–371.84

–389.15

–400.53

–470.80

–3.64

–73.91

3

3

–34.70

–88.44

–371.22

–388.50

–405.92

–476.95

2.08

–68.95

4

53

–21.33

–74.84

–369.57

–386.78

–390.90

–461.63

2.67

–68.06

5

92

–22.77

–76.42

–371.21

–388.50

–393.98

–464.92

6.80

–64.14

6

447

–24.46

–77.59

–366.31

–383.38

–390.77

–460.98

2.12

–68.09

R with –pIC50

0.54

0.59

0.72

0.72

0.69

0.73

0.37

0.49

R2 with –pIC50

0.30

0.35

0.52

0.52

0.48

0.53

0.14

0.24

Fragment (1)

Fragment (2)

Fragment (1) + (2)

Structure

Cpd

IC50 (nM)

HF

MP2

HF

MP2

HF

MP2

QM/MM-opt

1

2

–24.64

–77.75

–370.40

–387.66

–395.04

2

2

–28.70

–80.98

–371.64

–388.95

3

3

–31.52

–83.83

–371.83

4

53

–21.51

–73.70

5

92

–23.40

6

447

R with –pIC50 2

R with –pIC50

HF

MP2

Fragment (1) + (2) HF

MP2

–465.41

–2.48

–72.85

–400.34

–469.93

–7.27

–76.86

–389.13

–403.34

–472.96

–6.78

–76.40

–369.47

–386.68

–390.98

–460.38

–1.09

–70.49

–75.50

–371.41

–388.71

–394.81

–464.21

–0.06

–69.46

–24.14

–76.08

–366.12

–383.19

–390.26

–459.27

1.69

–67.32

0.65

0.71

0.73

0.73

0.78

0.82

0.87

0.92

0.42

0.51

0.54

0.54

0.60

0.67

0.76

0.85

48


Page 49 of 59

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


Table 2.

Interaction energies between Cpds 1–6 and Pim1. The FMO3-IFIEs are

included in both the MP2 total energies and the electron correlation energies (CORR). For FMO2-PIEDA, each energy component is shown. The CORR is assumed to be identical to the DI term in this paper. int of Fragment (1) (kcal/mol) ∆E FMO

FMO3-IFIE Cpd

IC50 (nM)

MP2

1

2

–77.75

2

2

3

FMO2-PIEDA

CORR

Total

ES

–53.11

–75.91

–60.94

–17.18

–52.23

–80.98

–52.28

–79.86

–63.34

–18.11

–51.60

3

–83.83

–52.31

–81.73

–63.47

–16.44

–51.50

4

53

–73.70

–52.19

–72.02

–57.15

–17.06

–51.40

5

92

–75.50

–52.10

–73.61

–59.28

–16.94

–51.26

6

447

–76.08

–51.94

–74.32

–58.62

–16.74

–51.12

R with –pIC50

0.71

0.71

0.73

0.81

0.40

0.79

R2 with –pIC50

0.51

0.50

0.53

0.65

0.16

0.62

∆E

int FMO

IC50 (nM)

MP2

1

2

–465.41

2

2

3

DI

of Fragment (1) + (2) (kcal/mol)

FMO3-IFIE Cpd

CT+mix

FMO2-PIEDA Total

ES

–70.36

–466.00

–447.75

–36.86

–69.51

–469.93

–69.59

–471.24

–451.30

–37.83

–68.93

3

–472.96

–69.62

–473.31

–451.58

–36.17

–68.83

4

53

–460.38

–69.40

–461.12

–443.08

–36.62

–68.63

5

92

–464.21

–69.40

–464.73

–446.92

–36.67

–68.58

6

447

–459.27

–69.01

–460.01

–441.56

–36.12

–68.23

0.82

0.79

0.83

0.86

0.54

0.86

0.67

0.63

0.68

0.74

0.29

0.74

R with –pIC50 2

R with –pIC50

CORR

CT+mix

DI

49



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 50 of 59

Page 51 of 59


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



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


Page 52 of 59

Page 53 of 59

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




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 54 of 59

Page 55 of 59


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



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


Page 56 of 59

Page 57 of 59


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



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


Page 58 of 59

Page 59 of 59

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19



Theoretical Analysis of Activity Cliffs among Benzofuranone Class

Recommend Documents