A molecular interaction field describing nonconventional intermolecular interactions and its application to protein–ligand interaction prediction

Journal of Molecular Graphics and Modelling 96 (2020) 107515

Contents lists available at ScienceDirect

Journal of Molecular Graphics and Modelling journal homepage: www.elsevier.com/locate/JMGM

A molecular interaction field describing nonconventional intermolecular interactions and its application to proteineligand interaction prediction Daichi Hayakawa*, Nahoko Sawada, Yurie Watanabe, Hiroaki Gouda** School of Pharmacy, Showa University, 1-5-8 Hatanodai, Shinagawa-ku, Tokyo, 142-8555, Japan

a r t i c l e i n f o

a b s t r a c t

Article history: Received 19 August 2019 Received in revised form 18 December 2019 Accepted 18 December 2019 Available online 23 December 2019

Nonconventional noncovalent interactions such as CH/O, CH/p, and halogen bonds play important roles in molecular recognition in biological systems and have been increasingly exploited in structure-based drug design. In silico approaches that consider these interactions would be an effective strategy for drug discovery. The computation of the molecular interaction field (MIF), which is a three-dimensional (3D) potential map describing the interactions formed around a target compound, would assist the design of molecules that bind to the target biomolecule via nonconventional interactions. In this study, we developed a novel MIF calculation method that describes nonconventional interactions. This method evaluates the MIF as the interaction energy between the target ligand molecule and probe molecule. To describe the nonconventional interactions, our method employs quantum chemical calculations with four types of probe molecules. The calculated MIFs for casein kinase 2 (CK2) inhibitors correctly identify the halogen bond, CH/p, and CH/O interactions formed in the CK2/inhibitor complexes. Additionally, we have developed a method for calculating the proteineligand interaction energy (Eint) based on the MIF and a coarse-grained protein model. The calculated interaction energies for CK2 inhibitors correlate with the experimental log(Ki) values. Thus, MIF and Eint obtained by our method show promise as descriptors for protein-ligand interaction prediction by considering nonconventional noncovalent interactions. © 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Keywords: Molecular interaction field Nonconventional hydrogen bonds Halogen bonds Proteineligand binding Structure-based drug design

1. Introduction The intermolecular interactions between a drug candidate and a target biopolymer are the main consideration in structure-based drug design (SBDD) [1,2]. In particular, classical hydrogen bonds and ionic bonds are usually the dominant factors in ligand binding to biopolymers. These interactions are, thus, the main targets in SBDD. However, in addition to these interactions, nonconventional noncovalent interactions such as halogen bonds, CH/n (n ¼ O, N, F, Cl, etc.), CH/p, and XH/p (X ¼ O, N, etc.) interactions are also suggested to play an important role in molecular recognition in biological systems (Fig. 1) [3,4]. A halogen bond is a noncovalent bond between the s-hole on a halogen atom in a molecular entity and an electron-rich region in

* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (D. Hayakawa), godah@ pharm.showa-u.ac.jp (H. Gouda).

another, or the same, molecular entity [4]. In fact, interactions between iodine and aromatic hydrocarbons were suggested based on spectroscopic studies as early as the 1940s and 1950s [5,6]. Later, intermolecular halogen bonds between oxygen and bromine in bromine 1,4-dioxanate crystals were explicitly observed by X-ray crystallography experiments in the 1950s by Hassel and Hvoslef and Hassel and Strømme [7,8]. XH/p, CH/n, and CH/p interactions are classified as weak hydrogen bonds. XH/p interactions, such as OH/p and NH/p, are weak hydrogen bonds between hydroxyl or NH groups and p-electrons. This type of interaction was studied intensively by infrared spectroscopy in the 1950s and 1960s [3]. Additionally, CH/n interactions (n ¼ O, N, F, Cl, etc.) are another type of weak hydrogen bond between CH groups such as alkyl group and lone pair electrons [3]. Sutor first suggested the CH/O interaction, which was the first example of a CH/n-type hydrogen bond, based on crystallographic data in 1962 [9,10]. CH/p interactions, which were first suggested by Nishio et al., in 1977, are weak hydrogen bonds between CH groups and p electrons [11e14]. The evidence for the existence of halogen bonds [15], XH/p [16], CH/n [17e20], and CH/p [16,21] interactions in small molecules has been obtained

https://doi.org/10.1016/j.jmgm.2019.107515 1093-3263/© 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

2

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

Fig. 1. Examples of CH/p interactions, halogen bonds, and CH/O interactions formed in a proteineligand complex. This figure is based on the structure of the casein kinase 2 (CK2)/ K25 complex (PDB ID: 1ZOE).

by statistical analyses using the Cambridge Structural Database (CSD) [22,23]. Further, analysis of the Protein Data Bank (PDB) (www.rcsb.org) [24] has revealed the formation of halogen bonds [25e27], CH/O [28e30], XH/p [30,31], and CH/p [30,32e34] interactions in biopolymers such as proteins and nucleic acids. Recently, these interactions in biopolymers have been investigated in detail using computational chemistry [3,35e43]. Because nonconventional interactions can be formed between a drug molecule and a biomolecule, they are becoming increasingly important in drug discovery. For instance, the replacement of a proton of the Phosphodiesterase type 5(PDE5) inhibitor with bromine leads to the formation of a halogen bond between PDE5 and the bromine-substituted inhibitor [44]. The formation of the halogen bond results in the four-times-higher inhibitory activity than that of the unmodified PDE5 inhibitor [44]. Additionally, recent quantum mechanics (QM) level calculations have revealed the importance of CH/p interactions for ligand binding to proteins [35e40]. Considering these precedents, nonconventional interactions are promising tools for the expansion of the structural diversity of the drug candidates in SBDD. In the lead optimization stage of SBDD, for example, polar functional groups such as hydroxyl group can be introduced into the lead compound to form hydrogen bonds with the target. However, the introduction of these functional groups is sometimes unfavorable for drug absorption because of their negative influence on drug lipophilicity. In this case, a nonpolar functional group that forms nonconventional interactions could be an alternative to the polar functional group. Therefore, the consideration of the nonconventional interactions could extend the molecular design strategy in drug discovery. A molecular interaction field (MIF) describing nonconventional interactions would be useful for such a molecular design. A MIF is a three-dimensional (3D) potential map describing the interactions formed around a compound of interest [45e47]. It is evaluated as the interaction energy between a compound molecule and a

chemical probe (a molecule or an atom) set in a 3D grid around the compound [45e47]. The probe molecules reflect physicochemical characteristics of a binding partner of the target compound [46]. MIF is widely used as a descriptor to construct quantitative structureeactivity relationship (QSAR) models [48,49] and has been applied to the prediction of the ligand binding sites of proteins [50,51]. MIF is, thus, a very important concept for in silico drug design. The MIF is usually calculated using the standard Coulomb and 6e12 Lennard-Jones potentials. To obtain detailed descriptions of physical properties, various potential functions for MIF have been proposed, including the (1) molecular lipophilicity potential [52e57], (2) polarizability potential [58], and (3) hydrogen bond potential [59]. However, there is no computational method to calculate MIFs focused on nonconventional interactions. The development of MIFs describing nonconventional interactions is, thus, desirable. The main aim of this study is to develop an effective sampling method for MIF that considers nonconventional interactions as well. Herein, we propose a novel method for making this consideration. QM treatment is generally required to evaluate nonconventional interactions, because, for example, the strength of halogen bonds correlates with the s-hole volume on the halogen atom and dispersion forces are dominant in CH/p interactions [3,4,60e63]. Therefore, our method evaluates the MIF at the QM level of the interaction energy between the target compound and the probe molecule. The use of QM calculations also makes it possible to explore the interaction area of the nonconventional interactions [64e66]. In addition, our method employs a set of four small molecules, H2O, NH3, CH4, and benzene, as probe molecules, which are considered to be representative molecules possessing hydrogen-bond donor, hydrogen-bond acceptor, lipophilic, and aromatic properties, respectively. Utilizing these four types of probes makes it possible to analyze the various types of nonconventional interactions.

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

In the present study, MIFs for casein kinase 2 (CK2) inhibitors were first calculated using our method. The obtained MIFs adequately identified the hydrogen bonds, halogen bonds, and CH/ p and CH/O interactions observed in the crystal structures of CK2/ inhibitor complexes. In addition, we developed a method for calculating protein-ligand interaction energies based on an MIF and coarse-grained protein model, for ranking the affinity of multiple ligands for a protein. Its availability for SBDD was investigated based on the proof-of-concept calculation for CK2 inhibitors. The calculated interaction energies for CK2/inhibitor complexes correlated well with the experimental inhibition constants of the inhibitors. These results also support the availability of MIF and Eint obtained by our method as descriptors for the prediction of proteinligand binding affinities taking into consideration the nonconventional noncovalent interactions. 2. Methods 2.1. Generation of MIF using probe molecules and QM calculations The generation of a MIF requires knowledge of the interaction energy between a target ligand molecule and a chemical probe (molecule or atom) at 3D grid points defined around the target [46,47]. Our method evaluates the interaction energies using QM calculations with four types of probe molecules: H2O, NH3, CH4, and benzene. The reason for using these four probe molecules is described in the next section. The calculation procedure is shown in Fig. 2. First we must prepare active conformations (protein-bound conformations) of the target ligand molecules. These can be obtained from proteineligand complex structures obtained by experimental methods such as X-ray crystallography or by computational methods such as docking calculations. In this work, MIF calculations were performed for the ligand conformations extracted from the crystal structures of proteineligand complexes, which are available from the PDB. The calculation method proposed here can be applied to both neutral and charged ligands. Next, 3D grids are generated around the target ligand molecule (Fig. 2a). A total of 9261 grid points (21  21  21) with an interval of 1 Å were generated in this work. The probe molecule is set on a grid point, as shown in Fig. 2b. Vectors Gn are then defined from the grid point to the atoms constituting the ligand molecule (Fig. 2b). Additionally, a dummy atom(s) is/are defined at the center of the aromatic ring(s) if the ligand molecule contains aromatic ring(s), and a vector(s) Gn is/are also defined from the grid point to the dummy atom(s) (Fig. 2b). The number of vectors (N) is equal to the sum of numbers of atoms constituting the target ligand molecule and dummy atoms (Fig. S1a). Vectors Gn (n ¼ 1, 2, 3, …, N) are numbered in order of their proximity to the grid point (Fig. S1b). Following the probe arrangement process, the probe molecule is oriented such that the unit vector M defined on the probe is coincident with vector Gn (Fig. 2c). The unit vector M is defined from the center of mass of the probe molecule along its molecular axis. The definition of vector M will be described in the following section in detail. The interaction energy between the target ligand molecule and the probe molecule is then calculated for the obtained coordinate (Fig. 2d). The interaction energy is evaluated by single point calculation at the uB97XD/6-31G(d,p) level with basis set superposition error (BSSE) correction [67,68]. Because MIF calculation requires a lot of calculations of intermolecular interactions, the DFT calculation method is desirable. It is known that the uB97XD density functional describes intermolecular interactions appropriately [69,70]. Additionally, the result of the calculation using uB97XD is certainly in good agreement with the results of MP2 and CCSD(T) in the test calculations using model systems (Fig. S2). In this work, uB97XD is thus applied for MIF

3

calculations. The definition of N Gn vectors generates N coordinate sets that differ in the orientation of the probe. The interaction energies are calculated for these coordinate sets and the most stable is selected as the MIF energy of the grid point. However, interactions with closer atoms are generally dominant. The interaction energy for the orientations generated by Gn vectors defined between the grid point and the remote atoms, thus, affect the obtained MIF feature negligibly. In other words, it is not necessary to perform the calculations for the orientations generated by Gn defined between the grid point and the remote atom. The optimum number of vector Gn (¼ Ncut) depends on the structure of the target ligand molecule. In the present work, the optimum Ncut values were determined based on the pre-calculations of the MIF for K25 (a CK2 inhibitor) (Figs. S3 and S4). As a result, the calculation using an H2O probe considering one orientation (Ncut ¼ 1) results in an incomplete MIF, whereas the calculation considering three orientations (Ncut ¼ 3) completely describes the interaction-forming area around K25. The calculations using NH3, CH4, and benzene probes considering one orientation (Ncut ¼ 1) led to MIFs with an interaction-forming area of sufficient resolution. Therefore, Ncut was set at 3 for the H2O probe, whereas Ncut was set at 1 for the others in the MIF calculations demonstrated in this work. Additionally, the calculation was not performed for the orientations generated by Gn defined for atoms 4.0 Å or more apart from the grid point. The full MIF is obtained by repeating the above procedures (Fig. 2e). Intermolecular interactions usually diverge at the short range and converge to zero at the long range (Fig. S2). Therefore, the energy calculation was not performed for grid points within 1.5 Å or further than 5.0 Å from the target ligand molecule to reduce computational cost. The MIF energy of the grid point where the calculation was not performed was set at 0.0 kcal/mol. As a result, 1500 to 3500 energy calculations were performed per target ligand in this work. All QM calculations were performed in Gaussian09 [71]. The arrangement and the orientation of the molecules for MIF calculation were obtained using in-house code. The in-house programs used in this work are currently not available to the public because these are not made expecting utilization by third parties. 2.2. Probe molecules for MIF calculations This section explains our reasons for selecting four probe molecules (H2O, NH3, CH4, and benzene). The probe molecules reflect the physico-chemical properties of chemical groups constituting target bio-macromolecules. The chemical groups constituting proteins are generally classified as hydrogen bond donors (HD), hydrogen bond acceptors (HA), lipophilic (or hydrophobic) groups (LP), and aromatic groups (AR). This classification rule is known to be appropriate for determining the similarity of protein surfaces and ligand binding pockets [72,73]. Therefore, H2O, NH3, CH4, and benzene were selected as probe molecules in this work because these are representative of HD, HA, LP, and AR molecules, respectively. Each probe molecule possesses a unit vector M. The definition of unit vector M of the probe affects the calculated MIFs. The definition of M and the evaluable interactions are described here. H2O and NH3 can be both a hydrogen bond donor and acceptor. However, in this work, the unit vector M of NH3 was defined along its C3 axis, as shown in Table 1. Consequently, its lone pair is consistently oriented toward the target ligand molecule at the grid points during the calculation process (Fig. 2d). Therefore, the MIF obtained using the NH3 probe describes the interaction areas formed by the hydrogen bond acceptor groups. On the other hand, unit vector M of the H2O probe was defined along the OeH axis of H2O (Table 1). Thus, the MIF obtained with the H2O probe describes the interaction areas formed by the hydrogen bond donor groups. Similarly,

4

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

Fig. 2. MIF calculation procedure. The definition of grid points around the target ligand molecule (a). The position of a probe molecule at a grid point and the definition of the vector Gn from the grid point to the atom constituting the molecule or the dummy atom (b). A dummy atom defined in the center of the ring structure is shown in orange dots. Orientation of the probe molecule (c). The obtained structure for interaction energy calculation (d). Example of the obtained MIF (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

unit vectors M of CH4 and benzene were defined along the C3 axis of CH4 and the C6 axis of benzene, respectively (Table 1). The protons of CH4 and the p-electrons of benzene are, thus, consistently oriented toward the target ligand molecule during the MIF calculation process. The MIFs with CH4 and benzene probes describe the interaction areas formed by lipophilic and aromatic groups, respectively. For example, the MIF obtained with the CH4 probe would reveal the CH/p interaction-forming area around the aromatic groups of the target ligand molecule. In the case of the benzeneeCH4 complex, the structure where one hydrogen atom of CH4 point towards the center of benzene is the most stable [60]. Therefore, dummy atoms are placed at the center of the aromatic rings of the target ligand molecule to treat CH/p interactions adequately. The MIF obtained with the NH3 probe describes the CH/N interactionforming areas around the alkyl groups of the target ligand. However we used NH3 as a hydrogen-bond acceptor probe molecule. Therefore, the interactions around the alkyl groups obtained by NH3 probe are generally regarded as “CH/n (n ¼ O, N, F, Cl, etc.)” interaction-forming areas. Because halogen bonds are noncovalent interactions between the s-hole and Lewis acid, the halogen-bondforming areas around the halogen atoms of the target ligand molecule can also be explored by the NH3 probe [4]. The evaluable interactions of the four probes are listed in Table 1. In this work, the

MIFs obtained using the H2O, NH3, CH4, and benzene probes are denoted MIF(HD), MIF(HA), MIF(LP), and MIF(AR), respectively.

2.3. Computation of the proteineligand interaction energies based on MIFs and coarse-grained protein models Applying the MIF to SBDD requires the evaluation of the proteineligand interaction energies based on the MIFs. Thus, we developed a method of calculating the proteineligand interaction energy based on the calculated MIF and coarse-grained protein model. The outline of the calculation is shown in Fig. 3. This calculation utilizes the 3D structure of the proteineligand complex and the calculated MIF for the active conformation of ligand. In this work, the structures of the complexes obtained from PDB were utilized to avoid errors that arise from the failure of ligand docking. First, the complex is divided into the ligand and receptor that consists of protein, water, and ions (Fig. 3a). The MIF calculation is performed for the ligand molecule as described above. To simplify the protein structure, this method represents the chemical groups of proteins as five types of spheres based on their physicochemical properties: hydrogen-bond donors (HD), hydrogen-bond acceptors (HA), hydrogen-bond donors/acceptors (DA), lipophilic (or hydrophobic) groups (LP), and aromatic groups (AR) (Fig. 3a) [72,73]. To prevent confusion, HD, HA, DA, LP, and AR

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

5

Table 1 The list of the probe molecules. Probe moleculea

Physicochemical propertiesb

Evaluable interactionsc

H2O

H-bond donor (HD)

H-bonds, XH/p

NH3

H-bond acceptor (HA)

H-bonds, CH/n, halogen bonds

CH4

Lipophilic (LP)

CH/n, CH/p

Benzene

Aromatic (AR)

p/p, CH/p, Halogen/p

a b c

Structures of probe molecules. The molecular axes defined on these probe molecules are representexd by arrows. Classes of the probes in physicochemical property. Example of the evaluable interactions by these probes.

Fig. 3. Outline of the calculation of proteineligand interaction energy (a). (i) The ligand with a binding conformation is extracted from the complex. (ii) A MIF calculation is performed for the ligand. (iii) The protein structure is coarse grained using five types of spheres [72,73]. (iv) The calculated MIF is superimposed onto the coarse-grained-proteineligand complex. The enlarged view of the binding site (b). The Rec(HA) spheres (blue spheres) indicated by orange arrows are in the interaction area of the calculated MIF(HA). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

spheres defined on protein structure are denoted as Rec(HD), Rec(HA), Rec(DA), Rec(LP), and Rec(AR) spheres, respectively. An example of the allocation of spheres in the coarse-grained model are shown in Fig. 4 and the details for all amino acid residues are shown in Figs. S5eS8. One property sphere is usually allocated for each chemical group. As shown in Fig. 4, the property sphere is usually placed at the center of the heavy atom, which is highlighted. In the case of aromatic ring groups, the sphere is arranged at the center of the ring. Following the construction of the

coarse-grained model, the calculated MIF is overlaid on the protein model (Fig. 3a). Our method regards a grid point of the MIF that is closest to the coordinate of similar property sphere defined on protein as the interaction point between the protein and the ligand. For instance, Fig. 3b shows the superposition of the MIF(HA) around the active conformation of K25, a CK2 inhibitor, (MIF energy < 2 kcal/mol) and the property spheres defined on the ligand binding site in CK2. There are several Rec(HA) spheres within the MIF(HA). Therefore, we can consider that the chemical

6

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

2.4. The protein and its inhibitors utilized as the test set for MIF calculations

Fig. 4. Example of the coarse-grained amino acid residue obtained using five types of property spheres: hydrogen bond donor (HD), hydrogen bond acceptor (HA), hydrogen bond donor/acceptor (DA), lipophilic (LP), and aromatic (AR) [72,73].

groups corresponding to these Rec(HA) spheres in CK2 (that is, hydrogen-bond acceptor groups in CK2) would form favorable interactions with K25 having energies of at least 2 kcal/mol. More precise interaction energies could be estimated from the MIF energies of the grid points that are extremely close to the coordinates of these Rec(HA) spheres defined on CK2. Therefore, we consider these MIF grid points as the interaction points of the K25-CK2 complex, and we assume that the MIF energies of these grid points approximately represent the intermolecular interaction energies between the hydrogen bond acceptor groups in CK2 and K25. Based on this assumption, the interaction energy between the ligand and the property spheres defined on protein can be expressed by eq (1).

ES;i ¼

X   ES ðjÞd RS;i  Rj ; S ¼ HD; HA; DA; LP; AR

(1)

j

Here, Es(j) is the MIF energy of probe S at grid point j. Rs,i and Rj are the position vectors of the ith property sphere S defined on protein and grid point j. d is the delta function, which is set at d ¼ 1 for grid point j within 1.0 Å from ith property sphere S and d ¼ 0 for the others. d is set at d ¼ 1 for the grid point with the most stable MIF energy if there are some grid points within 1.0 Å. For Rec(DA) spheres, their ES,i values are evaluated using both MIF(HD) and MIF(HA), and the energetically most stable value is defined as EDA,i. Finally, the total proteineligand interaction energy is given by eq (2).

X X X X X Eint ¼ EHD;i þ EHA;i þ EDA;i þ ELP;i þ EAR;i i

i

i

i

(2)

i

Each term of eq (2) corresponds to the HD, HA, DA, LP, and AR parts of the total interaction energy. This method represents the proteineligand interaction energy as the addition of the interaction energies between the ligand molecules and chemical groups of the protein; this corresponds to the two-body approximation. That is, in this calculation, we mainly have two approximations, namely pairwise additive and coarse-graining [74]. The allocation of the property spheres was achieved using an in-house code. This code was individually remade with slight modification on the basis of previous reports [72,73]. Additionally, the calculation of the interaction energy based on eqs (1) and (2) was performed using inhouse code. We specify that the in-house programs are currently not open to the public.

MIF calculations were performed for seven casein kinase 2 (CK2) inhibitors that contain halogen atoms whose inhibitory activities and complex structures with CK2 are known [75e78]. The chemical structures, PDB IDs, and the experimentally determined inhibition constant (Ki) values for the seven inhibitors are shown in Table 2 [75e78]. There are CH/O, CH/p, and halogen bond interactions in the CK2/inhibitor complexes (Fig. 1). Therefore, these complex structures are suitable as the test sets for the verification of the interaction-forming areas obtained by the calculation of the MIF. Additionally, the correlation between the interaction energies calculated based on MIF and the Ki values can be evaluated. In this work, the test set for assessment of MIF calculation is required to be a series of inhibitors with known inhibitory activities for the same protein. Additionally, the inhibitors for the test set are required to form intermolecular halogen bonds with the target protein. However, there are few inhibitor sets that satisfy all the above conditions. The CK2/inhibitor complex is one of the few protein/inhibitor systems that can be used for the test set for MIF calculation assessment. 2.5. Preparation of CK2/inhibitor complex structures and probe molecules Because the positions of protons were not determined by X-ray crystallography, protons were added to the proteins, inhibitors, and water molecules of the X-ray crystal structures using the Protein €dinger, LLC) [79]. The protons were Preparation Wizard (Schro energetically minimized using Prime with the OPLS3 force field [80e84]. The heavy atoms were fixed during the minimization. The MIF calculations were performed for the inhibitors extracted from the prepared complex structures. The protonation states of inhibitors were determined using Epik [85,86]. Hydrogen bond networks were analyzed using Protein Preparation Wizard program for the water molecules and the ions within 5 Å of the ligand. Crystal structures of CK2/inhibitor complexes revealed the existence of water molecules and ions that mediate the proteins and the inhibitor, which are defined as the bridging/interstitial waters and ions. Structural optimizations at the uB97XD/6-31G(d,p) level of theory were performed for the probe molecules. The vibrational frequencies were calculated at the same level to confirm that they were stationary points on the potential energy surface. Molecular €dinger Maestro and Visual Mographics were prepared in Schro lecular Dynamics (VMD) [87]. 3. Results and discussion 3.1. Visually verifying the MIF overlaid on the CK2/inhibitor complex The results in this section demonstrate the usefulness of MIFs obtained for CK2 inhibitors. Although MIFs cannot be experimentally observed, they can be indirectly verified by comparison between the interaction area obtained by MIF calculation and the interactions observed in the CK2/inhibitor complexes obtained by X-ray crystallography. The MIFs of the K25 inhibitor overlaid on the CK2/K25 complex are shown in Figs. 5 and 6 as examples. Additionally, the interactions observed in the crystal structure of the CK2/K25 complex and their MIF energies, which correspond to the Es,i values defined by eq (1), are listed in Table 3. Details of the geometrical parameters of the interactions observed in CK2/K25 complex are listed in

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515 Table 2 The list of the CK2 inhibitors applied to the MIF calculation. PDB ID

Ki [nM]d

K25

1ZOEa

40

2

K37

1ZOGa

70

3

K44

1ZOHa

100

4

K32

2OXDb

150

5

K22

2OXXb

200

6

K17

2OXYb

300

7

TBS

1J91c

400

Entry

Compound

1

a b c d

Structure

Structure reported in reference [75]. Structure reported in reference [76]. Structure reported in reference [77]. Data from reference [78].

Table S1. The MIF(HD) for K25 is shown in Fig. 5a. An attractive interaction-forming area with MIF energy of 3 to 6.6 kcal/mol is formed around the unprotonated nitrogen(N3) of the benzimidazole ring of K25 (Fig. 5a).

7

This corresponds to the hydrogen-bond-forming area with N3. The interaction energy of classical hydrogen bonds is typically 3 to 7 kcal/mol [3]. Therefore, the MIF energy of the interaction area is consistent with this value. When we examined the crystal structure of the CK2/K25 complex, a water molecule, w4, was found to be located within this MIF. This water molecule actually forms a hydrogen bond with the N3 group of K25. The coincidence between the calculated interaction area and the interaction site determined by X-ray crystallography suggests the usefulness of the MIF calculation. Additionally, the interaction area with MIF energies of 2 to 3 kcal/mol above and below the benzimidazole plane corresponds to the XH/p interaction-forming area. XH/p interactions are not formed between CK2 and K25 because CK2 has no chemical group that has Rec(HD) properties either above or below the benzimidazolium plane. However, this MIF energy is consistent with the general value of the XH/p interaction reported in previous work [3]. The MIF(HA) around K25 is shown in Fig. 5b. The interactionforming area with MIF energies of 4 to 9.6 kcal/mol around the eNH group of inhibitor K25 corresponds to the hydrogen bond forming area. Within this area, we found a chloride ion Cl derived from the crystal structure of CK2/K25 complex. An ionic hydrogen bond is formed between the NH group of K25 and this Cl ion (Table 3) [75]. These results suggest that the MIF(HA) correctly describes the hydrogen-bond-forming area around the NH group. The interaction area with MIF energy of 2 to 4.2 kcal/mol is formed around the eCH3 groups of inhibitor K25 (Fig. 5b). This corresponds to the CH/n interaction-forming area. The MIF energy of this area is consistent with typical CH/n interaction energies (2 to 4 kcal/mol) [3]. The crystal structure of the CK2/K25 complex suggests the formation of the CH/O interactions between the methyl groups of inhibitor K25 and the carbonyl oxygen of the Arg47 residue and carboxylic oxygen of the Asp175 residue (Table 3) [75]. These carbonyl and carboxyl oxygen atoms are in the CH/n interaction-forming area. These observations highlight the ability of MIF to identify the CH/n interaction-forming area correctly. Interaction-forming areas with interaction energies of 2 to 3 kcal/mol are formed in the direction along the CeBr bond and correspond to the halogen-bond-forming area. Additionally, the halogen bond forming area is limited to a narrow area, which reflects the directionality of the halogen bonds [4,64]. The MIF energies are consistent with the values of the interaction energies of Br/O and Br/N halogen bond (2 to 4 kcal/mol) proposed in the previous study using bromobenzeneeacetone and halomethaneeammonia models [62,88]. A previous X-ray crystallographic analysis revealed the formation of a halogen bond between the ring bromine (5Br) and the carbonyl oxygen of Val116 (Table 3) [75]. As shown in Fig. 5b, the halogen bond forming area is around this carboxyl oxygen. This result again demonstrates the relevance of the interaction-forming areas obtained by MIF calculations. The MIF(LP) around K25 is shown in Fig. 6a. The interactionforming areas are formed in the direction perpendicular to the benzimidazole plane, similar to the MIF(HD). These correspond to the CH/p interaction-forming areas. Its MIF energy is about 1 kcal/ mol weaker than that of the MIF(HD) describing the XH/p interactions and is consistent with typical CH/p interaction energies (1.5 to 2.5 kcal/mol) [3]. In the crystal structure of the CK2/K25 complex, CH/p interactions are formed between the benzimidazole part of K25 and the alkyl groups of the Val53, Ile66, and Ile174 (Table 3). CH/p interaction-forming areas formed around these alkyl groups (Fig. 6a), and this suggests that the calculated MIF validly identifies the CH/p interaction-forming areas. The MIF(AR) around K25 is shown in Fig. 6b. The interactionforming areas with interaction energies of 3 to 7.9 kcal/mol

8

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

Fig. 5. The calculated MIF(HD) (a) and MIF(HA) (b) of compound K25. These MIFs are overlaid onto the CK2/K25 complex structure. The attractive regions are shown as a heat map (surface model).

Fig. 6. Calculated MIF(LP) (a), and MIF(AR) (b) of compound K25. These MIFs are overlaid onto the CK2/K25 complex structure. The attractive regions are shown as a heat map (surface model).

are formed in the direction perpendicular to the benzimidazole plane of inhibitor K25. This represents the p/p interaction-forming area. The MIF energy is consistent with typical values for p/p interactions [89]. Additionally, the interaction-forming areas with interaction energies of 2 to 3 kcal/mol are formed in the direction along the CeBr bond and correspond to the halogen/p

interaction-forming areas. The range of energies is consistent with halogen/p interaction energy between benzene and hydrogen bromide reported previously (2.05 kcal/mol) [90]. Previous X-ray crystallographic studies have revealed that halogen/p interactions are formed between the ring bromine of K25 and the phenyl group of Phe113 in the CK2/K25 complex [75]. The halogen/p interaction-

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515 Table 3 Interactions formed in CK2/K25 complex and the corresponding MIF energies. Donor

Acceptor

Da

MIF energyb

Interaction typec

OH (w4) 1NH (K25) CH3 (K25) CH3 (K25) 5Br (K25) Cg2 (Il66) Cg1 (Val53) Cg2 (Val53) Cd (Ile174) 7Br (K25)

N3 (K25) Cl O (Arg47) Od (Asp175) O (Val116) Benzene(K25) Imidazole (K25) Imidazole (K25) Imidazole (K25) Phenyl group (Phe113)

2.86 3.37 3.60 3.83 3.24 4.17 3.98 3.96 3.33 3.93

2.60 8.89 2.88 2.64 2.69 1.67 2.00 2.05 0.11 2.31

H-bond H-bond CH/O CH/O X-bond CH/p CH/p CH/p CH/p Halogen/p

a D: Distance between atom R of donor (ReH) and atom A of acceptor in the case of hydrogen bonds. In the case of halogen bonds, D correspond to the distance between halogen X and acceptor Y. D is defined for the center of an aromatic ring when the aromatic ring acts as the acceptor. b [kcal/mol]. c Interactions are classified as H-bond (hydrogen bonds), CH/O, CH/p, X-bond (halogen bonds), and halogen/p based on the types of chemical groups of donor and acceptor and D value.

forming area is formed around the phenyl group of Phe113 (Fig. 6b), indicating that the calculated MIF correctly identifies the halogen/p interaction area. As mentioned above, the MIFs correctly identify the nonconventional-interaction-forming areas, as well as conventional hydrogen-bond-forming areas around the target ligand molecule. The results of the analyses for the other six inhibitors are listed in the supporting information (Tables S1eS7). 3.2. Prediction of CK2/inhibitor interaction energies based on MIFs and coarse-grained protein model The binding affinities of candidate compounds are usually estimated by force-field-based or knowledge-based scoring functions in SBDD [2,46]. The force-field-based scoring function is based on the proteineligand interaction energy (Eint) [46]. Applying the MIF to SBDD, thus, requires the evaluation of the proteineligand interaction energy (Eint) based on the MIF. This section shows the calculated Eint for the CK2/inhibitor complexes based on eqs (1) and (2). Additionally, the correlations between the calculated Eint values and experimental log(Ki) values are demonstrated here. We first calculated the values of Eint for the CK2/inhibitor complexes using only the property spheres defined on the amino acid residues, which are denoted Eint(pro) values. This calculation ignores the bridging/interstitial waters and ions. The calculated Eint(pro) values are listed in Table 4. Additionally, the correlation between the calculated Eint(pro) and log(Ki) is shown in Fig. 7a. The Eint(pro) values of K22 and TBS shift to the lower energy region, whereas the others show a good correlation with the experimental log(Ki) (coefficient of determination (r2) ¼ 0.15 (all inhibitors), r2 ¼ 0.70 (excluding Table 4 The calculated CK2/inhibitor interaction energies and the solvation energies of the inhibitors at the uB97XD/6-31G(d,p) level of theory. Entry

Compound

PDB ID

Ki [nM]

Eint(pro)

1 2 3 4 5 6 7

K25 K37 K44 K32 K22 K17 TBS

1ZOE 1ZOG 1ZOH 2OXD 2OXX 2OXY 1J91

40 70 100 150 200 300 400

53.42 45.39 48.23 43.81 58.65 43.33 69.96

a

Eint(all)

b

65.45 59.59 54.82 51.15 80.21 49.48 81.27

DEsolvc 11.17 9.04 11.79 13.01 49.38 10.62 48.51

a Calculated interaction energy without bridging/interstitial waters and ions (BW) ([kcal/mol]). b Calculated interaction energy with BW ([kcal/mol]). c Solvation energy of the ligand.

9

K22 and TBS)). Notably, K22 and TBS are negatively charged, whereas the other compounds are neutral. Therefore, the Eint(pro) value is a good index for the ranking of the affinity for the neutral ligands but is not appropriate for the ligand set containing both the neutral and charged ligands. We considered that the effects of the Coulomb energy of charged ligands may be the major cause of the shift to the lower energy region because the same trend in the interaction energies of charged ligands has been reported in previous studies [91]. Charged ligands interact with proteins more strongly than neutral ligands, but they also form strong interactions with water molecules in aqueous solution. Therefore, the desolvation of charged ligands upon complex formation would lead to a larger enthalpic loss. The ranking of the ligand affinity by Eint(pro) values alone, thus, overestimates of the affinity of the charged ligands because the effects of desolvation are ignored in this calculation. Therefore, the energy values obtained by subtracting the solvation energy of the ligands (DEsolv) from Eint(pro) would be more appropriate for ranking the ligand affinity. The solvation energies of the inhibitors (DEsolv) were calculated by the SMD method at the uB97XD/6-31G(d,p) level of theory (Table 4) [92]. In contrast, we assumed that the effects of desolvation for protein was negligible in our system because all inhibitors target the same protein. The correlation between the calculated Eint(pro) e DEsolv and the experimental log(Ki) is shown in Fig. 7b. As shown, the correlation was improved (r2 ¼ 0.44 (all inhibitors)), indicating that desolvation must be considered. The bridging/interstitial water molecules and ions observed in the crystal structure of the CK2/K25 complex (Fig. 5) are considered to affect the binding affinity of the inhibitors to proteins [93,94]. Therefore, Eint and Eint  DEsolv were next calculated for the coarsegrained models where the property spheres are defined on the bridging/interstitial water molecules and ions that directly interact with the inhibitors, as well as protein, denoted Eint(all) and Eint(all)  DEsolv, respectively. These models regard the bridging/ interstitial waters and ions as the parts of the protein. The bridging/ interstitial waters and ions are listed in Table S8 in detail. Rec(DA), Rec(HD), and Rec(HA) spheres are defined on the bridging/interstitial water molecules and positive and negative ions, respectively. We assumed that the interaction formed between the ligand and the Naþ cation can be explored by the HD probe because Naþ ion acts as Lewis acid as with protons, although cations cannot act as hydrogen bond donors. Therefore, we defined Naþ cations with a Rec(HD) sphere. Even when considering water and ions as part of the protein, equation (2) could be applied to calculate Eint(all) values in the same manner as the bare protein. The correlation between the calculated Eint(all) values and experimental log(Ki) values is shown in Fig. 7c. The correlation trend is similar to that observed in Fig. 7a because the desolvation energy is not corrected (r2 ¼ 0.06 (all inhibitors), r2 ¼ 0.92 (only neutral inhibitors)). The correlation between the calculated Eint(all)  DEsolv and the experimental log(Ki) values is shown in Fig. 7d. Good correlation (r2 ¼ 0.80) was obtained, indicating that the consideration of desolvation and bridging/interstitial water molecules and ions improves the correlation dramatically. To date, there have been several QM based studies that have investigated the correlation between the calculated proteineligand interaction energies and experimental inhibitory activity (log(Ki) and log(IC50)). These previous studies investigated different proteineligand system with 5e60 ligands using QM/molecular mechanics (MM) or fragment-based methods. The r2 values reported in these studies range from 0.58 to 0.94 [38,39,44,91,95e98]. The correlation between the proteineligand interaction energy and log(Ki) in the present study is comparable to that proposed in the previous studies. Finally, we performed additional calculations for Cyclin-dependent kinase 2/inhibitor

10

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

Fig. 7. Correlation between experimentally measured log(Ki) values and the calculated values of Eint and Eint e DEsolv for CK2 inhibitors. The correlation between log(Ki) and Eint(pro) (a) and log(Ki) and Eint(pro) e DEsolv (b) for the CK2/inhibitor complex. Correlation between log(Ki) and Eint(all) (c) and log(Ki) and Eint(all) e DEsolv (d) for CK2/inhibitor/bridginginterstitial-water molecules-and-ions multi-complex. r2 values calculated for only neutral inhibitors are noted in brackets. The charged ligands are underlined.

complexes to confirm the applicability of our method in another ligand-protein system. The results are shown in the Figiure S9. The calculated Eint(all) show good correlation with the experimental DGbind (r2 ¼ 0.80). These results suggest the validity of our method. Therefore, Eint calculated in the manner proposed here would be promising as a descriptor for protein-ligand interaction prediction considering nonconventional interactions. In the future, the combination between this new MIF calculation and QSAR modeling methods is expected to improve ligand affinity predictions not only in SBDD but also in ligand-based drug design (LBDD). Furthermore, the calculated MIFs for the compounds with active conformations can be diverted to the prediction of the proteineligand interactions for arbitrary proteins without additional QM calculations. It is also expected that the construction of the MIF database will make fast QM level estimations of the proteineligand interactions that explicitly describe the possible nonconventional interactions. 4. Conclusions We have developed a method for calculating MIFs using QM calculations with four types of probe molecule. The calculated MIFs for CK2 inhibitors correctly describe the hydrogen bond, halogen bond, CH/p, and CH/O interactions formed in CK2/inhibitor complexes. The method, thus, is promising for the construction of QSAR

and the other modeling method with consideration of nonconventional interactions. Additionally, we have developed a calculation method for proteineligand interaction energies based on the MIFs and coarse-grained protein models. The calculated interaction energies with considering the desolvation of ligands and the bridging/interstitial waters and ions correlated well with experimental log(Ki) values of CK2 inhibitors. These proof-of-concept calculations proposed in the present study suggest that MIF and Eint are promising as descriptors ranking ligand affinities. To confirm the universality of the methodology, the benchmark calculation using various types of protein families is desired. In addition, the combination of this method with chemoinformatics techniques is very promising to progress the accuracy of the prediction of protein-ligand affinity. These would be performed in our future work. In conclusion, the basic methodology to evaluate protein-ligand interactions based on MIF is established in the present study. Our methods for calculating MIF and Eint are promising for the prediction of protein-ligand interactions by considering the nonconventional noncovalent interactions. Funding This work was supported by JSPS KAKENHI Grant Number JP18K14887 (D.H.) and JP19K07029 (H.G.).

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.jmgm.2019.107515. References [1] A.C. Anderson, The process of structure-based drug design, Chem. Biol. 10 (2003) 787e797. [2] C.N. Cavasotto, In Silico Drug Discovery and Design. Theory, Methods, Challenges, and Applications, CRC Press, Boca Raton, 2016. [3] O. Takahashi, Y. Kohno, M. Nishio, Relevance of weak hydrogen bonds in the conformation of organic compounds and bioconjugates: evidence from recent experimental data and high-level ab initio MO calculations, Chem. Rev. 110 (2010) 6049e6076. [4] G. Cavallo, P. Metrangolo, R. Milani, T. Pilati, A. Priimagi, G. Resnati, G. Terraneo, The halogen bond, Chem. Rev. 116 (2016) 2478e2601. [5] H.A. Benesi, J.H. Hildebrand, A spectrophotometric investigation of the interaction of iodine with aromatic hydrocarbons, J. Am. Chem. Soc. 71 (1949) 2703e2707. [6] R.M. Keefer, L.J. Andrews, The interaction of iodine monochloride with benzene and certain of its derivatives, J. Am. Chem. Soc. 72 (1950) 5170e5173. [7] O. Hassel, J. Hvoslef, The structure of bromine 1,4-dioxantane, Acta Chem. Scand. 8 (1954) 873. [8] O. Hassel, K.O. Strømme, Halogen-molecule bridges in solid addition compounds, Nature 182 (1958) 1155e1156. [9] D.J. Sutor, The C-H,,,O hydrogen bond in crystals, Nature 195 (1962) 68e69. [10] D.J. Sutor, Evidence for the existence of C-H,,,O hydrogen bonds in crystals, J. Chem. Soc. 1105 (1963). [11] Y. Kodama, K. Nishihata, M. Nishio, Attractive interaction between aliphatic and aromatic systems, Bull. Chem. Soc. Jpn. 24 (1977) 2105e2108. [12] Y. Kodama, K. Nishihata, S. Zushi, M. Nishio, J. Uzawa, K. Sakamoto, H. Iwamura, The conformation of a diastereoisomeric pair of 2,2-Dimethyl-4phenyl-3-pentanols, Bull. Chem. Soc. Jpn. 52 (1979) 2661e2669. [13] M. Nishio, M. Hirota, CH/p interaction: implications in organic chemistry, Tetrahedron 45 (1989) 7201e7245. [14] M. Nishio, M. Hirota, Y. Umezawa, The CH/p Interaction. Evidence, Nature and Consequences, Wiley-VCH, New York, 1998. [15] A. Mukherjee, S. Tothadi, G.R. Desiraju, Halogen bonds in crystal engineering: like hydrogen bonds yet different, Acc. Chem. Res. 47 (2014) 2514e2524. [16] D. Braga, F. Grepioni, E. Tedesco, X-H,,,p(X ¼ O, N, C) hydrogen bonds in organometallic crystals, Organometallics 17 (1998) 2669e2672. [17] R. Taylor, O. Kennard, Crystallographic evidence for the existence of C-H,,,O, C-H,,,N, and C-H,,,Cl hydrogen bonds, J. Am. Chem. Soc. 104 (1982) 5063e5070. [18] G.R. Desiraju, Distance dependence of C-H,,,O interactions in some chloroalkyl compounds, J. Chem. Soc., Chem. Commun. (1989) 179e180. [19] G.R. Desiraju, Strength and linearity of C-H,,,O bonds in molecular crystals: a database study of some terminal alkynes, J. Chem. Soc., Chem. Commun. (1990) 454e455. [20] T. Steiner, G.R. Desiraju, Distinction between the weak hydrogen bond and the van der Waals interaction, Chem. Commun. (1998) 891e892. [21] Y. Umezawa, S. Tsuboyama, K. Honda, J. Uzawa, M. Nishio, CH/p interaction in the crystal structure of organic compounds. A database study, Bull. Chem. Soc. Jpn. 71 (1998) 1207e1213. [22] F.H. Allen, The Cambridge structural database: a quarter of a million crystal structures and rising, Acta Crystallogr. B 58 (2002) 380. [23] C.R. Groom, F.H. Allen, The Cambridge structural database in retrospect and prospect, Angew. Chem. Int. Ed. 53 (2014) 662e671. [24] H.M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T.N. Bhat, H. Weissig, I.N. Shindyalov, P.E. Bourne, The protein data bank, Nucleic Acids Res. 28 (2000) 235e242. [25] N. Muzet, B. Guillot, C. Jelsch, E. Howard, C. Lecomte, Electrostatic complementarity in an aldose reductase complex from ultra-high-resolution crystallography and first-principles calculations, Proc. Natl. Acad. Sci. U. S. A 100 (2004) 8742e8747. [26] E.I. Howard, R. Sanishvili, R.E. Cachau, A. Mitschler, B. Chevrier, P. Barth, V. Lamour, M. Van Zandt, E. Sibley, C. Bon, D. Moras, T.R. Schneider, A. Joachimiak, A. Podjarny, Ultrahigh resolution drug design I: details of interactions in human aldose reductase-inhibitor complex at 0.66 Å, Proteins Struct. Funct. Genet. 55 (2004) 792e804. [27] P. Auffinger, F.A. Hays, E. Westhof, P.S. Ho, Halogen bonds in biological molecules, Proc. Natl. Acad. Sci. U. S. A 101 (2004) 16789e16794. [28] Z.S. Derewenda, L. Lee, U. Derewenda, The occurrence of C-H,,,O hydrogen bonds in proteins, J. Mol. Biol. 252 (1995) 248e262.

11

[29] M.C. Wahl, M. Sundaralingam, C-H,,,O hydrogen bonding in biology, Trends Biochem. Sci. 22 (1997) 97e102. [30] M.S. Weiss, M. Brandl, J. Sühnel, D. Pal, R. Hilgenfeld, More hydrogen bonds for the (structural) biologist, Trends Biochem. Sci. 26 (2001) 521e523. [31] T. Steiner, G. Koellner, Hydrogen bonds with p-acceptors in proteins: frequencies and role in stabilizing local 3D structures, J. Mol. Biol. 305 (2001) 535e557. [32] Y. Umezawa, M. Nishio, CH/p interactions in the crystal structure of TATA-Bos binding protein/DNA complexes, Bioorg. Med. Chem. 8 (2000) 2643e2650. [33] M. Brandl, M.S. Weiss, A. Jabs, J. Sühnel, R. Hilgenfeld, C-H,,,p-Interactions in proteins, J. Mol. Biol. 307 (2001) 357e377. [34] M. Nishio, Y. Umezawa, J. Fantini, M. Weiss, P. Chakrabarti, CH-p Hydrogen bonds in biological macromolecules, Phys. Chem. Chem. Phys. 16 (2014) 12648e12683. [35] T. Ozawa, K. Okazaki, K. Kitaura, Importance of CH/p hydrogen bonds in recognition of the core motif in proline-recognition domains: an ab initio fragment molecular orbital study, J. Comput. Chem. 32 (2011) 2774e2781. [36] T. Ozawa, K. Okazaki, K. Kitaura, CH/p hydrogen bonds play a role in ligand recognition and equilibrium between active and inactive states of the b2 adrenergic receptor: an ab initio fragment molecular orbital (FMO) study, Bioorg. Med. Chem. 19 (2011) 5231e5237. [37] J. Koseki, H. Gouda, S. Hirono, Molecular orbital study of The formation of intramolecular hydrogen bonding of a ligand molecule in a protein aromatic hydrophobic pocket, Chem. Pharm. Bull. 64 (2016) 1031e1035. [38] C. Watanabe, H. Watanabe, K. Fukuzawa, L.J. Parker, Y. Okiyama, H. Yuki, S. Yokoyama, H. Nokano, S. Tanaka, T. Honma, 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, J. Chem. Inf. Model. 57 (2017) 2996e3010. [39] M. Ozawa, T. Ozawa, K. Ueda, Application of the fragment molecular orbital method analysis to fragment-based drug discovery of BET (bromodomain and extra-terminal proteins) inhibitors, J. Mol. Graph. Model. 74 (2017) 73e82. [40] M. Ozawa, T. Ozawa, M. Nishio, K. Ueda, The role of CH/p interactions in the high affinity binding of streptavidin and biotin, J. Mol. Graph. Model. 75 (2017) 117e124. [41] D. Hayakawa, K. Ueda, C. Yamane, H. Miyamoto, F. Horii, Molecular Dynamics simulation of the dissolution process of a cellulose triacetate-II nano-sized crystal in DMSO, Carbohydr. Res. 346 (2011) 2940e2947. [42] T. Kobayashi, D. Hayakawa, T. Khishigjargal, K. Ueda, Investigation of the structure and interaction of cellulose triacetate I crystal using ab initio calculations, Carbohydr. Res. 388 (2014) 61e66. [43] T. Ishikawa, D. Hayakawa, H. Miyamoto, M. Ozawa, T. Ozawa, K. Ueda, Ab initio studies on the structure of and atomic interactions in cellulose IIII crystals, Carbohydr. Res. 417 (2015) 72e77. [44] Z. Xu, Z. Liu, T. Chen, T.T. Chen, Z. Wang, G. Tian, J. Shi, X. Wang, Y. Lu, X. Yan, G. Wang, H. Jiang, K. Chen, S. Wang, Y. Xu, J. Shen, W. Zhu, Utilization of halogen bond in lead optimization: a case study of rational drug design of potent Phosphodiesterase type 5 (PDE5) inhibitors, J. Med. Chem. 54 (2011) 5607e5611. [45] P.J. Goodford, A computational procedure for determining energetically favorable binding sites on biologically important macromolecules, J. Med. Chem. 28 (1985) 849e857. €ltje, W. Sippl, D. Rognan, G. Folkers, Molecular Modeling. Basic [46] H.-D. Ho Principles and Applications Third, Revised and, Expanded Edition, Wiley-VCH, Weinheim, 2008. [47] J.P. Doucet, A. Panaye, Three Dimensional QSAR. Application in Pharmacology and Toxicology, CRC Press, Boca Raton, 2010. [48] R.D. Cramer, D.E. Patterson, J.D. Bunce, Comparative molecular field analysis (CoMFA). 1. Effect of shape on binding of steroids to carrier proteins, J. Am. Chem. Soc. 110 (1988) 5959e5967. [49] G. Klebe, U. Abraham, T. Mietzner, Molecular similarity indices in a comparative analysis (CoMSIA) of drug molecules to correlate and predict their biological activity, J. Med. Chem. 37 (1994) 4130e4146. [50] N. Yamaotsu, A. Oda, S. Hirono, Determination of ligandebinding sites on proteins using long-range hydrophobic potential, Biol. Pharm. Bull. 31 (2008) 1552e1558. [51] A. Oda, N. Yamaotsu, S. Hirono, Evaluation of the searching abilities of HBOP and HBSITE for binding pocket detection, J. Comput. Chem. 30 (2009) 2728e2737. [52] P. Broto, G. Moreau, C. Vandycke, Molecular structures: perception, autocorrelation descriptor and SAR studies. System of atomic contribution for the calculation of the n-octanol/water partition coefficients, Eur. J. Med. Chem. 19 (1984) 71e78. [53] P. Furet, A. Sele, N.C. Cohen, 3D molecular lipophilicity potential profiles: a new tool in molecular modeling, J. Mol. Graph. 6 (1988) 182e189. re, P. Quarendon, L. Kaetterer, Estimating and representing hy[54] J.-L. Fauche drophobicity potential, J. Mol. Graph. 6 (1988) 203e206. [55] G.E. Kellogg, S.F. Semus, D.J. Abraham, HINT: a new method of empirical hydrophobic field calculation for CoMFA, J. Comput. Aided Mol. Des. 5 (1991) 545e552. [56] P. Gaillard, P.-A. Carrupt, B. Testa, A. Boudon, Molecular lipophilicity potential, A tool in 3D QSAR: method and applications, J. Comput. Aided Mol. Des. 8 (1994) 83e96. [57] B. Testa, P.-A. Carrupt, P. Gaillard, F. Billois, P. Weber, Lipophilicity in molecular modeling, Pharm. Res. 13 (1996) 335e343.

12

D. Hayakawa et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107515

[58] M. Bradley, C.L. Waller, Polarizability fields for use in three-dimensional quantitative structureeactivity relationship (3D-QSAR), J. Chem. Inf. Comput. Sci. 41 (2001) 1301e1307. [59] K.H. Kim, G. Greco, E. Novellino, C. Silipo, A. Vittoria, Use of the hydrogen bond potential function in a comparative molecular field analysis (CoMFA) on a set of benzodiazepines, J. Comput. Aided Mol. Des. 7 (1993) 263e280. [60] S. Sakaki, K. Kato, T. Miyazaki, Y. Musashi, K. Ohkubo, H. Ihara, C. Hirayama, Structures and binding energies of benzeneemethane and benzeneebenzene complexes. An ab initio SCF/MP2 study, J. Chem. Soc., Faraday Trans. 89 (1993) 659e664. [61] T. Clark, M. Hennemann, J.S. Murray, P. Politzer, Halogen bonding: the s-hole, J. Mol. Model. 13 (2007) 291e296.   , M.C. Concha, F.M. Ramos, [62] K.E. Riley, J.S. Murray, J. Fanfrlík, J. Rez a c, R.J. Sola P. Politzer, Halogen bond tunability I: the effects of aromatic fluorine substitution on the strengths of halogen-bonding interactions involving chlorine, bromine, and iodine, J. Mol. Model. 17 (2011) 3309e3318. [63] P. Politzer, J.S. Murray, T. Clark, Halogen bonding and other s-hole interactions: a perspective, Phys. Chem. Chem. Phys. 15 (2013) 11178e11189. [64] R. Wilcken, M.O. Zimmermann, A. Lange, S. Zahn, B. Kirchner, F.M. Boeckler, Addressing methionine in molecular design through directed sulfurehalogen bonds, J. Chem. Theory Comput. 7 (2011) 2307e2315. [65] A. Lange, M.O. Zimmermann, R. Wilcken, S. Zahn, F.M. Boeckler, Targeting histidine side chains in molecular design through nitrogenehalogen bonds, J. Chem. Inf. Model. 53 (2013) 3178e3189. [66] R. Wilcken, M.O. Zimmermann, A. Lange, A.C. Joerger, F.M. Boeckler, Principles and applications of halogen bonding in medicinal chemistry and chemical biology, J. Med. Chem. 56 (2013) 1363e1388. [67] J.-D. Chai, M. Head-Gordon, Systematic optimization of long-range corrected hybrid density functionals, J. Chem. Phys. 128 (2008), 084106. [68] J.-D. Chai, M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections, Phys. Chem. Chem. Phys. 10 (2008) 6615e6620. [69] S. Kozuch, J.M.L. Martin, Halogen bonds: benchmarks and theoretical analysis, J. Chem. Theory Comput. 9 (2013) 1918e1931.  Va zquez-Mayagoitia, B.G. Sumpter, C.D. Sherrill, Density-func[70] L.A. Burns, A. tional approaches to noncovalent interactions: a comparison of dispersion corrections (DFT-D), exchange-hole dipole moment (XDM) theory, and specialized functionals, J. Chem. Phys. 134 (2011), 084107. [71] Gaussian 09, Revision D.01, Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dap€ Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; prich, S.; Daniels, A. D.; Farkas, O.; Fox, D. J. Gaussian, Inc., Wallingford CT, 2009. [72] K. Iwase, S. Hirono, Estimation of active conformations of drugs by a new molecular superposing procedure, J. Comput. Aided Mol. Des. 13 (1999) 499e512. [73] N. Yamaotsu, S. Hirono, In silico fragment-mapping method: a new tool for fragment-based/structure-based drug discovery, J. Comput. Aided Mol. Des. 32 (2018) 1229e1245. [74] S. Kmiecik, D. Gront, M. Kolinski, L. Wieteska, A.E. Dawid, A. Kolinski, Coarsegrained protein models and their applications, Chem. Rev. 116 (2016) 7898e7936. [75] R. Battistutta, M. Mazzorana, S. Sarno, Z. Kazimierczuk, G. Zanotti, L.A. Pinna, Inspecting the structureeactivity relationship of protein kinase CK2 inhibitors derived from tetrabromo-benzimidazole, Chem. Biol. 12 (2005) 1211e1219. [76] R. Battistutta, M. Mazzorana, L. Cendron, A. Bortolato, S. Sarno, Z. Kazimierczuk, G. Zanotti, S. Moro, L.A. Pinna, The ATP-binding site of protein kinase CK2 holds a positive electrostatic area and conserved water molecules, Chembiochem 8 (2007) 1804e1809. [77] R. Battistutta, E.D. Moliner, S. Sarno, G. Zanotti, L.A. Pinna, Structural features

[78] [79]

[80]

[81]

[82]

[83]

[84]

[85]

[86]

[87] [88]

[89]

[90]

[91]

[92]

[93]

[94] [95]

[96]

[97]

[98]

underlying selective inhibition of protein kinase CK2 by ATP site-directed tetrabromo-2-benzotriazole, Protein Sci. 10 (2001) 2200e2206. M. Mazzorana, L.A. Pinna, R. Battistutta, A structural insight into CK2 inhibition, Mol. Cell. Biochem. 316 (2008) 57e62. G.M. Sastry, M. Adzhigirey, T. Day, R. Annabhimoju, W. Sherman, Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments, J. Comput. Aided Mol. Des. 27 (2013) 221e234. M.P. Jacobson, D.L. Pincus, C.S. Rapp, T.J.F. Day, B. Honig, D.E. Shaw, R.A. Friesner, A hierarchical approach to all-atom protein loop prediction, Protein Struct. Funct. Bioinform. 55 (2004) 351e367. W.L. Jorgensen, J. Tirado-Rives, The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin, J. Am. Chem. Soc. 110 (1988) 1657e1666. W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids, J. Am. Chem. Soc. 118 (1996) 11225e11236. D. Shivakumar, J. Williams, Y. Wu, W. Damm, J. Shelley, W. Sherman, Prediction of absolute solvation free energies using molecular Dynamics free energy perturbation and the OPLS force field, J. Chem. Theory Comput. 6 (2010) 1509e1519. E. Harder, W. Damm, J. Maple, C. Wu, M. Reboul, J.Y. Xiang, L. Wang, D. Lupyan, M.K. Dahlgren, J.L. Knight, J.W. Kaus, D.S. Cerutti, G. Krilov, W.L. Jorgensen, R. Abel, R.A. Friesner, OPLS3: a force field providing broad coverage of drug-like small molecules and proteins, J. Chem. Theory Comput. 12 (2016) 281e296. J.R. Greenwood, D. Calkins, A.P. Sullivan, J.C. Shelley, Towards the comprehensive, rapid, and accurate prediction of the favorable tautomeric states of drug-like molecules in aqueous solution, J. Comput. Aided Mol. Des. 24 (2010) 591e604. J.C. Shelley, A. Cholleti, L. Frye, J.R. Greenwood, M.R. Timlin, M. Uchimaya, Epik: a software program for pKa prediction and protonation state generation for drug-like molecules, J. Comput.-Aided Mol. Design 21 (2007) 681e691. W. Humphrey, A. Dalke, K. Schulten, Vmd d visual molecular Dynamics, J. Mol. Graph. 14 (1996) 33e38. K.E. Riley, P. Hobza, The relative roles of electrostatics and dispersion in the stabilization of halogen bonds, Phys. Chem. Chem. Phys. 15 (2013) 17742e17751. M.O. Sinnokrot, E.F. Valeev, C.D. Sherrill, Estimates of the ab initio limit for pp interactions: the benzene dimer, J. Am. Chem. Soc. 124 (2002) 10887e10893. K.E. Riley, C.L. Ford Jr., K. Demouchet, Comparison of hydrogen bonds, halogen bonds, CH,,,p interactions, and C-X,,,p interactions using high-level ab initio methods, Chem. Phys. Lett. 621 (2015) 165e170. K. Maruyama, Y. Sheng, H. Watanabe, K. Fukuzawa, S. Tanaka, Application of singular value decomposition to the inter-fragment interaction energy analysis for ligand screening, Comput. Theor. Chem. 1132 (2018) 23e34. A.V. Marenich, C.J. Cramer, D.G. Truhlar, Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions, J. Phys. Chem. B 113 (2009) 6378e6396. J. Liu, X. He, J.Z.H. Zhang, Improving the scoring of proteineligand binding affinity by including the effects of structural water and electronic polarization, J. Chem. Inf. Model. 53 (2013) 1306e1314. G. Klebe, Applying thermodynamic profiling in lead finding and optimization, Nat. Rev. Drug Discov. 14 (2015) 95e110. A. Heifetz, G. Trani, M. Aldeghi, C.H. MacKinnon, P.A. McEwan, F.A. Brookfield, E.I. Chudyk, M. Bodkin, Z. Pei, J.D. Burch, D.F. Ortwine, Fragment molecular orbital method applied to lead optimization of novel interleukin-2 inducible T-cell kinase (ITK) inhibitors, J. Med. Chem. 59 (2016) 4352e4363. M.P. Mazanetz, O. Ichihara, R.J. Law, M. Whittaker, Prediction of cyclindependent kinase 2 inhibitor potency using the fragment molecular orbital method, J. Cheminf. 3 (2011) 2. R. Takeda, I. Kobayashi, R. Suzuki, K. Kawai, A. Kittaka, M. Takimoto-Kamimura, N. Kurita, Proposal of potent inhibitors for vitamin-D receptor based on ab initio fragment molecular orbital calculations, J. Mol. Graph. Model. 80 (2018) 320e326. B. Thapa, D. Beckett, J. Erickson, K. Raghavachari, Theoretical study of proteineligand interactions using the molecules-in-molecules fragmentationbased method, J. Chem. Theory Comput. 14 (2018) 5143e5155.