Theoretical research in structure characteristics of different inhibitors and differences of binding modes with CBP bromodomain

Bioorganic & Medicinal Chemistry 26 (2018) 712–720

Contents lists available at ScienceDirect

Bioorganic & Medicinal Chemistry journal homepage: www.elsevier.com/locate/bmc

Theoretical research in structure characteristics of different inhibitors and differences of binding modes with CBP bromodomain Xue-Song Wang b, Qing-Chuan Zheng a,b,⇑ a

Key Laboratory for Molecular Enzymology and Engineering of the Ministry of Education, Jilin University, Changchun 130012, People’s Republic of China Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, International Joint Research Laboratory of Nano-Micro Architecture Chemistry, Jilin University, Changchun 130023, People’s Republic of China b

a r t i c l e

i n f o

Article history: Received 25 September 2017 Revised 18 December 2017 Accepted 24 December 2017 Available online 26 December 2017 Keywords: Bromodomain Inhibitor MD simulations LUDI Rational design

a b s t r a c t The CBP (CREB (cAMP responsive element binding protein) binding protein) bromodomain (BRD) could recognize and bind with acetyl K382 of human tumor suppressor protein p53 which the mutation of encoding gene might cause human cancers. CBP-BRD serves as a promising drug target for several disease pathways and a series of effective drug have been discovered. In this study, molecular dynamics (MD) simulations and molecular mechanics generalized born surface area (MM-GB/SA) approaches were performed to investigate the different binding modes between five inhibitors with CBP-BRD. Based on the energy and conformation analyses, a potent core fragment is chosen to act as the starting point for new inhibitor design by means of LUDI and rational drug design approaches. Then, T.E.S.T and molinspirition were applied to evaluate oral bioavailability and drug promiscuity of the new molecules. These results shed light on the idea for further inhibitor design. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The CBP (CREB (cAMP responsive element binding protein) binding protein) acts as lysine acetyl transferase which plays a key role in gene posttranslational modification in human cells.1 CBP possesses a single bromodomain (BRD) that could recognize and bind acetyl lysine of histones and non-histone proteins.2 The lysine 382 (K382) of human tumor suppressor protein p53 acting as gene-specific transcription factors could be recruited by CBPBRD.3 The mutations of p53 encoding gene are frequent in most human cancers such as leukemias, lymphomas, ischemic brain injury, and auto-immune diseases, etc., causing the binding proAbbreviations: CBP, CREB (cAMP responsive element binding protein) binding protein; BRD, bromodomain; CPI-098, (4R)-4-methyl-1,3,4,5-tetrahydro-2H-1,5benzodiazepin-2-one; CPI-703, (4R)-6-(1-tert-butyl-1H-pyrazol-4-yl)-4-methyl1,3,4,5-tetrahydro-2H-1,5-benzodiazepin-2-one; CPI-637, (4R)-4-methyl-6-[1methyl-3-(1-methyl-1H-pyrazol-4-yl)-1H-indazol-5-yl]-1,3,4,5-tetrahydro-2H-1, 5-benzodiazepin-2-one; GNE-272, 1-[3-[[2-fluoranyl-4-(1-methylpyrazol-4-yl) phenyl]amino]-1-[(3
{S})-oxolan-3-yl]-6,7-dihydro-4
{H}-pyrazolo[4,3-c]pyridin-5-yl]ethanone; SGC-CBP30, , 2-[2-(3-chloro-4-methoxyphenyl)ethyl]-5(3,5-dimethyl-1,2-oxazol-4-yl)-1-[(2S)-2-(morpholin-4-yl)propyl]-1H-benzimidazole; MD, molecular dynamics; MM-GB/SA, molecular mechanics generalized born surface area; LD50, lethal dose 50%. ⇑ Corresponding author at: Institute of Theoretical Chemistry, Jilin University, Changchun 130023, People’s Republic of China. E-mail address: [email protected] (Q.-C. Zheng). https://doi.org/10.1016/j.bmc.2017.12.040 0968-0896/Ó 2017 Elsevier Ltd. All rights reserved.

gress between CBP-BRD and K382 over-expressed.4 The BRDs have been regarded as a promising drug target.5 BRD is made up of
110 residues.6 There are 61 known BRDs sharing the conserved fold, which contains a left-handed bundle of four a-helices (aZ, aA, aB, aC), in spite of the large sequence variations.7–9 The ZA-loop and BC-loop form a largely hydrophobic acetyl-lysine binding site.5,10 In CBP-BRD, the NH2 of N1168 (at the beginning of the BC-loop) binds the N-acetyl group of acetyl-lysine via hydrogen bond, and the phenolic hydroxyl group of Y1125 (in the ZA-loop) forms a water-mediated hydrogen bond.11 Aiming to modulate the connection and function of the acetyl lysine with BRDs, a large number of inhibitors have been developed.5,11 A class of BRD inhibitors directly engage the binding site in the way that competitively inhibit the binding of acetyl-lysine to BRDs.5 In this study, we selected five small molecular inhibitors (CPI-098,12 CPI-703,12 CPI-637,13 GNE-27214 and SGC-CBP3015) that every of them has its own characteristics. The CPI-098 is the template, and CPI-703 gets a small substituent group, while the CPI-637 receives a long substituent. As their study identified, the IC50 values are 14 mm, 0.47 mm and 0.03 mm in turn, decreasing order of magnitude. GNE-272 comes from another paper of the same research group, based on different synthetic method. SGCCBP30 has a better treatment effect as marketed drug. In order to clarify the reason that different structures of inhibitors result prominent inhibitory effect and the detail of binding modes,

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

molecular dynamics (MD) simulations and molecular mechanics generalized born surface area (MM-GB/SA) approaches are employed to investigate the detail interactions between CBP-BRD and five ligands. Based on these discussions, we look forward to get some new clues to design new inhibitors.

713

for constraining bonds involving hydrogen atoms, and sampling every 2 fs. The cluster analyses of five complexes were performed using the average linkage as the clustering algorithm.21 All the data given in the tables or figures were obtained from the last 15 ns (50 ns timescale) or last 5 ns (10 ns timescale) MD simulations and processed by PyMOL,22 Chimera23 and VMD24 softwares.

2. Method 2.3. MM-GB/SA calculations 2.1. Molecular system preparation In this study, the initial five crystal structures were taken from the Protein Data Bank (PDB ID: 4YK0,12 5DBM,12 5I8G,13 5KTX14 and 4NR715). Five inhibitors structures are shown in Fig. 1. Residues number of CBP-BRD starts from 1082 to 1197. As mentioned in Refs. 12–15, waters locating between inhibitors and Try1125 could form water-mediated hydrogen bonds. Thus several crystal water molecules deeply inside the binding pocket were remained and others were deleted using Discovery Studio 2.516 software. 2.2. Molecular dynamics simulations The MD simulations were carried out by AMBER 1117 software package with the ff99SB18 force field. The GAFF19 force field was employed for ligands. Sodium ions were added to keep the system neutral by tleap module in AMBER 11. Each complex was solvated with TIP3P20 waters in a truncted octahedral box added up to a distance of 10 Å from the solute. Both the protein and ligands were fixed with a 500 kcal mol1 Å2 positional restraint, and minimized the energy of the water and the ions for 3000 steps of steepest descent (SD) method and 5000 steps of conjugate gradient (CG) algorithms. Then the minimization was repeated for 3000 steps SD and 5000 steps CG without restraint. After that, the temperature was increased gently to 300 K with all restraints by a 100 kcal mol1 Å2 weights and then equilibrated for 200 ps. Finally, 50 ns MD simulations were performed for five complex systems under NPT ensemble conditions; on the other hand, the complexes with new design ligands performed 10 ns MD simulations. ParticleMesh Ewald technique was used to calculate the long range electrostatic interactions with a non-bonded cutoff of 10.0 Å to limit the direct space sum. During dynamics SHAKE algorithm was used

Binding free energy of each complex was calculated by molecular mechanics generalized born surface area (MM-GB/SA)25–28 approaches in AMBER 11. We extracted the last 7500 snapshots from each simulation trajectory to calculate energy. The energy calculation was based on the following equation:

DGbind ¼ Gcomplex  ðGprotein þ Gligand Þ

ð1Þ

Here, Gcomplex, Gprotein, and Gligand are the free energies of the complex, the protein, and ligands, respectively. The free energy (Gx = complex,protein,ligand) of each species can be estimated by using MMGB/SA methods:

G ¼ EMM þ Gsol  TS

ð2Þ

EMM ¼ Eele þ Evdw þ Eint

ð3Þ

Gsol ¼ GGB þ GSA

ð4Þ

Here, EMM is the gas phase molecular mechanical energy, Gsol is the solvation free energy, and Eele, Evdw, and Eint are the electrostatic energy, the van der Waals interaction energy, and the internal energy, respectively. The Gsol can be separated into an electrostatic solvation energy (GGB) and nonelectrostatic solvation energy (GSA). GGB can be calculated with Generalized Born method by AMBER 11. The GSA contribution is considered to be proportional to the molecular solvent accessible surface area buried on binding. Normalmode analysis was performed to estimate the change in conformational entropy upon ligands binding (TS) for all atoms by using the nmode module of AMBER 11. We took five hundred snapshots from the last 15 ns trajectory at an interval of 5 snapshots to calculate the entropic contribution.

Fig. 1. Structures of five CBP-BRD inhibitors used in this study.

714

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

For obtaining the detailed view of protein and ligands interaction, MM-GB/SA method was employed to calculate the binding free energy of each residue. Each residue can be separated into two parts in the calculation, backbone and side chain, both the backbone and side chain use the same calculation method. 2.4. Structure-based design LUDI29–31 is a fragment-based de novo drug design algorithm. Based on our computational results, the core fragment of CPI-637 was chosen as the starting point for drug design. LUDI was performed to add new fragment onto template and run on Discovery Studio 2.5. In two searching attempts, Gln1113 and Phe1177 were defined as the center of search respectively. The Linkage parameter was set to use all hydrogen atoms. The remaining parameters were set at default settings. Gaussian0932 software was used for molecule geometry optimization by density functional theory33,34 B3LYP35/6-31G.36 2.5. Toxicity test Oral rat LD50 of all ligands were measured by performing the consensus method of T.E.S.T.37 CLogP of all ligands were calculated by carrying Molinspirition38,39 mapping compounds structures on line.

Fig. 3. RMSF value (Å) at each residue position for five complexes.

3. Results and discussion

plex show higher flexibility as might be expected. These residues in 4YK0 complex fluctuate especially higher than in other complexes. CPI-098 is a short-chain inhibitor and could not well stabilize CBPBRD, causing the increase in the structural motion of partial residues.

3.1. Structural stability during MD simulations

3.2. Energetic analysis of complex

The detailed information of structural stability was investigated by calculating average root-mean-square deviation (RMSD) of backbone atoms (Fig. 2). RMSD values of four complexes (4YK0, 5DBM, 5I8G and 5KTX) oscillate ranging from 0.5 Å to 2.25 Å during the first 35 ns, and finally stabilize around 1.5 Å during the last 15 ns. For 4NR7 complex, RMSD values are apparently lower than others during the whole simulations, and finally get stable around 1.0 Å during the last 15 ns. This observation demonstrates that 4NR7 complex displays the best structural stability. For all systems, the last 15 ns trajectories are selected for further analyses. To analyze the variation of flexibility in five complexes, average root-mean-square fluctuations (RMSF) of backbone atoms were calculated (Fig. 3). A similar pattern of flexibility is seen in different complexes. Residues (1112-1123) located on ZA-loop of 4YK0 com-

3.2.1. Binding free energy analysis The binding free energies (DGbind) for five complexes were estimated from snapshots within the last 15 ns trajectories by MMGB/SA approach. The predicted DGbind for five complexes are summarized in Table 1 that are 8.05 kcal mol1, 13.13 kcal mol1, 17.76 kcal mol1, 14.37 kcal mol1 and 19.35 kcal mol1 respectively, in agreement with the variation tendency of experimental results.12–15 In five complexes, DGnonp provides the major favorable contribution to the binding free energy. DGnonp is the sum of DEvdw and DGSA, and the most of its contribution is from DEvdw. The values of DEvdw term in five complexes are decreasing in turn, with the increase of chain length and number of aromatic rings for five inhibitors. While DGpol which is the sum of DEele and DGSA makes the unfavorable contribution to binding free energy, owing to that strongly unfavorable DGSA counteracts all of DEele. The contribution of entropic changes (TS) impair the binding between inhibitors and CBP-BRD. In general, favorable van der Waals and nonpolar energy are the main driving force for the binding between ligands and CBP-BRD.

Fig. 2. RMSD value (Å) for five complexes during 50 ns simulations.

3.2.2. Per residue energy decomposition To further explore the key residues making major contributions to binding, free energy decomposition on per residue was calculated. The values of total interaction energy have been decomposed on a per residue basis into contributions from van der Waals energy, columbic interaction energy, polar solvation free energy, and nonpolar solvation free energy. The residues with free energy decomposition <0.5 kcal mol1 are listed in Table 2 and displayed in Fig. 4. The details of energy decomposition of key residues are listed in Tables S1–S5. Comparing these data, ten amino acid residues (Pro1110, Phe1111, Val1115, Leu1120, Ile1122, Tyr1125, Ala1164, Tyr1167, Asn1168 and Val1174) are observed to play favorable roles in five complexes. Asn1168 and Pro1110 have significant columbic inter-

715

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720 Table 1 Binding free energies between CBP-BRD and ligands in five complexes (kcal mol1).

DEele DEvdw DGGB DGSA DGpola DGnonpb DGMM-GB/SAc –TDS DGbindd a b c d

4YK0

5DBM

5I8G

5KTX

4NR7

17.97 ± 2.65 21.92 ± 2.05 20.04 ± 1.67 3.40 ± 0.09 2.07 25.32 –23.26 ± 2.32 15.20 ± 2.22 8.05 ± 3.21

13.88 ± 3.09 31.79 ± 2.55 18.89 ± 2.02 4.80 ± 0.22 5.01 36.59 31.58 ± 2.88 18.45 ± 1.50 13.13 ± 3.25

14.04 ± 3.43 39.98 ± 2.69 22.14 ± 2.77 5.69 ± 0.16 8.10 45.67 37.56 ± 2.69 19.80 ± 2.19 17.76 ± 3.47

14.57 ± 3.87 39.13 ± 2.58 23.85 ± 3.29 4.40 ± 0.23 9.28 43.53 34.25 ± 2.64 19.88 ± 2.28 14.37 ± 3.49

15.39 ± 3.29 45.09 ± 2.79 24.42 ± 3.12 5.33 ± 0.27 9.03 50.42 41.39 ± 2.94 22.05 ± 3.64 19.35 ± 4.68

Gpol = Eele + GGB. Gnonp = Evdw + GSA. DGMM-GB/SA = Eele + GGB + Evdw + GSA. DGbind = DGMM-GB/SA  TS.

Table 2 Free energy decomposition of several key residues in five complexes (kcal mol1).

Leu 1109 Pro 1110 Phe 1111 Val 1115 Leu 1120 Ile 1122 Tyr 1125 Ala 1164 Tyr 1167 Asn 1168 Arg 1173 Val 1174 Phe 1177

4YK0

5DBM

5I8G

5KTX

4NR7

0.21 1.54 1.50 0.55 1.09 1.20 0.59 0.86 0.84 2.47 0.05 1.11 0.05

0.81 2.09 1.12 0.91 2.45 1.05 0.80 0.82 0.76 3.00 0.14 1.77 0.05

2.69 2.93 1.10 0.96 1.86 1.04 0.73 0.79 0.61 2.84 0.09 2.46 0.98

2.42 2.05 1.11 1.14 2.50 1.22 0.69 0.91 0.77 1.90 0.08 2.43 0.22

2.68 2.08 1.17 1.59 2.39 0.47 0.64 0.90 0.60 1.04 1.06 2.27 0.73

Especially in 4NR7 complex, the free energy contribution of Arg1173 (1.06 kcal mol1) is larger than that of other complexes, which indicates a special contribution of Arg1173 to the binding stability between SGC-CBP30 and CBP-BRD.

3.3. Intermolecular interaction analysis To further investigate the major motivations for five inhibitors binding, the binding site is divided into three domains (Fig. 5), based on the differences of characteristic for side chains of differ-

Fig. 4. Free energy decomposition analyses for five complexes. Residues contributing significantly are highlighted.

action energy contribution. It could be hypothesized that these two residues may facilitate the binding by forming electrostatic interaction. Most of these key residues (Pro1110, Phe1111, Val1115, Leu1120, Ile1122, Tyr1125, Ala1164, Tyr1167 and Val1174) have large van der Waals energies and are involved in the main binding attractions by forming hydrophobic interactions. The free energy decompositions of Leu1109 and Phe1177 in 5I8G, 5KTX and 4NR7 complexes are remarkable lower than in other complexes with short-chain inhibitors. The high van der Waals contributions demonstrate that these two residues may help stabilizing the tail of long-chain inhibitors through the hydrophobic interactions.

Fig. 5. Hydrophobic surface of CBP-BRD. Orange surface represents hydrophobic and blue represents hydrophilic. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

716

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

ent ligands and contributions of hydrogen bonds and hydrophobic interactions from key residues. 3.3.1. Binding mode for domain1 As seen from Fig. 6, these key residues (Phe1111, Val1115, Ile1122, Tyr1125, Ala1164, Tyr1167, Asn1168 and Val1174) located in domain1 play a role of anchoring and stabling the benzodiazepinone template (in CPI-098, CPI-703 and CPI-637), dimethylisoxazole-benzimidazole (in SGC-CBP30), even pyrazolopiperidine (in GNE-272) deep into the pocket.

Hereinto, acting as a highly conserved residue, Asn1168 forms key hydrogen bond interactions with the lactam carbonyl and lactam NH (in CPI-098, CPI-703 and CPI-637), oxazole N and oxazole O (in SGC-CBP30), and amide carbonyl (in GNE-272) (Fig. 7). As shown in Fig. 8, these groups forming hydrogen bonds on five inhibitors display negative electrostatic surface potential. The detail information of hydrogen bonds is listed in Table S6. CPI-098, CPI703, and CPI-637 form four hydrogen bonds with Asn1168 respectively. Though the hydrogen bonds form with the same group, the occupancies for CPI-703 and CPI-637 are obviously larger than CPI-

Fig. 6. Different binding modes of five complexes. (A) CPI-098 binds with residues in domain1 of 4YK0 complex. (B) CPI-098 binds with residues in domain2 of 4YK0 complex. (C) CPI-703 binds with residues in domain1 of 5DBM complex. (D) CPI-703 binds with residues in domain2 of 5DBM complex. (E) CPI-637 binds with residues in domain1 of 5I8G complex. (F) CPI-637 binds with residues in domain2 and domain3 of 5I8G complex. (G) GNE-272 binds with residues in domain1 of 5KTX complex. (H) GNE-272 binds with residues in domain2 and domain3 of 5KTX complex. (I) SGC-CBP30 binds with residues in domain1 of 4NR7 complex. (J) SGC-CBP30 binds with residues in domain2 and domain3 of 4NR7 complex. Red dashed lines represent hydrogen bonds. Black dashed lines represent cation-p interaction.

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

717

Fig. 7. Hydrogen bonds occupancies for five complexes along the last 15 ns simulations.

098. It might mean long-chain structure is essential for inhibitor stability of CBP-BRD. One stable hydrogen bond is monitored between GNE-272 and Asn1168, and two stable hydrogen bonds are monitored between SGC-CBP30 and Asn1168, which occupancies are all greater than 95%. Besides, CPI-098, CPI-703, CPI-637 and GNE-272 also form water-mediated hydrogen bonds with Tyr1125 (Fig. 6A, C, E and G). As shown in Fig. 5, on the both sides of domain1, the dark orange expresses the strong hydrophobic ability of those residues. Hydrophobic residues Val1115, Ile1122, Tyr1125 and Tyr1167 are on one side, and Phe1111, Ala1164, and Val1174 are on the other side. The hydrophobic section forming by these residues also makes large contribution to stabilizing ligands. 3.3.2. Binding mode for domain2 Residues Pro1110, Gln1113 and Leu1120 of domain2 play important roles in stabilizing inhibitors. Pro1110 forms stable hydrogen bond interaction with CPI-098, CPI-703, CPI-637 and GNE-272. The hydrophobic contributions of Pro1110 and Leu1120 help to bind with tert-butyl-pyrazole of CPI-703, indazole of CPI-637, 2-fluorophenyl pyrazole of GNE-272, and benzimidazole with morpholino-ethylene moiety of SGC-CBP30, which could be clarified from van der Waals item of free energy decomposition (Tables S2–S5). Due to the short-chain structure, CPI-098 could not get into the deep position in domain2 and stabilize residues by providing hydrophobic contribution. As shown in Fig. 9C, it could be observed that Leu1120 in 4YK0 complex has obvious movement while it is tightened in other complexes. The columbic interaction energy of Gln1113 in 5KTX complex is 0.55 kcal mol1 (Table S4), and the acylamino group of

Gln1113 could form hydrogen bond with negative charged fluorine of GNE-272 (Fig. 6H), which does not exist in other complexes. As shown in Fig. 5, Gln1113 displays strong hydrophilic, which will be resistive for binding with hydrophobic group of inhibitors such as tert-butyl in CPI-703, methyl in CPI-637 and morpholino-ethylene moiety in SGC-CBP30. These results recommend us that when improving the binding affinity of ligands, it could be useful to add charged substitutes to form more stable interaction with Gln1113. 3.3.3. Binding mode for domain3 Due to the short-chain structures of CPI-098 and CPI-703, these two inhibitors could not bind with the residues in domain3. Thus, we will focus on the discussion about CPI-637, GNE-272 and SGCCBP30. Pro1106, Leu1109 and Phe1177 could form a hydrophobic section to stabilize the end hydrophobic groups such as methylpyrazole of CPI-637 and GNE-272 and chloro-aryl ring of SGCCBP30. The chlorine atom on the aryl ring of SGE-CBP30 inserts into the hydrophobic section between Phe1177 and Val1174, while methoxy group inserts into the hydrophobic section between Leu1109 and Pro1106, helping to stabilize aryl ring in position. These results are also observed in Hay’s study.14 It suggests that hydrophobic and bulky substituent may form more stable interactions with residues in domain3. In one terminal of BC-loop, Arg1173 displays strong positive charged (Fig. 8) and plays different functions in three complexes containing long-chain ligands. In 5I8G complex, Arg1173 swings around the exit of binding site, and can prevent ligand escaping (Fig. 6F). For 5KTX complex, Crawford et al. discovered that Arg1173 might form hydrogen bond with tetrahydrofuran of

718

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

Fig. 8. Electrostatic surface potential of (A) CPI-098, (B) CPI-703, (C) CPI-637, (D) GNE-272, (E) SGC-CBP30, and (F) CBP-BRD. The darker blue regions represent the higher positive electrostatic potential. The darker red regions represent the higher negative electrostatic potential.

Fig. 9. (A) The superimposed representative structures of all complexes. (B) Superimposed structures turn down 90° along the center axis. (C) Enlarge the variation region for clarity.

GNE-272.15 Considering from the distance between Arg1173 and tetrahydrofuran in crystal structure, this hydrogen bond has the possibility of existence. But regretfully in the simulations, electrostatic contribution of Arg1173 is negative and the occupancy of

this hydrogen bond is minimal during the whole simulations time (Fig. 6H). In 4NR7 complex, Arg1173 forms cation-p interaction with aryl ring (Fig. 6J). This is also discovered by Rooney et al40 through experiment.

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

These observations could provide a good strategy for inhibitor design. Forming favorable interaction with Arg1173 or enhancing the hydrophobic interaction within Pro1106, Leu1109 and Phe1177 will be helpful for inhibitor binding with CBP-BRD. 3.4. The structure-based inhibitor design Based on the results mentioned above, though SGC-CBP30 displays the best binding affinity than other inhibitors, the major contributions for binding are from morpholine and chloro-aryl substituents. Benzodiazepinone with methylindazole, the core fragment of CPI-637, could form more stable hydrogen bonds with key residues of domain1 and Pro1109 of domain2 than core fragments of other four inhibitors. Considering of these, the core fragment of CPI-637 was chosen to act as the starting structure for new inhibitor design. 3.4.1. De novo growth of core fragment LUDI was performed to add new fragments onto selected starting structure to increase binding affinity between ligand and CBP-

Table 3 Binding free energies between CBP-BRD and three potential molecules (kcal mol1).

DEele DEvdw DGGB DGSA DGpola DGnonpb DGMM-GB/SAc –TDS DGbindd a b c d

ludi-2

design-5

design-8

17.98 ± 2.54 39.08 ± 2.63 24.42 ± 1.83 4.29 ± 0.24 6.44 43.37 36.94 ± 2.72 19.84 ± 2.70 17.10 ± 3.83

19.50 ± 3.29 42.01 ± 2.62 28.20 ± 2.59 5.99 ± 0.17 8.70 48.00 39.30 ± 2.89 20.53 ± 2.35 18.77 ± 3.72

21.45 ± 2.81 40.17 ± 2.82 28.53 ± 2.17 5.97 ± 0.18 7.08 46.14 39.05 ± 2.95 20.81 ± 1.33 18.24 ± 3.24

Gpol = Eele + GGB. Gnonp = Evdw + GSA. DGMM-GB/SA = Eele + GGB + Evdw + GSA. DGbind = DGMM-GB/SA  TS.

Fig. 10. Hydrogen bonds interactions between design-5 with CBP-BRD. Red dashed lines represent hydrogen bonds.

719

BRD. First, we set the de novo link site containing key residues of domain2. Due to the space constraint, the first attempt does not give us a reasonable result. For the second attempt, the de novo link site is set to contain key residues of domain3. According to the LUDI-guided identification of interaction fragments, five molecules were obtained (Fig. S1). Then, 10 ns MD simulations were performed and the RMSD values (Fig. S2) were calculated to determine the structural stability of five complexes which were obtained by docking molecules (ludi-1, ludi-2, ludi-3, ludi-4 and ludi-5) into CBP-BRD respectively. As shown in Fig. S2, RMSD values of five molecules do not display obvious fluctuation and finally stabilize around 1.5 Å, which means that these five complexes are stable. To further distinguish the binding affinity of above five molecules, the last 5 ns trajectories were selected to calculate the binding free energy by MM-GB/ SA approach. DGbind value of ludi-2 (17.1 kcal mol1) is lower than that of other molecules and basically equals to the data of CPI-637. DGbind energies of other four molecules are much higher than that of CPI-637 (Table S7). These data demonstrate that the binding affinity for ludi-2 is the best in all design molecules. New fragments added onto the starting fragment could not form more favorable interaction in accordance with our design. 3.4.2. Rational design of inhibitors To obtain molecules with higher binding affinity, rational design of inhibitors was performed. Following are four main ideas to form more stable interaction with Gly1113 of domain2 and key residues of domain3. (1) Add negatively charged group on indazole to form more favorable interaction with Arg1173. (2) Change substitute of indazole to interact with Gln1113. (3) Add negatively charged group on the benzene ring of benzodiazepinone to stabilize Arg1173. (4) Replace methylpyrazole by hydrophobic and bulky groups to form more stable interaction with residues in domain3. Thus, eight rational ligands are obtained (Fig. S1) and respectively docked into CBP-BRD. Then, 10 ns MD simulations were performed and the RMSD values were calculated (Fig. S2). The DGbind was calculated from last 5 ns trajectories by MMGB/SA approach. Two better designs are selected according to the binding free energies compared with the value of CPI-637 (Table S8). Hereinto, DGbind of design-5 and design-8 are calculated to be 18.77 kcal mol1 and 18.24 kcal mol1 respectively. The van der Waals energy and electrostatic energy contributions of these two inhibitors are both increased (Table 3). As shown in Fig. 10, design-5 could form one more stable hydrogen bond with Leu1120 which plays an important role in ZA-loop. The occupancy of this hydrogen bond reaches 74.17% (Table S9). It could be considered that these two designs possess potential value for further discussion. 3.4.3. Toxicity test To evaluate the drug safety of our design results, two following method were performed to test toxicity of ludi-2, design-5 and design-8 (Fig. 11). (1) Oral rat LD50 of ligands is calculated by performing T.E.S.T.37 There is no limit to the value of oral rat LD50, so

Fig. 11. Molecular structures of three new potential ligands.

720

X.-S. Wang, Q.-C. Zheng / Bioorganic & Medicinal Chemistry 26 (2018) 712–720

Table 4 The results of MM-GB/SA and toxicity test for five inhibitors and three potent design molecules. Inhibitors

DG

CPI-098 CPI-703 CPI-637 GNE-272 SGC-CBP30 ludi-2 design-5 design-8

8.05 13.13 17.76 14.37 19.35 17.10 18.77 18.24

bind

(kcal mol1)

LD50 (mol kg1)

cLogP

1.90 2.45 2.46 2.34 3.10 2.76 N/A 2.43

3.28 4.59 4.81 2.24 5.08 6.89 3.89 6.96

the LD50 value of CPI-637 is used as the standard for comparing (Table 4). The LD50 value of design molecule will be acceptable if it is higher than the value of CPI-637. (2) CLogP (octanol/water partition coefficient) is calculated by Molinspiration ToxPredic38,39 which represents the tendency of binding between octanol molecular with multiple targets. This data is generally within 82 to 5, and the higher data means higher toxicity. All results are shown in Tables S7 and S8. The oral rat LD50 of ludi-2 and design-8 molecules are similar to CPI-637. While the data for design-5 is undefined. This indicates that ludi-2 and design-8 gain better oral bioavailability in rat than design-5. It could be concluded that substituents in R2 position will decrease the oral bioavailability. However, only cLogP of design-5 satisfies the limited range. The cLog P values of ludi-2 and design-8 are out of the maximum range. It means that ludi-2 and design-8 could bind with multiple targets and display higher promiscuity than design-5. These computational results may provide a potential and convenient way for further drug design. 4. Conclusion Molecular dynamics (MD) simulations and molecular mechanics generalized born surface area (MM-GB/SA) approaches were conducted to investigate the binding modes between five inhibitors (CPI-098, CPI-703, CPI-637, GNE-272 and SGC-CBP30) and CBP (CREB (cAMP responsive element binding protein) binding protein) bromodomain (BRD). To analyze the details of different binding modes, the binding site is separated into three parts (domain1, domain2 and domain3) according to the configurations of inhibitors and contributions of key residues. The computational results show that long-chain inhibitors (CPI-637, GNE-272 and SGC-CBP30) have higher binding affinity toward CBP-BRD than short-chain inhibitors (CPI-098 and CPI-703) does. Some key residues (Leu1109, Pro1110, Phe1111, Val1115, Leu1120, Tyr1125, Ala1164, Tyr1167, Val1174 and Phe1177) interact with inhibitors by providing hydrophobic contributions and several key residues (Pro1110, Gln1113, Asn1168 and Arg1173) stabilize inhibitors by electrostatic contributions. Residues in domain1 anchor inhibitors by forming the major hydrogen bonds. Inhibitor could express good inhibiting effect when it binds well with the key residues within the binding site. Charged substituents could form favorable interaction with Gln1113 in domain2 to further stabilize inhibitor. The presence of bulky and hydrophobic substituents interacting with residues in domain3 should improve binding affinity. Based on energetic and conformational analyses, the core fragment (benzodiazepinone with methylindazole) of CPI-637 is chosen to serve as the starting point for new inhibitor design by means of LUDI and

ration drug design approaches. In addition, two methods of toxicity test were performed to compare the oral bioavailability and drug promiscuity of the candidates. As a result, design-5 is the preferential molecule, as it could form more stable hydrogen bonds network and satisfy the rational range of cLogP. Our results may be effective for the design of new inhibitors and treatment for the mutations of p53 encoding gene related human cancers. Acknowledgement This work is supported by Natural Science Foundation of China [grant numbers 21273095, 21773084]. A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.bmc.2017.12.040. References 1. Xu M, Unzue A, Dong J, Spiliotopoulos D, Nevado C, Caflisch A. J Med Chem. 2016;59:1340–1349. 2. Popp TA, Tallant C, Rogers C, et al. J Med Chem. 2016;59:8889–8912. 3. Marmorstein R, Zhou MM. Cold Spring Harb Perspect Biol. 2014;6:a018762. 4. Unzue A, Zhao H, Lolli G, et al. J Med Chem. 2016;59:3087–3097. 5. Filippakopoulos P, Knapp S. Nat Rev Drug Discovery. 2014;13:337–356. 6. Sanchez R, Meslamani J, Zhou M-M. Biochim Biophys Acta. 2014;1839:676–685. 7. Magno A, Steiner S, Caflisch AJ. Chem Theory Comput. 2013;9:4225–4232. 8. Noguchi-Yachide T, Sakai T, Hashimoto Y, Yamaguchi T. Bioorg Med Chem. 2015;23:953–959. 9. Filippakopoulos P, Picaud S, Fedorov O, et al. Bioorg Med Chem. 2012;20:1878–1886. 10. Kuang M, Zhou J, Wang L, Liu Z, Guo J, Wu R. J Chem Inf Model. 2015;55:1926–1935. 11. Dancy BM, Cole PA. Chem Rev. 2015;115:2419–2452. 12. Ghosh S, Taylor A, Chin M, et al. J Biol Chem. 2016;291:13014–13027. 13. Taylor AM, Cote A, Hewitt MC, et al. ACS Med Chem Lett. 2016;7:531–536. 14. Crawford TD, Romero FA, Lai KW, et al. J Med Chem. 2016;59:10549–10563. 15. Hay DA, Fedorov O, Martin S, et al. J Am Chem Soc. 2014;136:9308–9319. 16. Studio D. Version 2.5. San Diego, CA, USA: Accelrys Inc.; 2009. 17. Case D, Darden T, Cheatham Iii T, et al. AMBER 11, 142. San Francisco: University of California; 2010. 18. Khoury GA, Thompson JP, Smadbeck J, Kieslich CA, Floudas CA. J Chem Theory Comput. 2013;9:5653–5674. 19. Khoury GA, Smadbeck J, Tamamis P, Vandris AC, Kieslich CA, Floudas CA. ACS Syn Biol. 2014;3:855–869. 20. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. J Chem Phys. 1983;79:926–935. 21. Shao J, Tanner SW, Thompson N, Cheatham TE. J Chem Theory Comput. 2007;3:2312–2334. 22. DeLano WL. PyMOL. San Carlos, CA: DeLano Scientific; 2002. 700. 23. Pettersen EF, Goddard TD, Huang CC, et al. J Chem Theory Comput. 2004;25:1605–1612. 24. Humphrey W, Dalke A, Schulten K. J. Mol. Graphics. 1996;14:33–38. 25. Sun H, Li Y, Tian S, Wang J, Hou T. PLoS Comput Biol. 2014;10:e1003729. 26. Sun H, Li Y, Shen M, et al. Phys Chem Chem Phys. 2014;16:22035–22045. 27. Chen F, Liu H, Sun H, et al. Phys Chem Chem Phys. 2016;18:22129–22139. 28. Xu L, Zhang Y, Zheng L, et al. J Med Chem. 2014;57:3737–3745. 29. Böhm H-J. J Comput Aided Mol Des. 1994;8:623–632. 30. Böhm H-J. J Comput Aided Mol Des. 1992;6:61–78. 31. Böhm H-J. J Comput Aided Mol Des. 1992;6:593–606. 32. Frisch M, Trucks G, Schlegel H, et al. Theor Chem Acc. 2008;120:215. 33. Becke AD. J Chem Phys. 1993;98:5648–5652. 34. Araújo JQ, de Mesquita Carneiro JW, de Araujo MT, Leite FHA, Taranto A. Bioorg Med Chem. 2008;16:5021–5029. 35. Ditchfield R, Hehre WJ, Pople JA. J Chem Phys. 1971;54:724–728. 36. Lee C, Yang W, Parr RG. Phys Rev B. 1988;37:785–789. 37. T.E.S.T. 3.0, available at EPA web site: http://www.epa.gov/. 38. Jarrahpour A, Motamedifar M, Zarei M, et al. Phosphorus. 2010;185:491–497. 39. Jarrahpour A, Fathi J, Mimouni M, et al. Med Chem Res. 2011;21:1984–1990. 40. Rooney TP, Filippakopoulos P, Fedorov O, et al. Angew Chem Int Ed. 2014;53:6126–6130.