Binding Mechanism of Inhibitors to BRD4 and BRD9 Decoded by Multiple Independent Molecular Dynamics Simulations and Deep Learning

Abstract

Bromodomain 4 and 9 (BRD4 and BRD9) have been regarded as important targets of drug designs in regard to the treatment of multiple diseases. In our current study, molecular dynamics (MD) simulations, deep learning (DL) and binding free energy calculations are integrated to probe the binding modes of three inhibitors (H1B, JQ1 and TVU) to BRD4 and BRD9. The MD trajectory-based DL successfully identify significant functional function domains, such as BC-loop and ZA-loop. The information from the post-processing analysis of MD simulations indicates that inhibitor binding highly influences the structural flexibility and dynamic behavior of BRD4 and BRD9. The results of the MM-GBSA calculations not only suggest that the binding ability of H1B, JQ1 and TVU to BRD9 are stronger than to BRD4, but they also verify that van der Walls interactions are the primary forces responsible for inhibitor binding. The hot spots of BRD4 and BRD9 revealed by residue-based free energy estimation provide target sites of drug design in regard to BRD4 and BRD9. This work is anticipated to provide useful theoretical aids for the development of selective inhibitors over BRD family members.

1. Introduction

Bromodomains (BRDs) are evolutionarily conserved epigenetic reader modules that play functional roles in identifying N-acetylated lysine (KAc) residues on histones and other proteins [1,2]. Histone tails undergo various post-translational modifications that participate in the regulation of chromatin accessibility [3,4]. The recognition motifs of BRDs on KAc residues cooperate with other chromatin factors to control gene transcription regulation [5]. Currently, BRDs are divided into eight families, consisting of at least 56 nuclear or cytoplasmic proteins in humans, based on their diverse structures and functions [1,6,7,8]. Previous studies have shown that BRD proteins, particularly those of the bromodomain and extraterminal (BET) family, are involved in tumorigenesis, cancers and inflammatory diseases. Moreover, there is growing evidence to suggest that BRDs play a role in the development of various diseases by regulating the transcription of several genes involved in cancer growth and inflammation [9,10,11]; therefore, BRD proteins have been an important target in regard to anti-cancer drug design. In addition, different BRD inhibitors have been designed, some of which are even undergoing clinical trials for both oncology and non-oncology indications [12,13,14,15]. Therefore, insights into the binding mechanism of inhibitors to BRDs are essential for drug development targeting BRDs.

In spite of large sequence variations, all BRD proteins share a common topology structure with a left-handed bundle of four alpha helices (αZ, αA, αB and αC) connected by functional loop regions (ZA and BC loops) that are responsible for substrate specificity (Figure 1A), and their binding pocket is depicted in Figure 1B. Among the family members of BRDs, BRD4 is one of the most widely studied proteins [14,16,17,18,19], as it binds to transcriptional start sites of genes expressed during the M/G1 transition, which affects mitotic progression [20]. Based on functional studies, abnormalities in BRD4 can result in disorders of gene expression, which play a key role in the occurrence and development process of tumors [10,21,22,23]. Disturbance of BRD4 activity mediated by inhibitors contributes potential therapeutic benefits in regard to cancer, inflammation, and immunology. However, mechanism information on inhibitor-mediated disturbances in the activity of BRD4 is still insufficient. Thus, it is of high significance to probe the atomic-level binding mechanisms of inhibitors to BRD4 for the development of highly secure small molecule inhibitors targeting BRD4.

As a member of the BRD family, BRD9 has received special attention as a druggable subunit of mSWI/SNF chromatin-remodeling complexes. Mutations from various subunits of these complexes have been found in nearly 20% of all human cancers [4,24,25,26]. Karim et al. investigated inhibitor selectivity in the BRD7/9 subfamily, and their results provided a new framework for the development of BRD7/9 inhibitors with improved selectivity or additional polypharmacologic properties [24]. Theodoulou et al. reported the first selective cellular chemical probe (I-BRD9) for BRD9, and this probe can be used for recognizing genes regulated by BRD9 in Kasumi-1 cells, which are involved in oncology and immune response pathways [4]. To date, small-molecule inhibitors targeting KAc sites have been designed for the mSWI/SNF complexes in cancer [27,28,29,30,31]. Despite these advances, the mechanism information regarding the binding of inhibitors to BRD9 is currently limited. Therefore, further exploration of the binding modes of inhibitors to BRD9 is highly significant in terms of the development of clinically available anti-cancer drugs targeting BRDs.

Recently, there has been a focus on anti-cancer drugs regarding insights into the binding selectivity of inhibitors towards proteins of the BRD family [13,31,32,33,34,35,36]. An inhibitor called LP99, reported by Clark et al., has shown a preference for inhibiting the association of BRD7 and BRD9 with acetylated histones in vitro and in cells [27]. Two inhibitors, BT0 and BS6, proposed by Jiang et al., have a selectivity of more than 100-fold for BRD4(1) over BRD4(2). Moreover, their work has indicated that selectively modulating individual BRDs could be a strategy for the treatment of acute gouty arthritis [37]. Su et al. utilized molecular dynamics (MD) simulations and binding free energy calculations to investigate the binding selectivity of inhibitors to BRD9 and BRD4. Their results clarified that two highly flexible loops, the ZA-loop and BC-loop, significantly contribute to the binding selectivity of inhibitors to BRD9 and BRD4 [38]. Wang et al. conducted multiple short MD simulations to study the binding selectivity of inhibitors to BRD9 and TAF1(2). Their work revealed that enthalpy contributions play a critical role in the selective recognition of inhibitors towards BRD9 and TAF1(2) [39]. Considering the important roles of BRD4 and BRD9 in diseases, it is necessary to explore the molecular mechanisms underlying the binding selectivity of inhibitors to BRD4 and BRD9 in order to develop highly selective inhibitors targeting BRDs.

Among the various computational tools used for investigating inhibitor−protein binding, MD simulations [40,41,42,43,44,45,46], Gaussian accelerated molecular dynamics (GaMD) simulations [47,48,49] and binding free energy predictions [50,51,52,53,54] have been popular methods. Furthermore, conventional MD (cMD) and GaMD simulations were applied to successfully explore inhibitor-induced conformational changes [55,56,57,58,59]. Recently, MD simulations have shown a power in applications in regard to discovering cysteine protease inhibitors and molecular-level information on drug design [60]. To extract essential information about protein conformational dynamics from trajectories, both cMD and GaMD simulations have been combined with machine learning and deep learning (DL). Plante et al. utilized MD simulations and machine learning to successfully identify the molecular determinants involved in the classification of G protein-coupled receptor (GPCR) structures and functions [61,62]. Miao’s group developed a method called GaMD, involving deep learning (DL) and a free energy profiling workflow (GLOW), to gain further mechanistic insights into GPCR function, and their results align well with previous experimental and computational studies [63,64]. We will adopt a cMD trajectory-based DL to enhance our understanding of the binding selectivity of inhibitors towards BRD4 and BRD9.

For this study, three inhibitors, namely JQ1, H1B and TVU, were selected to examine their binding selectivity towards BRD4 and BRD9. The structures of the three inhibitors are illustrated in Figure 1C–E. JQ1 is a novel thieno-triazolo-1,4-diazepine compound that possesses a bulky t-butyl ester functional group at C6. Additionally, binding of (+)-JQ1 significantly enhances the thermal stability of all BET family BRDs [36]. H1B, also known as I-BRD9, exhibits more than 700-fold selectivity towards the BET family and over 200-fold selectivity over the highly homologous BRD7 [4,24]. The IC50 value of the TVU against BRD9 is 10 nM, indicating its stronger inhibitory activity on BRD9 [4]. Therefore, gaining insights into the binding selectivity of these three inhibitors towards BRD4 and BRD9 is crucial for designing highly selective inhibitors targeting the BRD family. To achieve our objective, multiple independent MD (MIMD) simulations, DL, principal component analysis (PCA) [65,66,67,68,69] and binding free energy prediction were integrated in this study to decipher the binding selectivity of inhibitors in regard to BRD4 over BRD9. It is anticipated that this study will provide valuable theoretical guidance for the development of drugs targeting the BRD family.

2. Results and Discussion

2.1. Difference in Contacts of Structural Domains Revealed by Deep Learning

The process of the DL is described in Figure S1. Firstly, the conformations recorded in the MD trajectories are transformed into the images used for the DL. Secondly, the images are randomly divided into a training set and a validation set to perform the DL. Lastly, the trained results are visualized in Figure 2, Figure 3 and Figure S2. The details for the DL are clarified in the later method section. The inhibitor-bound BRD4 and BRD9, respectively, performed classification through DL on images extracted from the SMT using the MDTraj program. For the inhibitor-bound BRD4 (Supporting Information Figure S2A,B), the overall accuracy achieves 88.47% on the validation set after 30 epochs, while the overall loss reaches 0.281. Among 4000 snapshots for the validation of each system, three systems can be basically accurately identified (Figure 2A). The 3836 snapshots of H1B-BRD4, the 3859 snapshots of JQ1-BRD4 and the 2911 snapshots of TVU-BRD4 are accurately recognized. The 29 and 135 snapshots of H1B-BRD4 are inaccurately identified as JQ1-BRD4 and TVU-BRD4, respectively. From 4000 snapshots of the validation set of JQ1-BRD4, the 118 snapshots inaccurately recognized as H1B-BRD4 and 23 are inaccurately identified as TVU-BRD4. The 1047 and 42 snapshots from 4000 TVU-BRD4 are erroneously recognized as H1B-BRD4 and JQ1-BRD4, individually. The pixel-attributed residue contact gradient maps of the most populated BRD4 structures are bounded by H1B, JQ1 and TVU (Figure 2B–D), respectively. On the whole, the characteristic residue contacts of H1B-bound BRD4 are located between helixes αB and ZA-loop, as well as αC and ZA-loop (Figure 2B). The characteristic residue contacts of JQ1-bound BRD4 are situated between αC and αZ as well as αB and the ZA-loop (Figure 2C). Compared to the H1B-bound BRD4, binding of JQ1 leads to the appearance of contacts between αC and the ZA-loop. The characteristic residue contacts of TVU-bound BRD4 occur between αC and αZ, αA and αZ, αC and the ZA-loop, as well as αB and the ZA-loop (Figure 2D). By comparison with the H1B-bound BRD4, binding of TVU induces the appearance of two new contacts of αZ with αC and αA. The previous experimental work indicated that the ZA-loop, BC-loop, αC and αZ play important roles in the binding of inhibitors to BRD4 [17,36], which supports our current findings.

With respect to the inhibitor-bound BRD9, the overall accuracy achieves 87.20% on the validation set after 30 epochs, while the overall loss reaches 0.259 (Figure S2C,D). Among 4000 snapshots for the validation of each system, three systems can be accurately identified (Figure 3A). In detail, the 3255 frames of H1B-BRD9, 3582 frames of JQ1-BRD9 and 3616 frames of TVU-BRD9 are accurately identified (Figure 3A). The 40 and 705 frames of H1B-BRD9 are inaccurately recognized as JQ1-BRD9 and TVU-BRD9, separately. Among the 4000 frames of the validation set of JQ1-BRD9, 241 frames are inaccurately identified as H1B-BRD9 and 177 are inaccurately identified as TVU-BRD9. The 367 and 17 snapshots from 4000 TVU-BRD9 are erroneously recognized as H1B-BRD9 and JQ1-BRD9, respectively. The pixel-attributed residue contact gradient maps of the most populated BRD9 structures are associated by H1B, JQ1 and TVU (Figure 3B–D), respectively. In total, the characteristic residue contacts of H1B-bound BRD9 are situated between helixes αB and ZA-loop, as well as αC and αZ (Figure 3B). The characteristic residue contacts of JQ1-bound BRD9 are located between αZ and the ZA-loop, as well as the ZA-loop and BC-loop (Figure 3C). Compared to the H1B-bound BRD9, binding of JQ1 induces the appearance of contacts between αZ and the ZA-loop as well as the ZA-loop and BC-loop, but it leads to the disappearance of contacts between αB and the ZA-loop as well as αC and αZ (Figure 3A,D). The characteristic residue contacts of TVU-bound BRD9 primarily occur between αC and αZ, the BC-loop and ZA-loop, αC and the ZA-loop and αC and the BC-loop (Figure 3D). By comparison with the H1B-bound BRD9, the presence of TVU induces the appearance of the contacts between the BC-loop and ZA-loop, αC and the ZA-loop and αC and the BC-loop, but binding of TVU also leads to the disappearance of the contacts of αC with the ZA-loop and αB (Figure 3B,D). In a previous experimental study, it was observed that the ZA-loop, BC-loop, αC and αZ play important roles in the binding of inhibitors to BRD9, which is in basic agreement with our results [4,70].

According to the above analyses, the trajectory-based DL can rationally identify significant structure domains of BRD4 and BRD9, which is involved in the binding of inhibitors. Meanwhile, the DL also captures the difference between contacts of structural domains between BRD9 and BRD4. Compared to the inhibitor-bound BRD4, more structural contacts were found in the inhibitor-bound BRD9, which may imply the functional difference between BRD4 and BRD9.

2.2. Free Energy Profiles and Structural Dynamics of BRD4 and BRD9

To understand the dynamics behavior of BRD4 and BRD9, PCA was performed on the inhibitor-bound BRD4 and BRD9 by using the coordinates of the Cα atoms recorded at an MD trajectory. The function of eigenvalues as eigenvector indexes was depicted in Figure S3. The eigenvalues are usually used to describe the structural fluctuation of proteins along the first eigenvectors. The first six eigenvalues account for 69.15, 79.58 and 69.87% of the total movements for the H1B-, JQ1- and TVU-bound BRD4 (Figure S3A), respectively, while the first six eigenvalues account for 77.89, 75.99 and 73.49% of the total movements of the H1B-, JQ1- and TVU-bound BRD9 (Figure S3B), individually. On the whole, the first eigenvalues of the inhibitor-bound BRD9 are greater than that of the inhibitor-bound BRD4, indicating that the fluctuation amplitude of BRD9 along the first eigenvector is stronger than that of BRD4.

To probe the difference in free energy profiles of BRD4 and BRD9, the projections of MD trajectories onto the first two eigenvectors were used as reaction coordinates to build free energy landscapes (FELs), with the results being depicted in Figure 4 and Figure S4. MD simulations identify two, three and three energy basins (EBs) in the H1B-, JQ1- and TVU-bound BRD4 (Figure 4A,C,E), respectively, in which JQ1 produces the most significant effect on the free energy profiles of BRD4. By structural comparison of H1B, JQ1 and TVU, it is found that JQ1 brings a net positive charge, which may produce strong electrostatic interactions with those charged residues in BRD4. To check the structural difference in different EBs, the representative structures falling into the EBs were superimposed together (Figure 4B,D,F). The ZA-loop of BRD4 produces the mostly evident difference among different energy states, in particular the JQ1-bound BRD4 (Figure 4D). Meanwhile, two inhibitors, H1B and TVU, hardly yield big deviations in the binding pocket of BRD4 (Figure 4B,F), and only H1B slides slightly. In contrast to H1B and TVU, JQ1 generates large deviations among three representative structures and has three different binding poses.

On the other hand, these results indicate that strong electrostatic interactions caused by a net positive charge of JQ1 produce an obvious effect on the binding poses of JQ1 in the binding pockets of BRD4 and BRD9. For the BRD9, MD simulations detect three, four and two EBs in the H1B-, JQ1- and TVU-bound BRD9 (Figure S4A,C,E) individually. Similar to BRD4, JQ1 also produces the most significant influences on free energy profiles of BRD9. The representative structures trapped in different EBs were aligned together to examine structural differences among different energy states (Figure S4B,D,F). The ZA-loop and BC-loop of BRD9 generate evident deviations between different EBs; moreover, the structural deviations of the ZA-loop are much higher than the BC-loop, and the ZA-loop shows highly disordered states. Three inhibitors, H1B, JQ1 and TVU, yield obvious deviations in the binding pocket of BRD9, in particular JQ1. By comparison, the structural stability of H1B, JQ1 and TVU in the binding pocket of BRD9 is weaker than in BRD4, which may exert significant impacts on inhibitor binding.

In order to probe the dynamics behavior of BRD4 and BRD9, the first eigenvector was visualized and the results were depicted in Figure S5. It is observed that binding of three inhibitors brings about a different effect on the conformations of BRD4 and BRD9, but the common structural domains of BRD4 and BRD9 affected by inhibitor binding are primarily located at the ZA-loop and BC-loop. Compared to H1B and TVU, binding of JQ1 highly strengthens the structural fluctuations of the ZA-loop from BRD4 and BRD9 along the first eigenvector (Figure S5B,E). On the other hand, binding of three inhibitors enhances the structural fluctuation of the ZA-loop in BRD9 by comparison with that of the BRD4 (Figure S5). In addition, binding of JQ1 to BRD9 also strengthens the structural fluctuation of the helix αA relative to binding of JQ1 to BRD4 (Figure S5B,E). Based on the aforementioned analyses, the electrostatic property of JQ1 yields more obvious impacts on conformations of BRD4 and BRD9, especially BRD4.

Based on the current PCA, the binding of H1B, JQ1 and TVU produces a different effect on the dynamics behavior of BRD4 and BRD9, in particular JQ1, indicating that (1) the fluctuation amplitude of BRD9 bounded by three inhibitors along the first eigenvector is stronger than that of BRD4, (2) the structural stability of three current inhibitors in the binding pocket of BRD9 is weaker than in BRD4 and (C) that the presence of three inhibitors strengthens the structural fluctuation of the ZA-loop in BRD9 relative to BRD4. These results indicate that electrostatic properties should be given special attention in future designs of selective inhibitors in regard to the BRD family.

2.3. Structural Property of BRD4 and BRD9

In order to uncover the structural fluctuation, the root-mean-square deviations (RMSDs) of backbone atoms from BRD4 and BRD9 were calculated through the entire MD simulations by referencing the initially optimized structures (Figure 5). The RMSDs of the inhibitor-bound BRD4 fluctuate from 1.17 to 3.33 Å (Figure 5A), while those of the inhibitor-bound BRD9 are located at a range of 1.26–3.96 Å (Figure 5C), indicating that the inhibitor-bound BRD9 shows bigger fluctuations than the inhibitor-bound BRD4. The RMSD of the H1B-BRD4 is distributed at 1.69 and 2.30 Å, while that of the JQ1-BRD4 is populated at 1.69 and 1.30 Å and that of the TVU-BRD4 is located at 1.72 and 2.42 Å (Figure 5C). The RMSD of H1B-BRD9 is situated at 1.70 and 2.61 Å, while that of the JQ1-BRD9 is distributed at two bigger peak values of 2.09 and 3.07 Å and that of the TVU-BRD4 is populated at 1.60 and 2.31 Å (Figure 5D). Based on these results, the presence of three inhibitors produces a more evident influence on the structural stability of BRD9 than BRD4.

The root-mean-square fluctuations (RMSFs) of BRD4 and BRD9 were estimated by using the coordinates of the Cα atoms (Figure 6A,B). It is found that the ZA-loop and BC-loop of BRD4 and BRD9 are highly flexible, especially for the ZA-loop. Compared to H1B and TVU, binding of JQ1 increases the structural flexibility of the ZA-loop from BRD4 and BRD9, but the structural flexibility of the ZA-loop in the JQ1-bound BRD4 is stronger than that in the JQ1-bound BRD9 (Figure 6A,B). Additionally, it is found that the difference in the structural flexibility of the other structural regions is tiny (Figure 6A,B). These results suggest that the electrostatic property of JQ1 highly affect the structural flexibility of the ZA-loop in two BRD proteins.

Molecular surface areas (MSAs) were estimated to examine the impacts of inhibitor associations on the solvent accessible extents of BRD4 and BRD9. The function of MSAs as the simulation time was displayed in Figure S6. The MSAs of the inhibitor-bound BRD4 are located at a range from 5651 to 7378 Ų, while those of the inhibitor-bound BRD9 fluctuate from 5118 to 6862 Ų. The MSAs of the H1B-, JQ1- and TVU-bound BRD4 are distributed at 6266, 6467 and 6265 Ų (Figure 6C), respectively, while the MSAs of the H1B-, JQ1- and TVU-bound BRD9 are situated at the peak values of 5993, 6022 and 5831 Ų (Figure 6C), individually. The results indicate that the contact extent of BRD4 with solvent is higher than BRD9, which shows the difference in hydrophily between BRD4 and BRD9.

The Radius of gyrations (RGs) were computed to check the effect of inhibitor binding on the structural compacts of the BRD4 and BRD9. The time course of RGs for BRD4 and BRD9 was exhibited at Figure S7. The RGs of the inhibitor-bound BRD4 fall into a range of between 13.71 and 14.74 Å, while those of the inhibitor-bound BRD9 fluctuate from 13.25 to 14.46 Å. The RGs of the H1B-, JQ1- and TVU-bound BRD4 are situated at ~14.08 Å (Figure 6D), and the RGs of the H1B-, JQ1- and TVU-bound BRD9 are populated at the peak values of ~13.79 Å (Figure 6D), indicating that the structural compact extent of the inhibitor-bound BRD4 is slightly lower than the inhibitor-bound BRD9. The aforementioned analyses suggested that the electrostatic property of JQ1 exerts certain influences on the structural stability of BRD4 and BRD9.

2.4. MM-GBSA Calculations

The binding affinity of H1B, JQ1 and TVU to BRD4/BRD9 was estimated through the MMPBSA.py in Amber to understand their binding strength, and the results are shown in Table 1 and Table 2. In our current calculation, 400 snapshots extracted from the equilibrated MD trajectories were adopted. It is found that the rank of our calculated results is of good consistence with that determined by the known experimental data. More interestingly, the binding ability of H1B, JQ1 and TVU to BRD9 is stronger than that of those to BRD4.

According to Table 1 and Table 2, the EIs ($∆E_{ele}$) of H1B and TVU with BRD4 are −29.97 and −34.45 kcal/mol, respectively, while that of these two inhibitors with BRD9 are −43.82 and −39.64 kcal/mol, individually, which is favorable for the binding of H1B and TVU to BRD4/BRD9. Additionally, the EIs of H1B and TVU with BRD9 are stronger than that with BRD4. The EI of JQ1 with BRD4 is −22.91 kcal/mol, which is favorable for the JQ1-BRD4 binding, while that of JQ1 with BRD9 is 15.43 kcal/mol, impairing the JQ1-BRD9 binding.

In comparison to the EIs, the polar solvation free energies ($∆G_{gb}$) of the inhibitor-BRD4/BRD9 complexes provide opposite contributions for inhibitor binding. On the whole, the sum ($∆G_{pol}$) of the $∆E_{ele}$ and $∆G_{gb}$ for six current complexes are unfavorable for the binding of inhibitors. It is noted that the VDWIs ($∆E_{vdW}$) of inhibitors with BRD4 and BRD9 are advantageous for the associations of inhibitors. The VDWIs of H1B, JQ1 and TVU with BRD9 are strengthened by 6.96, 11.86 and 5.67 kcal/mol compared to that of H1B, JQ1 and TVU with BRD4, respectively (Table 1 and Table 2). The non-polar solvation free energies ($∆G_{surf}$) also provide advantageous contributions to the associations of inhibitors, and the $∆G_{surf}$ of the H1B-, JQ1- and TVU-BRD9 complexes are individually enhanced by 0.59, 1.08 and 0.36 kcal/mol relative to that of the $∆G_{surf}$ of the H1B-, JQ1- and TVU-BRD4 complexes (Table 1 and Table 2). The entropy contributions ($-T∆S$) of the H1B-, JQ1- and TVU-BRD4 complexes are 18.65, 15.01 and 19.34 kcal/mol, separately, while that of the H1B-, JQ1- and TVU-BRD9 complexes are 25.11, 20.24 and 22.58 kcal/mol, individually (Table 1 and Table 2), which greatly impairs the binding of inhibitors to BRD4/BRD9. By referencing the inhibitor-BRD4 complexes, the $-T∆S$ of the H1B-, JQ1- and TVU-BRD9 complexes are increased by 6.46, 5.23 and 3.24 kcal/mol, respectively, indicating that the entropy contributions produce a greater reduction in the binding ability of inhibitors to BRD9 than BRD4. The sum of two favorable forces ($∆E_{vdW} \mathrm{a}\mathrm{n}\mathrm{d} ∆G_{surf}$) of the H1B-, JQ1- and TVU-BRD9 complexes are strengthened 7.55, 12.94 and 6.03 kcal/mol compared to that of the H1B-, JQ1- and TVU-BRD9 complexes, individually. On the whole, the binding ability of H1B-, JQ1- and TVU to BRD9 is enhanced by 5.29, 5.56 and 6.32 kcal/mol relative to their binding abilities to BRD4. It is concluded that two favorable forces $∆E_{vdW} \mathrm{a}\mathrm{n}\mathrm{d} ∆G_{surf}$ play important roles in the increase in binding ability for H1B-, JQ1- and TVU to BRD9 in comparison to BRD4. Thus, VDWIs and non-polar solvation free energies should be given special attention in the design of highly selective inhibitors in regard to BRD4 over BRD9.

2.5. Interaction Network of Inhibitors with BRD9 and BRD4

To understand the contributions of separate residues to inhibitor-BRD4/BRD9 binding, inhibitor−residue interactions were estimated by means of a residue-based free energy decomposition method, and the results are provided in Figure 7 and Figure S8. The hydrogen bonding interactions (HBIs) of H1B, JQ1 and TVU with BRD4 and BRD9 were examined through the CPPTRAJ program in Amber, and the data is listed in Table 3 and Table 4. The geometric information on hydrophobic interactions and HBIs were displayed in Figure 8 and Figure 9. It was noted that three inhibitors (H1B, JQ1 and TVU) produce binding differences with separate residues in BRD4 and BRD9.

For the H1B-BRD4 complex, H1B produces interactions stronger than 0.8 kcal/mol with six residues, including PRO82, VAL87, LEU92, LEU94, TYR139 LYS141 and ILE146 (Figure 7A,B and Figure S8A). According to Figure 8A, the hydrophobic groups of PRO82, VAL87, LEU92, LEU94 and ILE146 are located near the hydrophobic ring of H1B. Thus, PRO82 structurally forms the π-π interaction of −1.95 kcal/mol with H1B. The alkyls of the residues VAL87, LEU92, LEU94 and ILE146 yield the CH-π interactions and their corresponding interaction energies are −1.51, −1.77, −1.59 and −2.74 kcal/mol, respectively (Figure 7A,B and Figure S8A). Structurally, the phenyl ring of TYR139 is next to the hydrophobic group of H1B, which leads to a π-π interaction of −0.85 kcal/mol. In addition, H1B also forms two HBIs with ASN140, and their occupancy is higher than 87.02%, indicating that these two hydrogen bonds are stable through the entire MD simulations (Table 3 and Figure 9A). A hydrogen bond with an occupancy of 61.73% appears between H1B and LYS141 (Table 3 and Figure 9A), implying that this hydrogen bond is also stable. As for the H1B-BRD9, eight residues generate interactions stronger than 0.8 kcal/mol, involving PHE160, VAL165, ILE169, ALA170, TYR215, ASN216, ARG217 and TYR222 (Figure 7C,D and Figure S7B). The phenyl groups of PHE160, TYR215 and TYR222 are adjacent to the hydrophobic ring of H1B, and, as a result, they are easy to yield to the π-π interactions (Figure 8D). Thus, PHE160, TYR215 and TYR222 contribute the interaction energies of −2.6, −1.07 and −2.76 kcal/mol to the H1B-BRD9 binding (Figure 7C,D and Figure S7B). The alkyls of VAL165, ILE169 and ALA170 are close to the hydrophobic ring of H1B, which leads to the formation of the CH-π interactions between them. Additionally, ILE169 forms a hydrogen bond with H1B, and its occupancy is 54.8% (Table 4 and Figure 9D). Hence, VAL165, ILE169 and ALA170 provide energy contributions of −1.43, −2.13 and −1.19 kcal/mol for the H1B-BRD9 association, respectively (Figure 7C,D and Figure S8B). The positive charge group of ARG217 is situated near the hydrophobic ring of H1B, and they form a cation-π interaction (Figure 8D). Meanwhile, ARG217 also forms a hydrogen bond with an occupancy of 96.82% with H1B (Table 4 and Figure 9D), and therefore ARG217 contributes an interaction energy of −2.12 kcal/mol to the H1B-BRD9 association. Furthermore, ASN216 form two hydrogen bonds with an occupancy of 99.16 and 99.64% with H1B (Table 4 and Figure 9D); as a result, ASN216 yields an interaction of −3.96 kcal/mol with H1B (Figure 7C,D and Figure S8B). Compared to BRD4, H1B produces more interactions with BRD9, which leads to a tighter binding to BRD9 than BRD4.

As for the JQ1-BRD4 complex, the interaction of JQ1 with five residues in BRD4 is stronger than 0.8 kcal/mol, and these five residues include TRP81, PRO82, LEU92, ILE146 and MET149 (Figure 7A,B and Figure S8B). The hydrophobic sidechains of TRP81 and PRO82 are located near the hydrophobic ring of JQ1, and they tend to form an π-π interaction (Figure 8B) that, respectively provides energy contributions of −1.07 and −1.35 kcal/mol for the JQ1-BRD4 binding. It is also observed that the hydrophobic sidechains of LEU92, ILE146 and MET149 are next to the hydrophobic ring of JQ1 (Figure 8B), and thus, structurally, they make it easy to produce the CH-π interactions. As a result, LEU92, ILE146 and MET149 individually contribute the interaction energy of −1.32, 1.4 and −1.17 kcal/mol to the JQ1-BRD4 associations (Figure 7A,B). With respect to the JQ1-BRD9 complex, four residues (PHE160, PHE163, ILE169 and TYR222) are involved in the interactions stronger than 0.8 kcal/mol with JQ1 (Figure 7C,D and Figure S8E). The interaction energies of PHE160, PHE163 and TYR222 with JQ1 are −1.42, −1.01 and −2.43 kcal/mol, which structurally agrees with the π-π interactions between JQ1 and these three residues (Figure 8E). The interaction strength of ILE169 with JQ1 is −1.98 kcal/mol, which structurally matches the CH-π interaction between JQ1 and ILE169 (Figure 8E).

Concerning the TVU-BRD4 complex, nine residues yield interactions stronger than 0.8 kcal/mol with TVU, and these eight residues include PRO82, VAL87, LEU92, LEU94, CYS136, TYR139, ASN140, LYS141 and ILE146 (Figure 7A,B and Figure S8C). According to Figure 8E, the alkyl groups of VAL87, LEU92, LEU94 and ILE146 are close to the hydrophobic ring of TVU, and thus the CH-π interactions are easy to form between them, and they provide energy contributions of −1.86, −1.96, −1.67 and −2.51 kcal/mol for the TVU-BRD4 binding (Figure 7A,B, Figure 8C and Figure S8C). The hydrophobic sidechains of PRO82 and TYR139 are situated near the hydrophobic rings of TVU, which structurally produces the π-π interactions (Figure 8C), and these two π-π interactions contribute the interaction energies of −1.88 and −1.23 kcal/mol to the TVU-BRD4 associations (Figure 7A,B and Figure S8C). The positive charge group of LYS141 structurally forms the cation-π interactions with the hydrophobic ring of TVU (Figure 8C). Meanwhile, LYS141 also forms a hydrogen bond with TVU, and its occupancy is 65.18% (Table 3 and Figure 9B); hence LYS141 provides an energy contribution of −1.17 kcal/mol to the TVU-BRD4 association (Figure 7A,B). In addition, ASN140 yields two HBIs with TVU and their occupancy is higher than 78.73% (Table 3 and Figure 9B), which leads to an interaction energy of −3.18 kcal/mol with TUV (Figure 7A,B). With regard to the TVU-BRD9 complex, TVU produces interactions stronger than 0.8 kcal/mol with eight residues, including PHE160, VAL165, ILE169, ALA170, TYR215, ASN216, ARG217 and TYR222 (Figure 7C,D and Figure S8F). The alkyls of VAL165, ILE169 and ALA170 are adjacent to the hydrophobic ring of TVU, and it is easy to form the CH-π interactions between them (Figure 8F). As a result, VAL165, ILE169 and ALA170, respectively, contribute the interaction energies of −1.98, −2.13 and −1.08 kcal/mol to the TVU-BRD9 associations (Figure 7C,D and Figure S8F). The hydrophobic rings of PHE160, TYR215 and TYR222 are next to that of TVU, which leads to the π-π interactions between them (Figure 8F). PHE160, TYR215 and TYR222 provide energy contributions of −2.39, −1.41 and −3.00 kcal/mol for the binding of TVU to BRD9 (Figure 7C,D and Figure S8F). According to Figure 8F, the positive charged group of ARG217 is located near the hydrophobic ring of TVU, which induces a cation-π interaction between them. In addition, ARG217 forms a hydrogen bond with TVU, and its occupancy is 96.72% (Table 4 and Figure 9D). In total, ARG217 contributes an energy contribution of −1.71 kcal/mol to the TVU-BRD9 association (Figure 7C,D). ASN216 generates two HBIs with TVU, and their occupancy is higher than 69.30%; this residue gives an energy contribution of −3.78 kcal/mol for binding of TVU to BRD9 (Table 4 and Figure 9D).

Based on the aforementioned analyses, it is found that residues (PRO82, PHE160), (LEU92, ILE169) and (ILE146, TYR222) in (BRD4, BRD9) play important roles in the binding of H1B, JQ1 and TVU to BRD4/BRD9. The CH-π, π-π, cation-π interactions and HBIs are responsible for the primary forces for the binding of inhibitors to BRD4 and BRD9; therefore, these common factors involved in inhibitor−residue interactions should be given special attention in future drug designs in regard to the BRD family.

3. Materials and Methods

3.1. Preparation of Simulation Systems

The initial atom coordinates of the JQ1-BRD4, H1B-BRD9 and TVU-BRD9 complexes were taken from the protein data bank (PDB), and they correspond to PDB entry 3MXF [36], 6V1B [24] and 4UIU [4], respectively. Due to differences in the residue sequences in BRD4 and BRD9, residues K59-E162 in BRD4 and residues E138-S239 were adopted to construct our simulation systems. Because of unavailable crystal structure of the JQ1-BRD9 complex in PDB, the JQ1-BRD9 structure was obtained by deleting BRD4 and H1B from the superimposed structures of 3MXF with 6V1B. The H1B-BRD4 and TVU-BRD4 structures were produced by removing BRD9 and JQ1 from the aligned structures of 3MXF with 6V1B and 4UIU. Apart from the inhibitors, the other ligands were cut from the initial model. The protonated states of residues from BRD4 and BRD9 were checked by means of the program H++ 3.0 [71], and rational protonation states were given to each residue of BRD4 and BRD9. The missing hydrogen atoms in the crystal structures were connected to their corresponding heavy atoms by using the Leap module in Amber 22 [72,73]. The ff19SB force field [74] was used to derive the parameters of BRD4 and BRD9. The structures of three inhibitors (JQ1, H1B and TVU) were optimized at a semi-empirical AM1 level, and, subsequently, the BCC charges [75,76] were assigned to each atom of three inhibitors by using the Antechamber module in Amber [77]. The general Amber force field (GAFF2) [78,79] was utilized to get the force field parameters of JQ1, H1B and TVU. The JQ1-, H1B- and TVU-bound BRD4/BRD9 were solved at an octahedral periodic box of water with a buffer of 10.0 Å to embody the solvent environment, while the force field parameters of water molecules were derived from the TIP3P mode [80]. The appropriate number of sodium ions (Na⁺) and chloride ions (Cl⁻) were placed in the water box in 0.15 M NaCl salt concentration to form the neutral simulation systems and the parameters of the Na⁺ and Cl⁻ ions were extracted from the studies of Joung et al. [81,82].

3.2. Multiple Independent Molecular Dynamics

Initialization of six BRD4 or BRD9-related systems possibly results in high energy contacts and orientations between atoms of the entire systems, which may disturb the stability of system simulations. In order to overcome this issue, all of six BRD4- or BRD9-related systems were subjected to a 5000-cycle steepest descent minimization followed by a 10,000-cycle conjugate gradient minimization. The temperature of six optimized systems was softly increased from 0 to 300 K within 2 ns in the canonical ensemble (NVT), and, in this process, all non-hydrogen atoms of the BRD4- or BRD9-related systems were restrained in a weak harmonic restriction of 2 kcal·mol⁻¹·Å². Then, a 3-ns equilibrium process was conducted on six current systems at 300 K under the isothermal–isobaric ensemble (NPT). Then, the 15-ns NPT simulation was implemented to keep the density of the system at 1.01 g/cm³. Subsequently, three 300-ns independent MD simulations were separately performed on six current systems at the NVT with periodic boundary conditions and the particle mesh Ewald method (PME). In each independent MD simulation, the initial atomic velocities were randomly assigned with the Maxwell distribution. To facilitate the post-processing analysis, three independent MD trajectories were joined into a single MD trajectory (SMT). Through all MD simulations, the chemical bonds linking with hydrogen atoms were restricted by adopting the SHAKE algorithm [83]. The Langevin dynamics with a collision frequency of 2.0 ps⁻¹ was executed to control the temperature of six BRD4- or BRD9-related systems [84]. The particle mesh Ewald (PME) method [85] cooperating with an appropriate cutoff value of 12 Å was aided to perform calculations of electrostatic interactions (EIs), and this cutoff was also utilized for treatment of the van der Waals interactions (VDWIs). The program pmemd.cuda [86,87] inlayed in Amber 22 was used to run cMD simulations.

3.3. Deep Learning

To probe impacts of inhibitor binding on internal structures of BRD4 and BRD9, DL was applied to identify the differences in residue contacts. The residue contact map in each MIMD trajectory snapshot of BRD4 and BRD9 was calculated by utilizing a Python packages MDTraj and contact map explorer. A contact definition of ≤4.5 Å between any Cα atoms of two proteins was employed. The derived 104 × 104 (102 × 102 for BRD9) residue contacts were transformed into the 104 × 104 (102 × 102 for BRD9) pixel grayscale images for the analysis by a two-dimensional (2D) convolutional neural network (CNN). In total, 120,000 images were produced for each BRD4-related or BRD9-related system, 80% of which were randomly selected for training, while the rest were used for validation. In this work, the 2D-CNN was built by using the PyTorch package (https://pytorch.org/). This CNN model was mainly composed of 2 convolutional layers of 1 × 1 kernel size, with 16, 32, filters, individually, followed by three fully connected layers, the first two of which included 512 and 128 filters with a dropout rate of 0.5 each. The final fully connected layer was the classification layer for the inhibitor-bound BRD4 or BRD9. Excluding the classification layers, the “ReLu” activation function was utilized in all layers of the 2D-CNN, in which the “softmax” activation function was adopted. A maximum pooling layer with 2 × 2 kernel size was added after each convolutional layer. Lastly, the backpropagation by vanilla gradient-based pixel attribution was employed to estimate the attention map of residue contact gradients so as to discriminate the function difference of BRD4 or BRD9 induced by inhibitor binding, in which the residue contact map was reflected using the most populated structural cluster of each BRD4- or BRD9-related system. We primarily referenced the work of Miao’s group and rewrote the program using PyTorch.

3.4. MM-GBSA Calculations

It is well known that binding free energy is regarded as an indicator to understanding binding strength of ligands to targets. The enthalpy changes (∆H) and entropy changes (−T∆S) are two important factors accompanying the associations of ligands. In this study, binding free energies of H1B, JQ1 and TVU to BRD4/BRD9 were calculated based on the following equation to evaluate binding difference of inhibitors to BRD4 and BRD9

∆G = ∆H − T∆S(1)

By now, the molecular mechanics Poisson−Boltzmann surface area (MM-PBSA) and MM-GBSA are two efficient methods for quickly computing binding-free energies [88,89]. Sun et al. compared the performance of MM-PBSA and MM-GBSA in predictions of binding free energies, and their works uncovered that MM-GBSA obtained a more favorable results than MM-PBSA [90,91]. As for MM-GBSA calculations, ∆H is further divided into two separate components that are included in Equation (2)

(2)$\begin{array}{c}∆\mathrm{H}=∆G_{pol}+∆G_{hydro}\hfill \\ =∆E_{ele}+∆E_{vdW}+∆G_{gb}+∆G_{surf}\hfill \end{array}$

where $∆G_{pol}=∆E_{ele}+∆G_{gb}$ is polar interactions between inhibitors and proteins, while $∆G_{hydro}=∆E_{vdW}+∆G_{surf}$ indicates hydrophobic interactions of inhibitors with proteins. In details, EIs ($∆E_{ele})$ and VDWIs ($∆E_{vdW})$ can be estimated with molecular mechanics and the ff19SB force field. Polar solvation free energies ($∆G_{gb})$ are estimated through the GB model proposed by Onufriev et al. [92]. Non-polar solvation free energies ($∆G_{surf}$) are calculated according to the following empirical Equation (3):

(3)$∆G_{surf}=\gamma \times ∆SASA+\beta$

where the term $∆SASA$ implies the difference in the solvent accessible surface area caused by the associations of inhibitors. The entropy changes $-T∆S$ are calculated by means of the MMPBSA.py program in Amber [93]. The two parameters $\gamma$ and $\beta$ are set as 0.0072 kcal·mol·Å⁻² and 0.0 kcal·mol⁻¹, separately [94].

3.5. Principal Component Analysis

PCA is an important method for deciphering the conformational dynamics of proteins. In this study, PCA was achieved by performing diagonalization on the covariance matrix built using the coordinates of the atoms Cα of BRD4 and BRD9 recoded in the SMT based on Equation (4)

(4)$\mathrm{C}=<(q_i-<q_i>){(q_j-<q_j>)}^T>$

In which $q_i$ and $q_j$ are the Cartesian coordinates of the ith and jth Cα atoms in BRD4 and BRD9, respectively, and the terms $<q_i>$ and $<q_j>$ are their averaged positions over conformational ensembles resulting from MIMD simulations. The eigenvector and the eigenvalue generated by diagonalization, respectively, characterize concerted movement of the structural domains and the fluctuation amplitude along an eigenvector. In this work, the PCA was realized by using the CPPTRAJ program in Amber.

4. Conclusions

BRD4 and BRD9 have been identified critical targets of drug development for treatment of various human diseases. Studying of the atomic-level binding mechanism of inhibitors to BRD4 and BRD9 has proven to be significant. Three independent MD simulations, each for running 400 ns, were carried out to improve conformational samplings of BRD4 and BRD9, and the post-processing analysis based on MD trajectory indicates that the electrostatic property of inhibitors highly affects the structural flexibility of the ZA-loop of BRD4 and BRD9. The MD trajectory-based DL reveals that the ZA-loop and BC-loop play critical roles in the function of BRD4 and BRD9. The calculated binding free energies through the MM-GBSA method not only indicate that VDWIs are responsible primary forces in bindings of inhibitors to BRD4/BRD9, but also suggest that the binding abilities of H1B, JQ1 and TVU to BRD9 are stronger than BRD4. The results from PCA implies that the electrostatic property of JQ1 greatly alters the dynamics behavior of the ZA-loop in BRD4 and BRD9. The interaction network analysis reveals key residues in inhibitor binding and provides targeting sites of drug design in regard to BRD4 and BRD9.

Author Contributions

Conceptualization, J.W., B.W. and J.C.; methodology, J.W., B.W. and J.C.; software, J.W.; validation, J.W., W.Y. and L.Z.; formal analysis, J.W., W.Y. and L.Z.; investigation, J.C.; data curation, J.W. and W.Y.; writing—original draft preparation, J.W.; writing—review and editing, J.W. and J.C.; visualization, W.Y and L.Z.; supervision, W.Y. and L.Z.; project administration, J.W. and J.C.; funding acquisition, J.W. and J.C. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article and Supplementary Materials.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Molecular structures: (A) structural superimposition of inhibitor-bound BRD4/BRD9 complexes, (B) binding pocket of two BRD proteins, (C) H1B, (D) JQ1 and (E) TUV. In this figure, BRD4 and BRD9 are shown in cartoon modes and inhibitors are displayed in stick modes.

Figure 2. Salience maps of key residues detected by DL from MD trajectory of inhibitor-bound BRD4: (A) classification of different BRD4 systems bound by H1B, JQ1 and TVU, (B) salience map of H1B-bound BRD4, (C) salience map of JQ1-bound BRD4 and (D) salience map of TVU-bound BRD4.

Figure 3. Salience maps of key residues detected by DL from the MD trajectory of inhibitor-bound BRD9: (A) classification of different BRD9 systems bound by H1B, JQ1 and TVU, (B) salience map of H1B-bound BRD9, (C) salience map of JQ1-bound BRD9 and (D) salience map of TVU-bound BRD9.

Figure 4. Free energy profiles and representative structures of the inhibitor-bound BRD4: (A,C,E) corresponding to FELs of the H1B-, JQ1- and TVU-bound BRD4, respectively, and (B,D,F) indicating superimposition of representative structures for BRD4 and inhibitors trapped in different EBs. The free energy is scaled in kcal/mol and BRD4 is shown in cartoon modes.

Figure 5. RMSDs of backbone atoms in BRD4 and BRD9: (A) the time course of RMSDs for BRD4, (B) the frequency distribution of RMSDs for BRD4, (C) the time course of RMSDs for BRD9 and (D) the frequency distribution of RMSDs for BRD9.

Figure 6. RMSFs of the Cα atoms from BRD4 and BRD9, molecular surface area and radius of gyration: (A) the RMSFs of BRD4, (B) the RMSFs of BRD9, (C) the frequency distribution of molecular surface area and (D) the frequency distribution for radius of gyrations.

Figure 7. Key residues playing important roles in binding of inhibitors to BRD4 and BRD9: (A) BRD4, (B) the radar representation of inhibitor−residue interactions in BRD4, (C) BRD9 and (D) the radar representation of inhibitor−residue interactions in BRD9.

Figure 8. Hydrophobic interactions of inhibitors with residues in BRD4 and BRD9: (A) H1B-BRD4, (B) JQ1-BRD4, (C) TVU-BRD4, (D) H1B-BRD9, (E) JQ1-BRD9 and (F) TVU-BRD9.

Figure 9. Hydrogen bonding interactions of inhibitors with BRD4 and BRD9: (A) H1B-BRD4, (B) TVU-BRRD4, (C) H1B-BRD9 and (D) TVU-BRD9.

Table 1

Binding free energies of inhibitors to BRD4 obtained by MM-PBSA method.

Complex	H1B		JQ1		TVU
	Average	Std	Average	Std	Average	Std
$∆E_{ele}$	−29.97	9.39	−22.91	15.93	−34.45	9.84
$∆E_{vdw}$	−41.78	3.38	−28.51	7.13	−44.36	2.92
$∆G_{gb}$	42.71	7.61	33.53	15.38	48.59	7.69
$∆G_{surf}$	−5.70	0.44	−3.50	0.90	−5.97	0.33
${}^a∆G_{pol}$	12.74	12.09	10.62	22.14	14.14	12.49
$-T∆S$	18.65	5.52	15.01	6.13	19.34	6.15
${}^b∆G_{bind}$	−16.09	5.81	−6.38	6.54	−16.85	6.44
${}^c∆G_{exp}$			−10.00

Table 2

Binding free energies of inhibitors to BRD9 obtained by MM-PBSA method.

Complex	H1B		JQ1		TVU
	Average	Std	Average	Std	Average	Std
$∆E_{ele}$	−43.82	7.69	15.43	9.03	−39.46	4.91
$∆E_{vdw}$	−48.74	2.89	−40.37	4.25	−50.03	3.28
$∆G_{gb}$	52.35	6.15	−2.66	8.85	50.06	4.22
$∆G_{surf}$	−6.29	0.30	−4.58	0.41	−6.33	0.29
${}^a∆G_{pol}$	3.62	6.80	12.77	12.65	10.60	6.48
$-T∆S$	25.11	4.35	20.24	3.75	22.58	4.15
${}^b∆G_{bind}$	−21.38	6.08	−11.94	5.55	−23.17	5.26
${}^c∆G_{exp}$	−10.16				−11.35

Note: Standard errors are given in parentheses. ${}^a∆G_{pol}=∆E_{ele}+∆G_{gb}$; ${}^b∆G_{bind}=∆E_{ele}+∆G_{pol}+∆E_{vdw}+∆G_{surf}-T∆S$; ${}^c∆G_{exp}$ The experimental values were derived from the experimental Ki values in Reference using the equation ${∆G}_{exp}=-RT\mathrm{l}\mathrm{n}{IC}_{50}$.

Table 3

The hydrogen bonds formed between key residues and inhibitors for BRD4.

Inhibitor	Donor	Acceptor	a Distance(Å)	a Angle(°)	b Occupied(%)
H1B	ASN140:ND2-HD21	H1B:O33	3.10	158.21	87.02
	H1B:N37-H13	ASN:140:OD1	3.96	157.19	90.26
	LYS141:N-H	H1B:O49	2.96	160.86	61.73
TVU	ASN140:ND2-HD21	TVU:O2	2.94	158.70	99.62
	TVU:N1-H4	ASN140:OD1	3.20	152.23	78.73
	LYS141:N-H	TVU:O5	2.94	159.11	65.18

^a The hydrogen bonds are determined by the donor…acceptor atom distance of <3.5 Å and acceptor… H-donor angle of >120°. ^b Occupancy is used to evaluate the stability and strength of the hydrogen bond.

Table 4

The hydrogen bonds formed between key residues and inhibitors for BRD9.

Inhibitor	Donor	Acceptor	a Distance(Å)	a Angle(°)	b Occupied(%)
H1B	ASN216:ND2-HD21	H1B:O33	2.96	155.69	99.16
	H1B:N37-H13	ASN216:OD1	2.93	160.29	99.64
	ARG217:N-H	H1B:O49	2.96	164.08	96.82
	H1B:N35-H12	ILE169:O	3.00	153.78	54.80
TVU	ASN216:ND2-HD21	TVU:O2	2.86	157.25	99.96
	TVU:N1-H10	ASN216:OD1	3.26	152.05	69.30
	ARG217:N-H	TVU:O5	2.96	160.72	96.72

^a The hydrogen bonds are determined by the donor…acceptor atom distance of <3.5 Å and acceptor… H-donor angle of >120°. ^b Occupancy is used to evaluate the stability and strength of the hydrogen bond.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules29081857/s1, Figure S1. Workflow for deep learning on MD trajectories; Figure S2. Learning curves of the training and validation datasets for inhibitor-bound BRD4 and BRD9; Figure S3. The function of eigenvalues over eigenvector indexes; Figure S4. Free energy profiles and representative structures of the inhibitor-bound BRD9; Figure S5. Concerted motions of BRD4 and BRD9 revealed by the first eigenvector from PCA; Figure S6. The time course of molecular surface area; Figure S7. The time evolution for radius of gyrations; Figure S8. Interactions of inhibitors with separate residues in BRD4 and BRD9.

References

1. Filippakopoulos, P.; Picaud, S.; Mangos, M.; Keates, T.; Lambert, J.-P.; Barsyte-Lovejoy, D.; Felletar, I.; Volkmer, R.; Müller, S.; Pawson, T. et al. Histone Recognition and Large-Scale Structural Analysis of the Human Bromodomain Family. Cell; 2012; 149, pp. 214-231. [DOI: https://dx.doi.org/10.1016/j.cell.2012.02.013] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22464331]

2. Zou, Z.; Huang, B.; Wu, X.; Zhang, H.; Qi, J.; Bradner, J.; Nair, S.; Chen, L.F. Brd4 maintains constitutively active NF-κB in cancer cells by binding to acetylated RelA. Oncogene; 2014; 33, pp. 2395-2404. [DOI: https://dx.doi.org/10.1038/onc.2013.179] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23686307]

3. Gardner, K.E.; Allis, C.D.; Strahl, B.D. OPERating ON Chromatin, a Colorful Language where Context Matters. J. Mol. Biol.; 2011; 409, pp. 36-46. [DOI: https://dx.doi.org/10.1016/j.jmb.2011.01.040] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21272588]

4. Theodoulou, N.H.; Bamborough, P.; Bannister, A.J.; Becher, I.; Bit, R.A.; Che, K.H.; Chung, C.-W.; Dittmann, A.; Drewes, G.; Drewry, D.H. et al. Discovery of I-BRD9, a Selective Cell Active Chemical Probe for Bromodomain Containing Protein 9 Inhibition. J. Med. Chem.; 2016; 59, pp. 1425-1439. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.5b00256] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25856009]

5. Chung, C.-W. Small Molecule Bromodomain Inhibitors: Extending the Druggable Genome; Lawton, G.; Witty, D.R. Elsevier: Amsterdam, The Netherlands, 2012; Volume 51, pp. 1-55.

6. Filippakopoulos, P.; Knapp, S. The bromodomain interaction module. FEBS Lett.; 2012; 586, pp. 2692-2704. [DOI: https://dx.doi.org/10.1016/j.febslet.2012.04.045] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22710155]

7. Hewings, D.S.; Rooney, T.P.C.; Jennings, L.E.; Hay, D.A.; Schofield, C.J.; Brennan, P.E.; Knapp, S.; Conway, S.J. Progress in the Development and Application of Small Molecule Inhibitors of Bromodomain–Acetyl-lysine Interactions. J. Med. Chem.; 2012; 55, pp. 9393-9413. [DOI: https://dx.doi.org/10.1021/jm300915b] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22924434]

8. Jahagirdar, R.; Zhang, H.; Azhar, S.; Tobin, J.; Attwell, S.; Yu, R.; Wu, J.; McLure, K.G.; Hansen, H.C.; Wagner, G.S. et al. A novel BET bromodomain inhibitor, RVX-208, shows reduction of atherosclerosis in hyperlipidemic ApoE deficient mice. Atherosclerosis; 2014; 236, pp. 91-100. [DOI: https://dx.doi.org/10.1016/j.atherosclerosis.2014.06.008]

9. Czerwinska, P.; Mackiewicz, A.A. Bromodomain (BrD) Family Members as Regulators of Cancer Stemness & mdash;A Comprehensive Review. Int. J. Mol. Sci.; 2023; 24, 995.

10. Donati, B.; Lorenzini, E.; Ciarrocchi, A. BRD4 and Cancer: Going beyond transcriptional regulation. Mol. Cancer; 2018; 17, 164. [DOI: https://dx.doi.org/10.1186/s12943-018-0915-9]

11. Sabnis, R.W. BRD9 Bifunctional Degraders for Treating Cancer. ACS Med. Chem. Lett.; 2021; 12, pp. 1879-1880.

12. Karim, R.M.; Bikowitz, M.J.; Chan, A.; Zhu, J.-Y.; Grassie, D.; Becker, A.; Berndt, N.; Gunawan, S.; Lawrence, N.J.; Schönbrunn, E. Differential BET Bromodomain Inhibition by Dihydropteridinone and Pyrimidodiazepinone Kinase Inhibitors. J. Med. Chem.; 2021; 64, pp. 15772-15786. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.1c01096] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34710325]

13. Guan, X.; Cheryala, N.; Karim, R.M.; Chan, A.; Berndt, N.; Qi, J.; Georg, G.I.; Schönbrunn, E. Bivalent BET Bromodomain Inhibitors Confer Increased Potency and Selectivity for BRDT via Protein Conformational Plasticity. J. Med. Chem.; 2022; 65, pp. 10441-10458. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.2c00453] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35867655]

14. Davison, G.; Martin, M.P.; Turberville, S.; Dormen, S.; Heath, R.; Heptinstall, A.B.; Lawson, M.; Miller, D.C.; Ng, Y.M.; Sanderson, J.N. et al. Mapping Ligand Interactions of Bromodomains BRD4 and ATAD2 with FragLites and PepLites—Halogenated Probes of Druglike and Peptide-like Molecular Interactions. J. Med. Chem.; 2022; 65, pp. 15416-15432. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.2c01357] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36367089]

15. Xu, Y.; Vakoc, C.R. Targeting Cancer Cells with BET Bromodomain Inhibitors. Cold Spring Harb. Perspect. Med; 2017; 7, [DOI: https://dx.doi.org/10.1101/cshperspect.a026674] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28213432]

16. Duffy, B.C.; Liu, S.; Martin, G.S.; Wang, R.; Hsia, M.M.; Zhao, H.; Guo, C.; Ellis, M.; Quinn, J.F.; Kharenko, O.A. et al. Discovery of a new chemical series of BRD4(1) inhibitors using protein-ligand docking and structure-guided design. Bioorg. Med. Chem. Lett.; 2015; 25, pp. 2818-2823. [DOI: https://dx.doi.org/10.1016/j.bmcl.2015.04.107]

17. Myers, S.M.; Miller, D.C.; Molyneux, L.; Arasta, M.; Bawn, R.H.; Blackburn, T.J.; Cook, S.J.; Edwards, N.; Endicott, J.A.; Golding, B.T. et al. Identification of a novel orally bioavailable ERK5 inhibitor with selectivity over p38α and BRD4. Eur. J. Med. Chem.; 2019; 178, pp. 530-543. [DOI: https://dx.doi.org/10.1016/j.ejmech.2019.05.057]

18. Zhao, L.; Cao, D.; Chen, T.; Wang, Y.; Miao, Z.; Xu, Y.; Chen, W.; Wang, X.; Li, Y.; Du, Z. et al. Fragment-Based Drug Discovery of 2-Thiazolidinones as Inhibitors of the Histone Reader BRD4 Bromodomain. J. Med. Chem.; 2013; 56, pp. 3833-3851. [DOI: https://dx.doi.org/10.1021/jm301793a] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23530754]

19. Su, J.; Liu, X.; Zhang, S.; Yan, F.; Zhang, Q.; Chen, J. A computational insight into binding modes of inhibitors XD29, XD35, and XD28 to bromodomain-containing protein 4 based on molecular dynamics simulations. J. Biomol. Struct. Dyn.; 2018; 36, pp. 1212-1224. [DOI: https://dx.doi.org/10.1080/07391102.2017.1317666]

20. Dey, A.; Nishiyama, A.; Karpova, T.; McNally, J.; Ozato, K. Brd4 Marks Select Genes on Mitotic Chromatin and Directs Postmitotic Transcription. Mol. Biol. Cell; 2009; 20, pp. 4899-4909. [DOI: https://dx.doi.org/10.1091/mbc.e09-05-0380]

21. Yin, M.; Guo, Y.; Hu, R.; Cai, W.L.; Li, Y.; Pei, S.; Sun, H.; Peng, C.; Li, J.; Ye, R. et al. Potent BRD4 inhibitor suppresses cancer cell-macrophage interaction. Nat. Commun.; 2020; 11, 1833. [DOI: https://dx.doi.org/10.1038/s41467-020-15290-0]

22. Andrieu, G.; Tran, A.H.; Strissel, K.J.; Denis, G.V. BRD4 Regulates Breast Cancer Dissemination through Jagged1/Notch1 Signaling. Cancer Res.; 2016; 76, pp. 6555-6567. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-16-0559] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27651315]

23. Torres, F.; Walser, R.; Kaderli, J.; Rossi, E.; Bobby, R.; Packer, M.J.; Sarda, S.; Walker, G.; Hitchin, J.R.; Milbradt, A.G. et al. NMR Molecular Replacement Provides New Insights into Binding Modes to Bromodomains of BRD4 and TRIM24. J. Med. Chem.; 2022; 65, pp. 5565-5574. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.1c01703] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35357834]

24. Karim, R.M.; Chan, A.; Zhu, J.-Y.; Schönbrunn, E. Structural Basis of Inhibitor Selectivity in the BRD7/9 Subfamily of Bromodomains. J. Med. Chem.; 2020; 63, pp. 3227-3237. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.9b01980] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32091206]

25. Middeljans, E.; Wan, X.; Jansen, P.W.; Sharma, V.; Stunnenberg, H.G.; Logie, C. SS18 Together with Animal-Specific Factors Defines Human BAF-Type SWI/SNF Complexes. PLoS ONE; 2012; 7, e33834. [DOI: https://dx.doi.org/10.1371/journal.pone.0033834] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22442726]

26. Kadoch, C.; Hargreaves, D.C.; Hodges, C.; Elias, L.; Ho, L.; Ranish, J.; Crabtree, G.R. Proteomic and bioinformatic analysis of mammalian SWI/SNF complexes identifies extensive roles in human malignancy. Nat. Genet.; 2013; 45, pp. 592-601. [DOI: https://dx.doi.org/10.1038/ng.2628] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23644491]

27. Clark, P.G.K.; Vieira, L.C.C.; Tallant, C.; Fedorov, O.; Singleton, D.C.; Rogers, C.M.; Monteiro, O.P.; Bennett, J.M.; Baronio, R.; Müller, S. et al. LP99: Discovery and Synthesis of the First Selective BRD7/9 Bromodomain Inhibitor. Angew. Chem. Int. Ed.; 2015; 54, pp. 6217-6221. [DOI: https://dx.doi.org/10.1002/anie.201501394] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25864491]

28. Crawford, T.D.; Vartanian, S.; Côté, A.; Bellon, S.; Duplessis, M.; Flynn, E.M.; Hewitt, M.; Huang, H.-R.; Kiefer, J.R.; Murray, J. et al. Inhibition of bromodomain-containing protein 9 for the prevention of epigenetically-defined drug resistance. Bioorg. Med. Chem. Lett.; 2017; 27, pp. 3534-3541. [DOI: https://dx.doi.org/10.1016/j.bmcl.2017.05.063] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28606761]

29. Martin, L.J.; Koegl, M.; Bader, G.; Cockcroft, X.-L.; Fedorov, O.; Fiegen, D.; Gerstberger, T.; Hofmann, M.H.; Hohmann, A.F.; Kessler, D. et al. Structure-Based Design of an in Vivo Active Selective BRD9 Inhibitor. J. Med. Chem.; 2016; 59, pp. 4462-4475. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.5b01865]

30. Picaud, S.; Strocchia, M.; Terracciano, S.; Lauro, G.; Mendez, J.; Daniels, D.L.; Riccio, R.; Bifulco, G.; Bruno, I.; Filippakopoulos, P. 9H-Purine Scaffold Reveals Induced-Fit Pocket Plasticity of the BRD9 Bromodomain. J. Med. Chem.; 2015; 58, pp. 2718-2736. [DOI: https://dx.doi.org/10.1021/jm501893k]

31. Clegg, M.A.; Bamborough, P.; Chung, C.-W.; Craggs, P.D.; Gordon, L.; Grandi, P.; Leveridge, M.; Lindon, M.; Liwicki, G.M.; Michon, A.-M. et al. Application of Atypical Acetyl-lysine Methyl Mimetics in the Development of Selective Inhibitors of the Bromodomain-Containing Protein 7 (BRD7)/Bromodomain-Containing Protein 9 (BRD9) Bromodomains. J. Med. Chem.; 2020; 63, pp. 5816-5840. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.0c00075]

32. Cui, H.; Divakaran, A.; Hoell, Z.J.; Ellingson, M.O.; Scholtz, C.R.; Zahid, H.; Johnson, J.A.; Griffith, E.C.; Gee, C.T.; Lee, A.L. et al. A Structure-based Design Approach for Generating High Affinity BRD4 D1-Selective Chemical Probes. J. Med. Chem.; 2022; 65, pp. 2342-2360. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.1c01779] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35007061]

33. Wu, T.-B.; Xiang, Q.-P.; Wang, C.; Wu, C.; Zhang, C.; Zhang, M.-F.; Liu, Z.-X.; Zhang, Y.; Xiao, L.-J.; Xu, Y. Y06014 is a selective BET inhibitor for the treatment of prostate cancer. Acta Pharmacol. Sin.; 2021; 42, pp. 2120-2131. [DOI: https://dx.doi.org/10.1038/s41401-021-00614-7]

34. Su, J.; Liu, X.; Zhang, S.; Yan, F.; Zhang, Q.; Chen, J. A theoretical insight into selectivity of inhibitors toward two domains of bromodomain-containing protein 4 using molecular dynamics simulations. Chem. Biol. Drug Des.; 2018; 91, pp. 828-840. [DOI: https://dx.doi.org/10.1111/cbdd.13148] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29139214]

35. Wang, L.; Wang, Y.; Yang, Z.; Xu, S.; Li, H. Binding Selectivity of Inhibitors toward Bromodomains BAZ2A and BAZ2B Uncovered by Multiple Short Molecular Dynamics Simulations and MM-GBSA Calculations. ACS Omega; 2021; 6, pp. 12036-12049. [DOI: https://dx.doi.org/10.1021/acsomega.1c00687] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34056358]

36. Filippakopoulos, P.; Qi, J.; Picaud, S.; Shen, Y.; Smith, W.B.; Fedorov, O.; Morse, E.M.; Keates, T.; Hickman, T.T.; Felletar, I. et al. Selective inhibition of BET bromodomains. Nature; 2010; 468, pp. 1067-1073. [DOI: https://dx.doi.org/10.1038/nature09504]

37. Jiang, F.; Hu, Q.; Zhang, Z.; Li, H.; Li, H.; Zhang, D.; Li, H.; Ma, Y.; Xu, J.; Chen, H. et al. Discovery of Benzo[cd]indol-2(1H)-ones and Pyrrolo[4,3,2-de]quinolin-2(1H)-ones as Bromodomain and Extra-Terminal Domain (BET) Inhibitors with Selectivity for the First Bromodomain with Potential High Efficiency against Acute Gouty Arthritis. J. Med. Chem.; 2019; 62, pp. 11080-11107. [DOI: https://dx.doi.org/10.1021/acs.jmedchem.9b01010]

38. Su, J.; Liu, X.; Zhang, S.; Yan, F.; Zhang, Q.; Chen, J. Insight into selective mechanism of class of I-BRD9 inhibitors toward BRD9 based on molecular dynamics simulations. Chem. Biol. Drug Des.; 2019; 93, pp. 163-176. [DOI: https://dx.doi.org/10.1111/cbdd.13398]

39. Wang, L.; Wang, Y.; Yu, Y.; Liu, D.; Zhao, J.; Zhang, L. Deciphering Selectivity Mechanism of BRD9 and TAF1(2) toward Inhibitors Based on Multiple Short Molecular Dynamics Simulations and MM-GBSA Calculations. Molecules; 2023; 28, 2583. [DOI: https://dx.doi.org/10.3390/molecules28062583]

40. Sun, J.; Liu, X.; Zhang, S.; Li, M.; Zhang, Q.; Chen, J. Molecular insights and optimization strategies for the competitive binding of engineered ACE2 proteins: A multiple replica molecular dynamics study. Phys. Chem. Chem. Phys.; 2023; 25, pp. 28479-28496. [DOI: https://dx.doi.org/10.1039/D3CP03392A]

41. Gao, Y.; Zhu, T.; Chen, J. Exploring drug-resistant mechanisms of I84V mutation in HIV-1 protease toward different inhibitors by thermodynamics integration and solvated interaction energy method. Chem. Phys. Lett.; 2018; 706, pp. 400-408. [DOI: https://dx.doi.org/10.1016/j.cplett.2018.06.040]

42. Sun, Z.; Gong, Z.; Xia, F.; He, X. Ion dynamics and selectivity of Nav channels from molecular dynamics simulation. Chem. Phys.; 2021; 548, 111245. [DOI: https://dx.doi.org/10.1016/j.chemphys.2021.111245]

43. Wang, Y.; Yang, F.; Yan, D.; Zeng, Y.; Wei, B.; Chen, J.; He, W. Identification Mechanism of BACE1 on Inhibitors Probed by Using Multiple Separate Molecular Dynamics Simulations and Comparative Calculations of Binding Free Energies. Molecules; 2023; 28, 4773. [DOI: https://dx.doi.org/10.3390/molecules28124773]

44. Xue, W.; Yang, F.; Wang, P.; Zheng, G.; Chen, Y.; Yao, X.; Zhu, F. What Contributes to Serotonin–Norepinephrine Reuptake Inhibitors’ Dual-Targeting Mechanism? The Key Role of Transmembrane Domain 6 in Human Serotonin and Norepinephrine Transporters Revealed by Molecular Dynamics Simulation. ACS Chem. Neurosci.; 2018; 9, pp. 1128-1140. [DOI: https://dx.doi.org/10.1021/acschemneuro.7b00490]

45. Zhang, H.; Ni, D.; Fan, J.; Li, M.; Zhang, J.; Hua, C.; Nussinov, R.; Lu, S. Markov State Models and Molecular Dynamics Simulations Reveal the Conformational Transition of the Intrinsically Disordered Hypervariable Region of K-Ras4B to the Ordered Conformation. J. Chem. Inf. Model.; 2022; 62, pp. 4222-4231. [DOI: https://dx.doi.org/10.1021/acs.jcim.2c00591]

46. Şenol, H.; Ghaffari-Moghaddam, M.; Alim Toraman, G.Ö.; Güller, U. Novel chalcone derivatives of ursolic acid as acetylcholinesterase inhibitors: Synthesis, characterization, biological activity, ADME prediction, molecular docking and molecular dynamics studies. J. Mol. Struct.; 2024; 1295, 136804. [DOI: https://dx.doi.org/10.1016/j.molstruc.2023.136804]

47. Miao, Y.; Feher, V.A.; McCammon, J.A. Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J. Chem. Theory Comput.; 2015; 11, pp. 3584-3595. [DOI: https://dx.doi.org/10.1021/acs.jctc.5b00436] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26300708]

48. Wang, J.; Miao, Y. Peptide Gaussian accelerated molecular dynamics (Pep-GaMD): Enhanced sampling and free energy and kinetics calculations of peptide binding. J. Chem. Phys.; 2020; 153, 154109. [DOI: https://dx.doi.org/10.1063/5.0021399] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33092378]

49. Wang, J.; Arantes, P.R.; Bhattarai, A.; Hsu, R.V.; Pawnikar, S.; Huang, Y.-M.M.; Palermo, G.; Miao, Y. Gaussian accelerated molecular dynamics: Principles and applications. WIREs Comput. Mol. Sci.; 2021; 11, e1521. [DOI: https://dx.doi.org/10.1002/wcms.1521] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34899998]

50. Yang, F.; Wang, Y.; Yan, D.; Liu, Z.; Wei, B.; Chen, J.; He, W. Binding Mechanism of Inhibitors to Heat Shock Protein 90 Investigated by Multiple Independent Molecular Dynamics Simulations and Prediction of Binding Free Energy. Molecules; 2023; 28, 4792. [DOI: https://dx.doi.org/10.3390/molecules28124792]

51. Sun, Z.; Huai, Z.; He, Q.; Liu, Z. A General Picture of Cucurbit[8]uril Host–Guest Binding. J. Chem. Inf. Model.; 2021; 61, pp. 6107-6134. [DOI: https://dx.doi.org/10.1021/acs.jcim.1c01208]

52. Chen, J.; Wang, X.; Pang, L.; Zhang, J.Z.H.; Zhu, T. Effect of mutations on binding of ligands to guanine riboswitch probed by free energy perturbation and molecular dynamics simulations. Nucleic Acids Res.; 2019; 47, pp. 6618-6631. [DOI: https://dx.doi.org/10.1093/nar/gkz499]

53. Xue, W.; Wang, P.; Tu, G.; Yang, F.; Zheng, G.; Li, X.; Li, X.; Chen, Y.; Yao, X.; Zhu, F. Computational identification of the binding mechanism of a triple reuptake inhibitor amitifadine for the treatment of major depressive disorder. Phys. Chem. Chem. Phys.; 2018; 20, pp. 6606-6616. [DOI: https://dx.doi.org/10.1039/C7CP07869B]

54. Jia, Z.-J.; Lan, X.-W.; Lu, K.; Meng, X.; Jing, W.-J.; Jia, S.-R.; Zhao, K.; Dai, Y.-J. Synthesis, molecular docking, and binding Gibbs free energy calculation of β-nitrostyrene derivatives: Potential inhibitors of SARS-CoV-2 3CL protease. J. Mol. Struct.; 2023; 1284, 135409. [DOI: https://dx.doi.org/10.1016/j.molstruc.2023.135409]

55. Miao, Y.; McCammon, J.A. Graded activation and free energy landscapes of a muscarinic G-protein–coupled receptor. Proc. Natl. Acad. Sci. USA; 2016; 113, pp. 12162-12167. [DOI: https://dx.doi.org/10.1073/pnas.1614538113]

56. Wang, J.; Miao, Y. Mechanistic Insights into Specific G Protein Interactions with Adenosine Receptors. J. Phys. Chem. B; 2019; 123, pp. 6462-6473. [DOI: https://dx.doi.org/10.1021/acs.jpcb.9b04867]

57. Chen, J.; Zeng, Q.; Wang, W.; Sun, H.; Hu, G. Decoding the Identification Mechanism of an SAM-III Riboswitch on Ligands through Multiple Independent Gaussian-Accelerated Molecular Dynamics Simulations. J. Chem. Inf. Model.; 2022; 62, pp. 6118-6132. [DOI: https://dx.doi.org/10.1021/acs.jcim.2c00961]

58. Chen, J.; Zhang, S.; Wang, W.; Pang, L.; Zhang, Q.; Liu, X. Mutation-Induced Impacts on the Switch Transformations of the GDP- and GTP-Bound K-Ras: Insights from Multiple Replica Gaussian Accelerated Molecular Dynamics and Free Energy Analysis. J. Chem. Inf. Model.; 2021; 61, pp. 1954-1969. [DOI: https://dx.doi.org/10.1021/acs.jcim.0c01470]

59. Chen, J.; Wang, W.; Sun, H.; He, W. Roles of Accelerated Molecular Dynamics Simulations in Predictions of Binding Kinetic Parameters. Mini-Rev. Med. Chem.; 2024; 24, [DOI: https://dx.doi.org/10.2174/0113895575252165231122095555]

60. dos Santos Nascimento, J.I.; Santana Gomes, J.N.; de Oliveira Viana, J.; de Medeiros e Silva, M.S.Y.; Barbosa, G.E.; de Moura, O.R. The Power of Molecular Dynamics Simulations and Their Applications to Discover Cysteine Protease Inhibitors. Mini-Rev. Med. Chem.; 2024; 24, pp. 1125-1146. [DOI: https://dx.doi.org/10.2174/1389557523666230901152257]

61. Plante, A.; Shore, D.M.; Morra, G.; Khelashvili, G.; Weinstein, H. A Machine Learning Approach for the Discovery of Ligand-Specific Functional Mechanisms of GPCRs. Molecules; 2019; 24, 2097. [DOI: https://dx.doi.org/10.3390/molecules24112097]

62. Plante, A.; Weinstein, H. Ligand-Dependent Conformational Transitions in Molecular Dynamics Trajectories of GPCRs Revealed by a New Machine Learning Rare Event Detection Protocol. Molecules; 2021; 26, 3059. [DOI: https://dx.doi.org/10.3390/molecules26103059]

63. Do, H.N.; Wang, J.; Bhattarai, A.; Miao, Y. GLOW: A Workflow Integrating Gaussian-Accelerated Molecular Dynamics and Deep Learning for Free Energy Profiling. J. Chem. Theory Comput.; 2022; 18, pp. 1423-1436. [DOI: https://dx.doi.org/10.1021/acs.jctc.1c01055]

64. Do, H.N.; Wang, J.; Miao, Y. Deep Learning Dynamic Allostery of G-Protein-Coupled Receptors. JACS Au; 2023; 3, pp. 3165-3180. [DOI: https://dx.doi.org/10.1021/jacsau.3c00503]

65. Ichiye, T.; Karplus, M. Collective motions in proteins: A covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations. Proteins; 1991; 11, pp. 205-217. [DOI: https://dx.doi.org/10.1002/prot.340110305]

66. Amadei, A.; Linssen, A.B.M.; Berendsen, H.J.C. Essential dynamics of proteins. Proteins; 1993; 17, pp. 412-425. [DOI: https://dx.doi.org/10.1002/prot.340170408]

67. Levy, R.M.; Srinivasan, A.R.; Olson, W.K.; McCammon, J.A. Quasi-harmonic method for studying very low frequency modes in proteins. Biopolymers; 1984; 23, pp. 1099-1112. [DOI: https://dx.doi.org/10.1002/bip.360230610]

68. Chen, J.; Wang, L.; Wang, W.; Sun, H.; Pang, L.; Bao, H. Conformational transformation of switch domains in GDP/K-Ras induced by G13 mutants: An investigation through Gaussian accelerated molecular dynamics simulations and principal component analysis. Comput. Biol. Med.; 2021; 135, 104639. [DOI: https://dx.doi.org/10.1016/j.compbiomed.2021.104639]

69. Bao, H.; Wang, W.; Sun, H.; Chen, J. Probing mutation-induced conformational transformation of the GTP/M-RAS complex through Gaussian accelerated molecular dynamics simulations. J. Enzym. Inhib. Med. Chem.; 2023; 38, 2195995. [DOI: https://dx.doi.org/10.1080/14756366.2023.2195995]

70. Weisberg, E.; Chowdhury, B.; Meng, C.; Case, A.E.; Ni, W.; Garg, S.; Sattler, M.; Azab, A.K.; Sun, J.; Muz, B. et al. BRD9 degraders as chemosensitizers in acute leukemia and multiple myeloma. Blood Cancer J.; 2022; 12, 110. [DOI: https://dx.doi.org/10.1038/s41408-022-00704-7]

71. Anandakrishnan, R.; Aguilar, B.; Onufriev, A.V. H++ 3.0: Automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations. Nucleic Acids Res.; 2012; 40, pp. W537-W541. [DOI: https://dx.doi.org/10.1093/nar/gks375]

72. Salomon-Ferrer, R.; Case, D.A.; Walker, R.C. An overview of the Amber biomolecular simulation package. WIREs Comput. Mol. Sci.; 2013; 3, pp. 198-210. [DOI: https://dx.doi.org/10.1002/wcms.1121]

73. Case, D.A.; Cheatham, T.E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K.M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R.J. The Amber biomolecular simulation programs. J. Comput. Chem.; 2005; 26, pp. 1668-1688. [DOI: https://dx.doi.org/10.1002/jcc.20290]

74. Tian, C.; Kasavajhala, K.; Belfon, K.A.A.; Raguette, L.; Huang, H.; Migues, A.N.; Bickel, J.; Wang, Y.; Pincay, J.; Wu, Q. et al. ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution. J. Chem. Theory Comput.; 2020; 16, pp. 528-552. [DOI: https://dx.doi.org/10.1021/acs.jctc.9b00591]

75. Jakalian, A.; Jack, D.B.; Bayly, C.I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem.; 2002; 23, pp. 1623-1641. [DOI: https://dx.doi.org/10.1002/jcc.10128]

76. Jakalian, A.; Bush, B.L.; Jack, D.B.; Bayly, C.I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. J. Comput. Chem.; 2000; 21, pp. 132-146. [DOI: https://dx.doi.org/10.1002/(SICI)1096-987X(20000130)21:2<132::AID-JCC5>3.0.CO;2-P]

77. Wang, J.; Wang, W.; Kollman, P.A.; Case, D.A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model.; 2006; 25, pp. 247-260. [DOI: https://dx.doi.org/10.1016/j.jmgm.2005.12.005]

78. Wang, J.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and testing of a general amber force field. J. Comput. Chem.; 2004; 25, pp. 1157-1174. [DOI: https://dx.doi.org/10.1002/jcc.20035]

79. He, X.; Man, V.H.; Yang, W.; Lee, T.-S.; Wang, J. A fast and high-quality charge model for the next generation general AMBER force field. J. Chem. Phys.; 2020; 153, 114502. [DOI: https://dx.doi.org/10.1063/5.0019056]

80. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys.; 1983; 79, pp. 926-935. [DOI: https://dx.doi.org/10.1063/1.445869]

81. Joung, I.S.; Cheatham, T.E., III. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B; 2008; 112, pp. 9020-9041. [DOI: https://dx.doi.org/10.1021/jp8001614]

82. Joung, I.S.; Cheatham, T.E., III. Molecular Dynamics Simulations of the Dynamic and Energetic Properties of Alkali and Halide Ions Using Water-Model-Specific Ion Parameters. J. Phys. Chem. B; 2009; 113, pp. 13279-13290. [DOI: https://dx.doi.org/10.1021/jp902584c]

83. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H.J.C. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys.; 1977; 23, pp. 327-341. [DOI: https://dx.doi.org/10.1016/0021-9991(77)90098-5]

84. Izaguirre, J.A.; Catarello, D.P.; Wozniak, J.M.; Skeel, R.D. Langevin stabilization of molecular dynamics. J. Chem. Phys.; 2001; 114, pp. 2090-2098. [DOI: https://dx.doi.org/10.1063/1.1332996]

85. Essmann, U.; Perera, L.; Berkowitz, M.L.; Darden, T.; Lee, H.; Pedersen, L.G. A smooth particle mesh Ewald method. J. Chem. Phys.; 1995; 103, pp. 8577-8593. [DOI: https://dx.doi.org/10.1063/1.470117]

86. Salomon-Ferrer, R.; Götz, A.W.; Poole, D.; Le Grand, S.; Walker, R.C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput.; 2013; 9, pp. 3878-3888. [DOI: https://dx.doi.org/10.1021/ct400314y]

87. Götz, A.W.; Williamson, M.J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R.C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput.; 2012; 8, pp. 1542-1555. [DOI: https://dx.doi.org/10.1021/ct200909j]

88. Wang, W.; Kollman, P.A. Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model. J. Mol. Biol.; 2000; 303, pp. 567-582. [DOI: https://dx.doi.org/10.1006/jmbi.2000.4057]

89. Wang, W.; Kollman, P.A. Computational study of protein specificity: The molecular basis of HIV-1 protease drug resistance. Proc. Natl. Acad. Sci. USA; 2001; 98, pp. 14937-14942. [DOI: https://dx.doi.org/10.1073/pnas.251265598]

90. Sun, H.; Li, Y.; Tian, S.; Xu, L.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys.; 2014; 16, pp. 16719-16729. [DOI: https://dx.doi.org/10.1039/C4CP01388C]

91. Sun, H.; Li, Y.; Shen, M.; Tian, S.; Xu, L.; Pan, P.; Guan, Y.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys.; 2014; 16, pp. 22035-22045. [DOI: https://dx.doi.org/10.1039/C4CP03179B]

92. Onufriev, A.; Bashford, D.; Case, D.A. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins; 2004; 55, pp. 383-394. [DOI: https://dx.doi.org/10.1002/prot.20033] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15048829]

93. Miller, B.R., III; McGee, T.D., Jr.; Swails, J.M.; Homeyer, N.; Gohlke, H.; Roitberg, A.E. MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput.; 2012; 8, pp. 3314-3321. [DOI: https://dx.doi.org/10.1021/ct300418h] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26605738]

94. Gohlke, H.; Kiel, C.; Case, D.A. Insights into Protein–Protein Binding by Binding Free Energy Calculation and Free Energy Decomposition for the Ras–Raf and Ras–RalGDS Complexes. J. Mol. Biol.; 2003; 330, pp. 891-913. [DOI: https://dx.doi.org/10.1016/S0022-2836(03)00610-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12850155]