Journal of Molecular Graphics and Modelling 83 (2018) 42e52 Contents lists available at ScienceDirect Journal of Molecular Graphics and Modelling journal homepage: www.elsevier.com/locate/JMGM Structural dynamics and quantum mechanical aspects of shikonin derivatives as CREBBP bromodomain inhibitors Sarmistha Mitra a, Raju Dash b, c, * a Department of Pharmacy, University of Chittagong, Chittagong, 4331, Bangladesh Molecular Modeling & Drug Design Laboratory (MMDDL), Pharmacology Research Division, Bangladesh Council of Scientific & Industrial Research (BCSIR), Chittagong, 4220, Bangladesh c Department of Biochemistry and Biotechnology, University of Science & Technology Chittagong, Chittagong, 4202, Bangladesh b a r t i c l e i n f o a b s t r a c t Article history: Received 29 November 2017 Received in revised form 16 April 2018 Accepted 23 April 2018 Available online 4 May 2018 The Proteins involved in the chemical modification of lysine residues in histone, is currently being excessively focused as the therapeutic target for the treatment of cell related diseases like cancer. Among these proteins, the epigenetic reader, CREB-binding protein (CREBBP) bromodomain is one of the most prominent targets for effective anticancer drug design, which is responsible for the reorganization of acetylated histone lysine residues. Therefore, this study employed an integrative approach of structure based drug design, in combination with Molecular Dynamics (MD) and QM/MM study to identify as well as to describe the binding mechanism of two shikonin derivatives, acetylshikonin and propionylshikonin as inhibitors of CREBBP bromodomain. Here induced fit docking strategy was employed to explore the important intrinsic interactions of ligands with CREBBP bromodomain, consistently molecular dynamics simulation with two different methods and binding energy calculations by MM-GBSA and MM-PBSA were adopted to determine the stability of intermolecular interactions between protein and ligands. The results showed that both these derivatives made direct contacts with the important conserved residues of the active site, where propionylshikonin demonstrated stronger binding and stability than acetylshikonin, according to molecular dynamics simulation and binding free energy calculations. Further, QM/MM energy calculation was employed to study the chemical reactivity of the propionylshikonin and also to describe the mechanism of non bonded interactions between the propionylshikonin and CREBBP bromodomain. Though this study demands in vitro and in vivo experiments to evaluate the efficiency of the compound, these insights would assist to design more potent CREBBP bromodomain inhibitor, guiding the site of modification of propionylshikonin moiety for designing selective inhibitors. © 2018 Elsevier Inc. All rights reserved. Keywords: CREBBP Bromodomain Molecular dynamics Propionylshikonin MM-PBSA MM-GBSA QM/MM 1. Introduction Chromatin, being a combination of histone protein and DNA, contains the entire genomic DNA in eukaryotic cells [1]. “Epigenetics” is a process which stands for heritable changes in gene expression but no amendment in DNA sequence [2]. The best known epigenetic processes are DNA methylation [3], histone protein's phosphorylation, acetylation, methylation, ubiquitination and sumoylation all of which can play a great role in cancer development [4]. Among them, lysine acetylation is one of the main * Corresponding author. Molecular Modeling and Drug Design Laboratory, Pharmacology Research Division, Bangladesh Council of Scientific & Industrial Research (BCSIR) Chittagong, 4220, Bangladesh. E-mail address: rajudash.bgctub@gmail.com (R. Dash). https://doi.org/10.1016/j.jmgm.2018.04.014 1093-3263/© 2018 Elsevier Inc. All rights reserved. modifications occurring in histone N-terminals tails, is known to be linked with activation of transcription through opening of chromatin architecture and with gene transcription [5]. The lysine acetyalation is established, modulated and interpreted by the several epigenetic enzymes and protein-protein interaction domains. Some enzymes are known as “writer”, catalyze the addition of post translational modifications. By the enzymes, called epigenetic “eraser”, post translational modification are removed. The domains that mediate post translational modifications recognition are called epigenetic “reader” [6e8]. Among the reader domains, Bromodomains (BRDs) have an important role in the targeting of chromatin-modifying enzymes and other proteins to specific sites, read the acetyl marks in histone tails and regulate the gene transcription [9]. These are proteins with acetyl-lysine binding modules; function as the principal mediators of molecular recognition S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 of acetylated chromatin [10]. “Bromodomain” is generally used to describe the bromodomain and extra-terminal (BET) family of human bromodomain proteins (BRD2, BRD3, BRD4). BRD-containing proteins have been reported to be involved in disease processes such as cancer, inflammation, cardiovascular diseases and viral replication [9,11]. From previous works, it was established that BET proteins are highly related to the genetical modification in cancer, directly regulating the expression of certain cancer-related genes, such as c-MYC. Preventing the BET proteins from binding at the MYC locus with BET inhibitors leads to a reduction in cell proliferation thus may be used as anti-cancer therapeutics [12e14]. Among the other bromodomains, one of the primary transcriptional regulators cAMP response element-binding protein (CREB) binding protein (CREBBP or CBP) bromodomain, which plays significant roles in the transcriptional activation of human cells [15]. The architecture of CREBBP protein consists of a HAT domain, a CREB binding domain, several zinc finger domains, a plant homology domain (PHD) and a bromodomain (BRD). Literature revealed that the bromodomain of CREBBP is one of the major considerations in anticancer drug target [16e19], involved several major biological functions like DNA replication and repair, genomic stability, cell cycle regulation and cell growth. Furthermore, various oncology-relevant transcription factors like cMYB, c-MYC and p53 are found to be amended by it [20e24]. Experiments done by Conery and Picaud group [15,25] showed that the bromodomain of CREBBP has better druggabilty than the other domains, inhibiting by small molecule leads to the enhancement the traditional chemotherapy by modulation intractable expression of oncogenic transcriptional factors [26]. In recent year, the shikonin derivatives, mainly isolated from the root extract of Lithospermum erythrorhizon have been reported to have anti-cancer, anti-inflammatory, wound healing, anti-viral and a wide range of biological effects [27,28]. Studies suggested that derivatives of shikonin are more effective than shikonin as well as less toxic and also have favorable pharmacokinetic and pharmacodynamic profiles in vivo [29]. Till to date, there are more published studies representing the anticancer activity of shikonin and its derivatives [30e35], however, little is known relating the activity of shikonin derivatives on bromodomains and other areas of study in epigenetics. The present study is concerned to evaluate CREBBP bromodomain inhibition by shikonin derivatives. For that, seventeen shikonin derivatives (Supplementary file 1) have been selected for virtual screening against CREBBP bromodomain using Glide extra precision docking and MM-GBSA binding energy calculation methods, where R-enantiomer of each compound has been considered (as described in Wang et al. [36]) for analysis. Among these derivatives, acetylshikonin and propionylshikonin have obtained the best binding energy (Supplementary file 1). The present study is therefore aimed to unveil the structural and molecular properties of acetylshikonin and propionylshikonin as inhibitors of CREBBP bromodomain. The investigations were accomplished by induced fit docking, classical MD simulations and two different binding energy calculations approaches. In addition, QM/MM calculation was also performed to describe the binding mode of the ligands at electronic level. 2. Methods 2.1. Dataset and ligand preparation The structures of two compounds were retrieved from PubChem database, i.e. acetylshikonin and propionylshikonin and then the 3D structures were built by using Ligand preparation wizard in Maestro 10.1 [39] with an OPLS_2005 force field. Their ionization states were €dinger Suite. generated at pH 7.0 ± 2.0 using Epik 2.2 in Schro 43 2.2. Protein preparation Three dimensional crystal structure of the bromodomain of human CREBBP (PDB id: 5I86, Chain: A,: Resolution: 1.05 Å, R-Value Free: 0.161 and R-Value Work: 0.139) was downloaded in pdb format from the protein data bank [37]. After that, the structure was prepared and refined using the Protein Preparation Wizard of €dinger-Maestro v10.1. Bond orders and charges were Schro imposed, hydrogens were combined to heavy atoms and all waters were removed. By optimizing the protein at neutral pH, readjustment of some thiol and hydroxyl groups, amide groups of asparagines, glutamines and imidazole ring of histidines, protonation states of histidines, aspartic acids and glutamic acids were done. Employing force field OPLS_2005, minimization was executed setting maximum heavy atom RMSD to 0.30 Å. 2.3. Induced fit docking Induced Fit Docking (IFD) was performed using the module €dinger-Maestro v10.1 [38]. In this Induced Fit Docking of Schro docking procedure, both ligands were docked into the target protein (PDB ID: 5I86) with a constrained minimization process and for generation of centroid of the residues 0.18 Å was chosen and the box size was automatically generated. After that, a soften potential glide docking was performed; in which, side chains were trimmed automatically based on B-factor and ligand and receptor with van der Waals scaling of 0.50 and 0.70, respectively. The generated number of poses were set to be 20. In the docking simulation, the closest residues to the ligand (within 5 Å of ligand pose) were kept flexible in prime refinement and during the process, the side chains were further optimized. Glide redocking process was further introduced for the ligand having the best pose within 30.0 kcal/mol. In the induced-fit receptor structure, the ligands were appropriately docked and for each output pose, the results yielded an IFD score. For further consideration, the pose having the lowest IFD score for each ligand was selected. 2.4. Molecular dynamics (MD) simulation MD simulations of the protein-ligand complexes obtained from induced fit docking; were performed by YASARA Dynamic software. Each complex was first cleaned and hydrogen bond network was optimized. A cubic simulation cell was generated with periodic boundary condition. Applying the AMBER14 [40,41] force field, the atoms of each complex were typed. Assignment of the parameters of each molecule was done through AutoSMILES [42] algorithms, where unknown organic molecules are fully automatically parameterized by the calculation of semi-empirical AM1 Mulliken point charges [42] with the COSMO solvation model, assigning of AM1BCC [43] atom and bond types, and also assigning GAFF [44] (General AMBER Force Field) atom types and remaining force field parameters. By using the transferable intermolecular potential3 points (TIP3P) water model the simulation box was solvated; maintaining the density of 0.997 g/L. During the solvation, pKa (acid dissociation constant) values of titratable amino acids of the protein were calculated. Each simulation system consisted of 62521 ± 10 atoms. An energy minimization protocol was used to remove the conformational stress, where the initial structure was energy minimized by steepest descent without electrostatic interaction. After that, the structure was subsequently relaxed by steepest descent minimization and subjected to full potential energy minimization for 5000 cycles, until convergence was reached. The energy was improved by less than 0.05 kJ/mol per atom during 5000 steps. MD simulations were carried out using the PME methods to describe long-range 44 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 electrostatic interactions at a cut off distance of 8 Å at physiological conditions (298 K, pH 7.4, 0.9% NaCl) [45]. A multiple time step algorithm together with a simulation time step interval of 2.50 fs was chosen [46]. At a constant pressure and Berendsen thermostat, MD simulations were performed for 100 ns long and MD trajectories were saved every 25 ps for further analysis. Using YASARA structure built in macros, VMD [47] and Bio3D [48] software, the resulting trajectories were subjected to analyze for the stability by various evaluative measures viz. RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuation), Radius of gyration and protein backbone comparisons. After that, MM-PBSA (Molecular mechanics-PoissoneBoltzmann Surface Area) binding free energies of all snapshots were calculated with YASARA software using the following formula, Binding Energy ¼ EpotRecept þ EsolvRecept þ EpotLigand þ EsolvLigand EpotComplex - EsolvComplex Here, YASARA built in macros was used to calculate MM-PBSA binding energy, using AMBER14 as a Force Field, where resulting more positive energies indicate the better binding. For calculating MM-GBSA binding free energy analysis, additional 30 ns MD simulations were performed with OPLS_3 force field by DESMOND software. The details of the method have been described in Supplementary file 2. 2.5. Prime MM-GBSA analysis Binding free energy was calculated to rescore and for choosing the top hits from the candidate ligands. In prime MM-GBSA method, the calculation of binding energy was performed by combining OPLSAA molecular mechanics energies (EMM), an SGB solvation model for polar solvation (GSGB) and a non-polar solvation term (GNP) composed of the non-polar solvent accessible surface area and van der Waals interactions [49]. Here, as the source in Prime MM-GBSA simulation, the glide pose viewer file of the best conformation was given. For modeling directionality of hydrogen-bond and p-stacking interactions, the dielectric solvent model VSGB 2.0 [50] was used to apply empirical corrections. Keeping the protein chain flexible [51e55], minimizing approach is applied as sampling methods. The analysis denotes more excellent binding by more negative binding energy. The overall free energy of binding: DGbind ¼ Gcomplex e (Gprotein þ Gligand), where G ¼ EMM þ GSGB þ GNP 2.6. QM/MM calculation In order to analyze the molecular interactions of ligand-protein complex at the electronic label, hybrid quantum mechanical/molecular mechanical (QM/MM) calculations had been employed. For QM/MM calculations, the QSITE program of Schrodinger suites was used. For the QM region, this module used JAGUAR [56] program and the IMPACT molecular modeling code for the MM region. As there were covalent connections between the QM and MM regions, frozen orbitals were used to mediate an interface between the two regions [57] during the calculation. The calculation of the interactions between QM and MM parts had been performed through van der Waals interactions among QM and MM atoms and the electrostatic interactions between MM point charges and the QM wave function. The MM region was treated with OPLS-AA force field [58], while for QM region density functional theory was used with the 6-31G*/LACVP* basis set, B3LYP density functional and “Ultrafine” SCF accuracy level (iacc ¼ 1, iacscf ¼ 2) [59]. All calculations were converged before reaching the limit; the highest number of iterations was set to 100 cycles and the root mean squared change in density matrix elements was less than the criterion of 5.0  10 6. 3. Results and discussion 3.1. Induced-fit docking analysis Many proteins undergo side chain or backbone movements, or both, upon ligand binding so the presumption of a rigid receptor always can seldom model the actual binding mechanism. Receptor alters its binding site due to these changes, so that is more closely conforms to the shape and binding mechanism of the ligand. This is one of the main complicating factors in structure-based drug design, which is often termed as “induced fit” [60]. To study the binding mechanism of the selected ligands, the IFD protocol was used. Bromodomains are composed with four a-helices (aZ, aA, aB, aC), constructed by ~110 residues, which are connected by interhelical loops, is also termed as bromodomain fold. Experimental studies reported that the N-acetyl group of lysine residues binds in most bromodomain by forming hydrogen bonds with ASN1168 residue of large hydrophobic pocket, which is attached by two loop regions named ZA and BC. Furthermore, the acetyl group of lysine residues is also seen to form water-mediated hydrogen bond with another conserved residue TYR1125, where the interaction is made through the phenolic hydroxyl group to acetyl carbonyl of lysine [26,61e64]. Other studies concluded that compound, which may act as an inhibitor is proved to have electrostatic attractions between the conserved ASN1168, TYR1125, PRO1110 side chains and also hydrophobic interactions with ILE1122, LEU1120 and VAL1174 residues are critical for inhibitor binding and inhibition [24,65,66]. From the docking results (Table 1), it was found that compound acetylshikonin, formed hydrogen bonds with the ASN1168 and PRO1110 residues in the active site of the protein. In these two bonds, the corresponding distances are 2.29 Å and 2.23 Å, respectively (Table 2), where both of these hydrogen bonds are contributed by naphthazarin ring of acetylshikonin. As shown in Fig. 1, this ring was also formed hydrophobic interactions with this ligand by forming pi-alkyl bonds with VAL1115 (4.21 Å), VAL1174 (4.42 Å), ALA1164 (5.33 Å), VAL1115 (4.31 Å), ALA1164 (4.58 Å) residues on the ZA loop. Two other alkyl hydrophobic interactions were also seen with the side chain of LEU1120 at the ZA loop with 4.47 Å and 3.96 Å bond distances respectively. While another ligand, propionylshikonin, here the naphthazarin ring was also found to have H-bonds with TYR1125 with bond distance of 2.55 Å; VAL1115 with 2.79 Å; ASN1168 with 2.10 Å; PRO1110 with 2.76 Å, GLN1113 with 2.77 Å and MET1160 with 2.18 Å (Fig. 2). Furthermore, another important residue ARG1173 formed hydrogen bond with oxygen atom of propionyl moiety having bond length of 2.72 Å. The ligand was also interacted with ILE1122 (3.96 Å), ILE1122 (4.22 Å), PRO1110 (4.49 Å) by forming alkyl hydrophobic bonds and with TYR1125 (4.39 Å), TYR1167 (5.21 Å), TYR1167 (4.21 Å), PRO1110 (5.45 Å), VAL1115 (4.66 Å), VAL1174 (4.85 Å), VAL1115 (4.62 Å) residues it formed pi-alkyl bonds. The side chain TYR1125 was also involved with two hydrophobic pi-pi interactions to ligand with the bond distance of 5.41 Å and 5.89 Å, respectively. From the above study, it is found that acetylshikonin binds with CREBBP bromodomain through interacting with ASN1168, PRO1110 residues by means of hydrogen bonds and the hydrophobic bond with VAL1174 residue. In case of propionylshikonin, the ligand formed polar interaction with ASN1168, TYR1125 residues and hydrophobic bonds with ILE1122 and VAL1174. These finding denotes 45 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 Table 1 Docking simulations analysis and protein-ligand binding free energy values of CREBBP bromodomain and the shikonin derivatives. Compound name Docking score (kcal/ mol) Acetylshikonin Propionylshikonin a b 8.14 9.85 IFD score (kcal/ mol) 273.02 273.48 Glide score (kcal/ mol) Glide emodel (kcal/ mol) 8.15 9.86 59.62 88.33 MM/PBSA Binding Energya (kcal/ MM/GBSA Binding Energyb (kcal/ mol) mol) 9.80 27.65 54.92 59.58 Average binding free energy calculated in 100ns MD simulation using AMBER14 force field. Average binding free energy calculated in 30ns MD simulation using OPLS_3 force field. Table 2 Binding site analysis of acetylshikonin and propionylshikonin with CREBBP Bromodomain. Compound name Hydrogen bonds (Distance in Å) Hydrophobic bonds (Distance in Å) Alkyl- Alkyl Pi-alkyl Pi-Pi Acetylshikonin ASN1168(2.29) PRO1110(2.23) LEU1120(4.47) PRO1120(5.25) LEU1120(3.96) Propionylshikonin VAL1115(2.79) TYR1125 (2.55) ASN1168 (2.10) ARG1173 (2.72) PRO1110 (2.76) GLN1113(2.77) MET1160(2.18) ILE1122(3.96) ILE1122(4.22) PRO1110(4.49) VAL1115(4.21) VAL1174(4.42) ALA1164(5.33) VAL1115(4.31) ALA1164(4.58) TYR1125(4.39) TYR1167(5.21) TYR1167(4.21) PRO1110(5.45) VAL1115(4.66) VAL1174(4.85) VAL1115(4.62) TYR1125(5.41) TYR1125(5.89) Fig. 1. Molecular interactions between acetylshikonin and CREBBP bromodomain. Here, the three dimensional representation is in hydrophobic surface view, and in 2D view, green dot lines represent hydrogen bond and pink for pi-alkyl interaction. that both the ligands bind to the active site of bromodomain and have the chance to act as inhibitors. 3.2. Molecular dynamics (MD) simulation In order to understand and evaluate the interactions of the ligand molecules with the CREBBP bromodomain, we carried out 100 ns of MD simulation for each protein-ligand complex and ligand free protein. Two main parameters, which were analyzed throughout the simulations, were Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuations (RMSF). RMSD value of the protein backbone of the complexes were calculated individually and plotted by comparing with the initial trajectory snapshot of protein structure. In the plot, the smaller difference of RMSDs represents the higher stability of protein conformation. As shown in Fig. 3a, the RMSD values for the protein backbone of acetylshikonin complex were found to fluctuate around 3 Å, while 2 Å for protein backbone of propionylshikonin. Interestingly, both these complexes were stabled during the simulations, after a structural drifting at 22 ns In a comparative manner, the protein backbone of acetylshikonin was less rigid than propionylshikonin, indicating the greater flexibility. While the ligand free CREBBP showed larger conformational changes; represented RMSD around 1.5 Å for 20 ns, however afterward significant fluctuations were observed. The major structural drifting was observed at 41 ns, which caused the RMSD raise to 3.6 Å and remained stable subsequently. In contrast, RMSD values of two ligands are also plotted together and from the Fig. 3b, it is clearly seen that the RMSD value of acetylshikonin is much higher than the RMSD value of 46 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 Fig. 2. Molecular interactions between propionylshikonin with CREBBP bromodomain. Here, the three dimensional representation is in hydrophobic surface view, and in 2D view, green dot lines represent hydrogen bond and pink for pi-alkyl interaction. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 3. Time evolution of root-mean-square deviations (RMSDs) of backbone atoms (C, Ca, and N) for a) protein and b) ligand of each docked complex. Here, orange and ash colors denote acetylshikonin, propionylshikonin complexes, respectively, while green line indicates apo-CREBBP bromodomain. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) propionylshikonin. Though the RMSD of acetylshikonin was ranged from 0.3 to 2 Å, this compound showed different level of fluctuations in several intervals. The Fig. 3b clearly describes that at initial position, the RMSD of acetylshikonin fluctuated to 2 Å and remained stable for 22 ns, after the curve sudden fall to 0.6 Å and continued to around 1.5 Å till to end of the simulation with several large fluctuations. On the contrary, RMSD of another ligand, propionylshikonin demonstrates stability in the active site during the simulation, as the level of RMSD within the range from 0.5 to 1.8 Å. For describing the local changes along the protein chain, RMSF (Root Mean Square Fluctuation) is used in this study (Fig. 4a). On this plot, peaks demonstrate the areas of the protein that fluctuated most in the entire simulation period. The overall fluctuations of the RMSF of the ligands were found from a range of 1e8 Å throughout Fig. 4. a) Ligand induced thermodynamic behaviour of protein during the simulations. The plots described structural changes of protein by means of Root Means Square Fluctuations (RMSF). Here, green, orange and ash colors denote apo-CREBBP bromodomain, acetylshikonin, propionylshikonin complexes respectively. The baseline below describes the secondary structure contents (Red color for alpha helix and green for loops) of CREBBP bromodoamin, calculated by VMD software, where the alpha helixes were seen in the regions of 1087e1094, 1096e1099, 1026e1028, 1135e1143, 1150e1167, 1174e1189, 1191e1196 respectively. b) Comparative Root Mean Square Fluctuations (RMSF) of the active site residues of acetylshikonin, propionylshikonin and apo-CREBBP bromodomain. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) the simulation. In all cases, the loop regions were found to show more fluctuation than the area of alpha helix of the proteins. Highest RMSF was observed at the regions of ZA loop, ALA1110 to ASP1120, which was around 4 Å. At this region, acetylshikonin showed large RMSF difference than the others, which may due to the interactions of ligand to active site residues, located in the loop. S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 To justify this statement, RMSF of ligand binding residues were further analyzed and rendered in Fig. 4b in a comparative manner. As can be seen in Fig. 4b, acetylshikonin interacting residues produced higher RMSF compared with other. Moreover, three main residues including Pro1110, Gln1113, Val1115 and TYR1125 showed high fluctuations, which confirms the flexibility of the loop ZA and it may undergoes different conformational changes. For better understanding, RMSD and Rg (Radius of Gyration) of the two loops that connected with the binding site, ZA and BC have been calculated and illustrated in Fig. 5a&b. The Rg value is an indicator of protein flexibility describes the compactness, where high value indicates loose packing. According to the Fig. 5a, the Rg of ZA loop was seen to fluctuate more in acetylshikonin complex, while it was 47 stabilized in propionylshikonin complex. In case of BC loop (Fig. 5b), both apo and propionylshikonin complex were seen to stable over the time, however, the fluctuations in Rg are more greater for acetylshikonin. The RMSD analysis for both these loops also supports the result from Rg, where the results of RMSD are directly correlated with the results of RMSD calculated from total protein. This observation also clearly concluded that flexibility of these loops directly modulate the overall flexibility of protein. Furthermore, significant differences in the fluctuations of ZA loop were also observed by the closer view of principal component analysis (PCA), where flexibility of this loop in acetylshikonin complex was associated with the changes in pocket shape and volume of CREBBP bromodomain. From the previous analysis of ligand fluctuation Fig. 5. Loop conformational changes, i.e. a) ZA loop (PRO1110 to TYR1125), b) BC loop (ASN1168 to SER1172), by means of Radius of gyration (Rg) and RMSD. The green line represents the apo-CREBBP, while orange and ash colors denote acetylshikonin, propionylshikonin complexes respectively. c) The first principal component of dynamics from bromodomain simulations calculated by Bio3D and imaged by VMD software. From left to right, the apo-CREBBP, acetylshikonin and propionylshikonin are shown in tube, colored from red to blue indicating increasing timestep. Here enlarge tube represent the flexibility. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 48 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 (Fig. 3b), the flexibility of acetylshikonin was seen to change dramatically, which may due to the rotation of acetyl moiety and this rotation may cause the conformational changes of the loops because of the molecular interactions between the ligand and residues of these loops. Thus, higher RMSD and RMSF were obtained for acetylshikonin. The Fig. 6 showed the number of H-bonds formed between the ligand and residues of ligand binding site of protein for each complex. As shown in Fig. 6 a&b, maximum hydrogen bonds were formed by propionylshikonin during the simulation, while acetylshikonin showed minimum interactions (Table 2). For the detailed understanding, we calculated percentage of H-bond occupancy, showed in Fig. 6b. Fig. 6c shows the percentage of H-bond occupancy, where each ligand involved acting as H-bond donor and the residues in the catalytic site as H-bond acceptor. As shown in Fig. 6c, acetylshikonin formed several H- bonds with ASN1168, TYR1125, GLN1113 residues during the simulation. Throughout the H-bond formation, acetylshikonin contributed as an acceptor and showed less than 15% occupancy with ASN1168, TYR1125, GLN1113 residues. However, as a donor, acetylshikonin shows less percentage of H-bond occupancy through the simulation, only 1.13% for GLN1113. Compared with acetylshikonin, propionylshikonin formed H-bonds with ASN1168, TYR1125 residues of the active site, where this ligand showed 40.88% of H-bond occupancy with ASN1168 as H-bond acceptor. In addition, water mediated hydrogen bonds were also observed for both complexes in 30 ns simulations by OPLS_3 force field (Figure S3, Supplementary file 2). In both complexes, active site residues including PRO1110, GLN1113, TYR1125, ARG1173, VAL1174 and ASN1163 were seen to form water mediated hydrogen bonds with the ligands. It is noteworthy to mention that the TYR1125 residue which formed water mediated hydrogen bonds with the acetyl group of lysine residue [26,61e64] was also formed similar interactions with both ligands, where propionylshikonin showed more interactions than acetylshikonin. In acetylshikonin complex, MET1133 and ASN1163 residues formed water mediated hydrogen bonds more than 20% and 30% of the simulation time. ARG1173 residue was also formed water bridges with both ligands, however, less than 20% of the simulation time. In addition, the results of RMSD, RMSF and molecular interactions were also calculated (Data are shown in Supplementary File 2), where the results were found to be consistent with the results of 100 ns simulations. The outcomes from both these simulations clearly conclude that propionylshikonin formed more stable complex with CREBBP bromodomain, however, binding energy is one of the main concern in protein ligand complex, which was analyzed and discussed in further section. 3.3. Binding free energy evaluation Fig. 6. Total number of hydrogen bonds formed between the protein and ligand in complex state, during the simulation. Here two figures at the upper position represent the number of intermolecular H-bond formed in a) propionylshikonin, b) acetylshikonin complexes respectively. And, the rest of the figure in the bottom depicts the percentage occupancy of hydrogen bond contacts formed between the ligands and active site residues. Figure c) shows the time evolution of the occupancy of the listed hydrogen bonds, where each ligand acts as acceptor. As stated before in methods sections, MD simulation of each complex was run for 100 ns in YASARA software, with the force field of AMBER14. During simulations, snapshots were saved every 25 ps, therefore total 4000 snapshots were generated in 100 ns simulation and all were considered for MM-PBSA calculation. The results of MM-PBSA binding energy for two complexes are plotted against time and shown in Fig. 7a. According to the analysis, propionylshikonin-CREBBP complex obtained highest average positive binding energy of 27.65 kcal/mol, while acetylshikonin showed the less binding energy of 9.80 kcal/mol (Table 1), which indicates weak binding. It should be noted that more positive MMPBSA binding energy designated as the better binding of protein ligand complex. In the same way, in MM-GBSA calculation, propionylshikonin showed tighter binding to the target structure compared to acetylshikonin, where the binding energy was 59.58 kcal/mol. However, acetylshikonin demonstrated less favored binding with the CREBBP bromodomain ( 54.92 kcal/mol). In MM-GBSA analysis, more negative binding energy represents strong binding. As can be seen in Fig. 7, the binding energy of acetylshikonin was fluctuated more in both different simulations (MM-PBSA and MMGBSA), while in both cases the binding energy of propionylshikonin was in steady level. It is further noteworthy to state that the acetylshikonin showed high ligand RMSD, which means that acetylshikonin had undergone to different conformations in the binding site of target protein during simulation that caused the formation and breakdown of molecular interactions with the cavity residues and thus resulted in fluctuation of binding energy. The results of binding energy calculation are also consistent with the ligand RMSD profile. 49 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 Fig. 7. : Binding free energy calculations of acetylshikonin and propionylshikonin complexes during the molecular dynamics simulation. Binding energy was evaluated with two different methods, a) MM-PBSA and b) MM-GBSA, respectively. Here, orange and ash colors denote acetylshikonin, propionylshikonin complexes respectively. In MM-PBSA analysis, binding energy calculated from trajectory of 100 ns molecular dynamics simulation with AMBER14 force field, where positive energy indicates better binding. In contrast, MM-GBSA binding energy calculated from the snapshots of 30 ns molecular dynamics simulation with OPLS_3 force field, where more negative binding energy denotes better binding. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 3.4. QM/MM calculation Concurrent with the investigation, QM/MM calculation has been performed on propionylshikonin and propionylshikonin-CREBBP bromodomain complex, as it demonstrated the better binding free energy in both MM-GBSA and MM-PBSA analysis. Application of Quantum Mechanics (QM) may provide a deeper understanding of molecular systems, where the basic physical quantities are described as a function of the electronic structure. With the development of advanced computational facilities (both in terms of software and hardware), QM based methods have become the center of attraction in the study of biomolecules. However, applications of QM methods are limited to small sized systems but the development of polarization force fields [67,68], variants of second generation pair-additive force fields [69], have improved the results of MM approach. In this study, QM/MM calculations have been executed on three different stages of propionylshikonin, i) QM/MM on propionylshikonin in TIP3P solvent system (Fig. 8a), where ligand was treated with QM and solvent by MM; ii) QM/MM on propionylshikonin in the binding with CREBBP protein in the TIP3P solvent system (Fig. 8b), here only ligand was subjected in QM region and rest of the portions were considered in MM region; iii) QM/MM on propionylshikonin in binding with CREBBP protein, considering VAL1115, TYR1125, ASN1168, ARG1173, PRO1110, GLN1113, MET1160 residues and the propionylshikonin in QM region (Fig. 10a) while resting of the portion in MM. During the QM optimization, the calculation of frontier molecular orbitals energies of propionylshikonin by means of HOMO and LUMO was carried out to elucidate the conformational energy of the ligand at the higher level of theory and also protein reorganization ability. Table 3 describes the HOMO and LUMO energies of the ligand in first two different stages; i.e. in solvent and in binding to the protein, the band gap between HOMO and LUMO was also calculated. Interestingly, the ligand showed larger band gap of 0.10552 eV in solvent system, however, the gap was seen to decrease when the ligand binds with the protein, as described in Table 3. The increase of band gap between HOMO and LUMO described the less chemical reactivity of a molecule; while the lower band gap denotes higher chemical reactivity, which is mostly considered in drug designing applications [70e72]. Furthermore, Fig. 9 represents the molecular orbital maps of HOMO and LUMO, describing the sites of ligand as electron donor and acceptor. Fig. 9 described the HOMO maps of propionylshikonin, where the ring containing hydroxyl group of naphthazarin part of propionylshikonin represented the most electron-sufficient area and also in the heteroatom as oxygen. In contrast, the electron acceptors regions are defined by LUMO maps indicate the electron-transfer ability of ketone groups of naphthazarin part mostly. From Table 3 HOMO and LUMO energies of the orbital, calculated form QM/MM simulation of two different systems. System Propionylshikonin Propionylshikonin - CREBBP HOMO (eV) LUMO (eV) 0.20501 0.23789 0.09949 0.13815 Band Gap (eV) 0.10552 0.09974 Fig. 8. Solvated system of propionylshikonin in two different stages, i.e. a) propionylshikonin, b) propionylshikonin in binding with CREBBP domain. Here, transferable intermolecular potential3 points (TIP3P) water model has been used for building solvation system for each stage. 50 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 Fig. 9. The geometries of propionylshikonin after QM/MM optimization. a) The electron density distributions of HOMO, and LUMO for propionylshikonin, and b) The electron density distributions of HOMO, and LUMO for propionylshikonin in binding with CREBBP bromodomain. The negative energy is represented by blue color, while red color for positive energy. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) induced fit docking, the maximum number of electrostatic interactions between protein and ligand were formed by this naphthazarin part and in both MD simulations similar phenomenon was also observed. These findings are in agreement with the results from molecular docking and MD simulations, describe how propionylshikonin contributes in binding to CREBBP bromodomain. To gain more structural insights, we subjected the propionylshikonin-CREBBP complex with solvent system for QM optimization, including the VAL1115, TYR1125, ASN1168, ARG1173, PRO1110, GLN1113, MET1160 residues along with the ligands (Fig. 10a). The final optimized complex was obtained and all possible hydrogen bonds between the residues have been analyzed through. It has been revealed that ASN1168 residue, which is mostly contributed in lignad binding in docking and MD simulations, also formed hydrogen bond after QM optimization with a distance of 2.21 Å (Fig. 10b). Another important residue TYR1125 formed strong hydrogen bonds with the oxygen atom of ketone group of naphthazarin ring with corresponding distance of 1.93 Å MET1160 residue in the active site represented medium hydrogen bonds with the ligand, where the respective distance was 2.33 Å, while PRO1110 formed two weak hydrogen bonds with 2.56 Å and 3.16 Å bond distances, respectively. Furthermore, the ketone group of propionyl Fig. 10. Propionylshikonin-CREBBP bromodomain complex for QM/MM calculation, where a) the QM part was consisted with the stick represented residues and the right side, Figure b) demonstrates the binding interactions of propionylshikonin-CREBBP bromodomain after QM/MM calculation. In Figure b), green dot lines describes H-bond between the only residues treated in QM region. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 moiety showed hydrogen bond with the guanidinium of the ARG1173 side chain with a bond distances of 2.08 Å. It should be noted that the selectivity CREBBP bromodomain inhibitor is served by the polar interaction with the ARG1173 side chain [65,66,73,74]. Though propionylshikonin contributed small portion of electrostatic interaction with the ARG1173 side chain throughout MD simulations, according to the QM/MM calculation, the modification of functional group at the same position of propionyl moiety could provide strong electrostatic interaction with this side chain and thereby selectivity can be achieved. 4. Conclusion In this study, an attempt has been made to study the binding potentiality of shikonin derivatives, as they are less toxic than shikonin. Using comprehensive molecular dynamics simulations and binding energy calculations, the study concluded that propionylshikonin has the strong binding to the CREBBP bromodomain. Furthermore, propionylshikonin is more reactive when binds with bromodomain, which is explained by the smaller band gap between the HOMO and LUMO. In addition, the QM/MM calculation revealed, modifications of propionyl moiety of propionylshikonin could be enabled to obtain selective CREBBP bromodomain inhibition. Therefore, these findings are valuable for the designing of target specific new CREBBP bromodomain inhibitor from propionylshikonin and shikonin derivatives and will provide a suitable starting point for further development as well as detailed in vitro and in vivo analysis. Acknowledgements We the authors are heartily grateful to Mr. S M Zahid Hosen, Scientific Officer, Molecular Modeling and Drug Design Laboratory, Pharmacology Research Division, BCSIR laboratories Chittagong, for providing the facilities to study Molecular Dynamics Simulations. Conflicts of interest All authors declare that they have no competing interests. Appendix A. Supplementary data Supplementary data related to this article can be found at https://doi.org/10.1016/j.jmgm.2018.04.014. References [1] D.E. Schones, K. Zhao, Genome-wide approaches to studying chromatin modifications, Nat. Rev. Genet. 9 (2008) 179. [2] E. Jablonka, M.J. Lamb, Epigenetic Inheritance and Evolution: the Lamarckian Dimension, Oxford University Press on Demand, 1999. [3] B. Weinhold, Epigenetics: the science of change, Environ. Health Perspect. 114 (2006) A160. [4] A.J. Bannister, T. Kouzarides, Regulation of chromatin by histone modifications, Cell Res. 21 (2011) 381. [5] T. Kouzarides, Chromatin modifications and their function, Cell 128 (2007) 693e705. [6] E. Ferri, C. Petosa, C.E. McKenna, Bromodomains: structure, function and pharmacology of inhibition, Biochem. Pharmacol. 106 (2016) 1e18. [7] A. Sadakierska-Chudy, M. Filip, A comprehensive view of the epigenetic landscape. Part II: histone post-translational modification, nucleosome level, and chromatin regulation by ncRNAs, Neurotox. Res. 27 (2015) 172e197. [8] V. Pande, Understanding the complexity of epigenetic target space: miniperspective, J. Med. Chem. 59 (2016) 1299e1307. rez-Salvia, M. Esteller, Bromodomain inhibitors and cancer therapy: [9] M. Pe from structures to applications, Epigenetics 12 (2017) 323e339. [10] S. Muller, P. Filippakopoulos, S. Knapp, Bromodomains as therapeutic targets, Expet Rev. Mol. Med. (2011) 13. [11] P. Filippakopoulos, S. Knapp, Targeting bromodomains: epigenetic readers of lysine acetylation, Nat. Rev. Drug Discov. 13 (2014) 337. 51 [12] L.-l. Fu, M. Tian, X. Li, J.-j. Li, J. Huang, L. Ouyang, et al., Inhibition of BET bromodomains as a therapeutic strategy for cancer drug discovery, Oncotarget 6 (2015) 5501. [13] D.M. Miller, S.D. Thomas, A. Islam, D. Muench, K. Sedoris, c-Myc and cancer metabolism, AACR 18 (2012) 5546e5553. [14] L. Ellis, P.W. Atadja, R.W. Johnstone, Epigenetics in cancer: targeting chromatin modifications, Mol. Canc. Therapeut. 8 (2009) 1409e1420. [15] S. Picaud, O. Fedorov, A. Thanasopoulou, K. Leonards, K. Jones, J. Meier, et al., Generation of a selective small molecule inhibitor of the CBP/p300 bromodomain for leukemia therapy, Canc. Res. 75 (2015) 5106e5119. [16] R.H. Goodman, S. Smolik, CBP/p300 in cell growth, transformation, and development, Genes Dev. 14 (2000) 1553e1577. [17] R.H. Giles, D.J. Peters, M.H. Breuning, Conjunction dysfunction: CBP/p300 in human disease, Trends Genet. 14 (1998) 178e183. [18] Y. Li, H.-X. Yang, R.-Z. Luo, Y. Zhang, M. Li, X. Wang, et al., High expression of p300 has an unfavorable impact on survival in resectable esophageal squamous cell carcinoma, Ann. Thorac. Surg. 91 (2011) 1531e1538. [19] M. Li, R.-Z. Luo, J.-W. Chen, Y. Cao, J.-B. Lu, J.-H. He, et al., High expression of transcriptional coactivator p300 correlates with aggressive features and poor prognosis of hepatocellular carcinoma, J. Transl. Med. 9 (2011) 1. [20] Q. Jin, L.R. Yu, L. Wang, Z. Zhang, L.H. Kasper, J.E. Lee, et al., Distinct roles of GCN5/PCAF-mediated H3K9ac and CBP/p300-mediated H3K18/27ac in nuclear receptor transactivation, EMBO J. 30 (2011) 249e262. [21] C. Das, M.S. Lucia, K.C. Hansen, J.K. Tyler, CBP/p300-mediated acetylation of histone H3 on lysine 56, Nature 459 (2009) 113e117. [22] A. Ito, C.H. Lai, X. Zhao, S.i. Saito, M.H. Hamilton, E. Appella, et al., p300/CBPmediated p53 acetylation is commonly induced by p53-activating agents and inhibited by MDM2, EMBO J. 20 (2001) 1331e1340. [23] A. Hammitzsch, C. Tallant, O. Fedorov, A. O'Mahony, P.E. Brennan, D.A. Hay, et al., CBP30, a selective CBP/p300 bromodomain inhibitor, suppresses human Th17 responses, Proc. Natl. Acad. Sci. Unit. States Am. 112 (2015) 10768e10773. ^ te , M.C. Hewitt, R. Pastor, Y. Leblanc, C.G. Nasveschuk, et al., [24] A.M. Taylor, A. Co Fragment-based discovery of a selective and cell-active benzodiazepinone CBP/EP300 bromodomain inhibitor (CPI-637), ACS Med. Chem. Lett. (2016). [25] A.R. Conery, R.C. Centore, A. Neiss, P.J. Keller, S. Joshi, K.L. Spillane, et al., Bromodomain inhibition of the transcriptional coactivators CBP/EP300 as a therapeutic strategy to target the IRF4 network in multiple myeloma, Elife 5 (2016), e10483. [26] L.R. Vidler, N. Brown, S. Knapp, S. Hoelder, Druggability analysis and structural classification of bromodomain acetyl-lysine binding sites, J. Med. Chem. 55 (2012) 7346e7359. [27] Y. Fujita, Y. Hara, C. Suga, T. Morimoto, Production of shikonin derivatives by cell suspension cultures of Lithospermum erythrorhizon, Plant Cell Rep. 1 (1981) 61e63. [28] X. Chen, L. Yang, J.J. Oppenheim, O. Howard, Cellular pharmacology studies of shikonin derivatives, Phytother Res. 16 (2002) 199e209. [29] I. Andujar, J.L. Rios, R.M. Giner, M.C. Recio, Pharmacological properties of shikonin - a review of literature since 2002, Planta Med. 79 (2013) 1685e1697. [30] J. Chen, J. Xie, Z. Jiang, B. Wang, Y. Wang, X. Hu, Shikonin and its analogs inhibit cancer cell glycolysis by targeting tumor pyruvate kinase-M2, Oncogene 30 (2011) 4297e4306. [31] G. He, G. He, R. Zhou, Z. Pi, T. Zhu, L. Jiang, et al., Enhancement of cisplatininduced colon cancer cells apoptosis by shikonin, a natural inducer of ROS in vitro and in vivo, Biochem. Biophys. Res. Commun. 469 (2016) 1075e1082. [32] S.H. Kim, I.C. Kang, T.J. Yoon, Y.M. Park, K.S. Kang, G.Y. Song, et al., Antitumor activities of a newly synthesized shikonin derivative, 2-hyim-DMNQ-S-33, Canc. Lett. 172 (2001) 171e175. [33] Y. Masuda, A. Nishida, K. Hori, T. Hirabayashi, S. Kajimoto, S. Nakajo, et al., Beta-hydroxyisovalerylshikonin induces apoptosis in human leukemia cells by inhibiting the activity of a polo-like kinase 1 (PLK1), Oncogene 22 (2003) 1012e1023. [34] U. Sankawa, Y. Ebizuka, T. Miyazaki, Y. Isomura, H. Otsuka, S. Shibata, et al., Antitumor activity of shikonin and its derivatives, Chem. Pharm. Bull. (Tokyo) 25 (1977) 2392e2395. [35] H. Wu, J. Xie, Q. Pan, B. Wang, D. Hu, X. Hu, Anticancer agent shikonin is an incompetent inducer of cancer drug resistance, PLoS One 8 (2013), e52706. [36] R. Wang, R. Yin, W. Zhou, D. Xu, S. Li, Shikonin and its derivatives: a patent review, Expert Opin. Ther. Pat. 22 (2012) 977e997. [37] H.M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T.N. Bhat, H. Weissig, et al., The protein data bank, Nucleic Acids Res. 28 (2000) 235e242. [38] T.N. Doman, S.L. McGovern, B.J. Witherbee, T.P. Kasten, R. Kurumbail, W.C. Stallings, et al., Molecular docking and high-throughput screening for novel inhibitors of protein tyrosine phosphatase-1B, J. Med. Chem. 45 (2002) 2213e2221. €dinger, Induced Fit Docking Protocol; glide Version 5.8, Prime Version [39] S. Schro €dinger, LLC, New York, 2012. 3.1, Schro [40] C.J. Dickson, B.D. Madej, Å.A. Skjevik, R.M. Betz, K. Teigen, I.R. Gould, et al., Lipid14: the amber lipid force field, J. Chem. Theor. Comput. 10 (2014) 865e879. [41] J.A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K.E. Hauser, C. Simmerling, ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB, J. Chem. Theor. Comput. 11 (2015) 3696e3713. 52 S. Mitra, R. Dash / Journal of Molecular Graphics and Modelling 83 (2018) 42e52 [42] J.J. Stewart, MOPAC: a semiempirical molecular orbital program, J. Comput. Aided Mol. Des. 4 (1990) 1e103. [43] A. Jakalian, D.B. Jack, C.I. Bayly, Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation, J. Comput. Chem. 23 (2002) 1623e1641. [44] J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, Development and testing of a general amber force field, J. Comput. Chem. 25 (2004) 1157e1174. [45] E. Krieger, J.E. Nielsen, C.A. Spronk, G. Vriend, Fast empirical pKa prediction by Ewald summation, J. Mol. Graph. Model. 25 (2006) 481e486. [46] E. Krieger, G. Vriend, New ways to boost molecular dynamics simulations, J. Computional. Chem 36 (2015) 996e1007. [47] W. Humphrey, A. Dalke, K. Schulten, VMD: visual molecular dynamics, J. Mol. Graph. 14 (1996) 33e38. [48] B.J. Grant, A.P. Rodrigues, K.M. ElSawy, J.A. McCammon, L.S. Caves, Bio3d: an R package for the comparative analysis of protein structures, Bioinformatics 22 (2006) 2695e2696. [49] B. Vijayakumar, A. Umamaheswari, A. Puratchikody, D. Velmurugan, Selection of an improved HDAC8 inhibitor through structure-based drug design, Bioinformation 7 (2011) 134e141. [50] J. Li, R. Abel, K. Zhu, Y. Cao, S. Zhao, R.A. Friesner, The VSGB 2.0 model: a next generation energy model for high resolution protein structure modeling, Proteins 79 (2011) 2794e2812. [51] F. Chen, H. Liu, H. Sun, P. Pan, Y. Li, D. Li, et al., Assessing the performance of the MM/PBSA and MM/GBSA methods. 6. Capability to predict proteineprotein binding free energies and re-rank binding poses generated by proteineprotein docking, Phys. Chem. Chem. Phys. 18 (2016) 22129e22139. [52] L. Xu, H. Sun, Y. Li, J. Wang, T. Hou, Assessing the performance of MM/PBSA and MM/GBSA methods. 3. The impact of force fields and ligand charge models, J. Phys. Chem. B 117 (2013) 8408e8421. [53] H. Sun, Y. Li, S. Tian, L. Xu, T. Hou, Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set, Phys. Chem. Chem. Phys. 16 (2014) 16719e16729. [54] H. Sun, Y. Li, M. Shen, S. Tian, L. Xu, P. Pan, et al., Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring, Phys. Chem. Chem. Phys. 16 (2014) 22035e22045. [55] T. Hou, N. Li, Y. Li, W. Wang, Characterization of domainepeptide interaction interface: prediction of SH3 domain-mediated proteineprotein interaction network in yeast by generic structure-based models, J. Proteome Res. 11 (2012) 2982e2995. [56] K. Senthilkumar, J.I. Mujika, K.E. Ranaghan, F.R. Manby, A.J. Mulholland, J.N. Harvey, Analysis of polarization in QM/MM modelling of biologically relevant hydrogen bonds, J. R. Soc. Interface 5 (2008) 207e216. [57] R.B. Murphy, D.M. Philipp, R.A. Friesner, A mixed quantum mechanics/molecular mechanics (QM/MM) method for large-scale modeling of chemistry in protein environments, J. Comput. Chem. 21 (2000) 1442e1457. [58] 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. [59] K.D. Singh, K. Muthusamy, Molecular modeling, quantum polarized ligand docking and structure-based 3D-QSAR analysis of the imidazole series as dual AT1 and ETA receptor antagonists, Acta Pharmacol. Sin. 34 (2013) 1592e1606. [60] W. Sherman, T. Day, M.P. Jacobson, R.A. Friesner, R. Farid, Novel procedure for modeling ligand/receptor induced fit effects, J. Med. Chem. 49 (2006) 534e553. [61] P. Filippakopoulos, S. Picaud, M. Mangos, T. Keates, J.-P. Lambert, D. BarsyteLovejoy, et al., Histone recognition and large-scale structural analysis of the human bromodomain family, Cell 149 (2012) 214e231. [62] P. Filippakopoulos, S. Knapp, The bromodomain interaction module, FEBS (Fed. Eur. Biochem. Soc.) Lett. 586 (2012) 2692e2704. [63] R. Sanchez, M.-M. Zhou, The role of human bromodomains in chromatin biology and gene transcription, Curr. Opin. Drug Discov. Dev 12 (2009) 659. [64] F. Vollmuth, W. Blankenfeldt, M. Geyer, Structures of the dual bromodomains of the P-TEFb-activating protein Brd4 at atomic resolution, J. Biol. Chem. 284 (2009) 36547e36556. [65] M. Xu, A. Unzue, J. Dong, D. Spiliotopoulos, C. Nevado, A. Caflisch, Discovery of CREBBP bromodomain inhibitors by high-throughput docking and hit optimization guided by molecular dynamics, J. Med. Chem. 59 (2015) 1340e1349. [66] A. Unzue, M. Xu, J. Dong, L. Wiedmer, D. Spiliotopoulos, A. Caflisch, et al., Fragment-based design of selective nanomolar ligands of the CREBBP bromodomain, J. Med. Chem. 59 (2015) 1350e1356. [67] P. Cieplak, F.-Y. Dupradeau, Y. Duan, J. Wang, Polarization effects in molecular mechanical force fields, J. Phys. Condens. Matter 21 (2009) 333102. [68] J.W. Ponder, C. Wu, P. Ren, V.S. Pande, J.D. Chodera, M.J. Schnieders, et al., Current status of the AMOEBA polarizable force field, J. Phys. Chem. B 114 (2010) 2549e2564. [69] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, et al., A second generation force field for the simulation of proteins, nucleic acids, and organic molecules, J. Am. Chem. Soc. 117 (1995) 5179e5197. [70] G. Lever, D.J. Cole, N.D. Hine, P.D. Haynes, M.C. Payne, Electrostatic considerations affecting the calculated HOMOeLUMO gap in protein molecules, J. Phys. Condens. Matter 25 (2013) 152101. [71] H.N. Banavath, O.P. Sharma, M.S. Kumar, R. Baskaran, Identification of novel tyrosine kinase inhibitors for drug resistant T315I mutant BCR-ABL: a virtual screening and molecular dynamics simulations study, Sci. Rep. (2014) 4. [72] F.A. Pasha, M.M. Neaz, Molecular dynamics and QM/MM-based 3D interaction analyses of cyclin-E inhibitors, J. Mol. Model. 19 (2013) 879e891. [73] J. Blaney, A.M. Davis, Structure-based design for medicinal chemists, in: The Handbook of Medicinal Chemistry, 2014, pp. 96e121. [74] W.A. Cortopassi, K. Kumar, R.S. Paton, Cationep interactions in CREBBP bromodomain inhibition: an electrostatic model for small-molecule binding affinity and selectivity, Org. Biomol. Chem. 14 (2016) 10926e10938.