Computational screening of potential bromodomain-containing protein 2 inhibitors for blocking SARS-CoV-2 infection through pharmacophore modeling, molecular docking and molecular dynamics simulation

2.1 Pharmacophore modeling and database screening

A pharmacophore model was generated using the ZINCPharmer web service (https://zincpharmer.csb.pitt.edu/) (Koes and Camacho, 2012) on the basis of crystal structure of the second bromodomain of human BRD2 bound with an inhibitor of quinazolinone (RVX-208) (Picaud et al., 2013). The PDB file of the crystal structure was downloaded from RCSB PDB database (PDB ID: 4MR5, Resolution: 1.63 Å) (https://www.rcsb.org/structure/4mr6). The crystal structure of RVX-208 in complex with the second bromodomain of BRD2 was used to construct pharmacophore model by ZINCPharmer server. The searching against ZINC15 Purchasable database (Last Updataed 12/20/14) with the generated pharmacophore model was then performed in ZINCPharmer server. The hits were matched the all pharmacophore features and ranked with respect to the fitness scores based on root mean squared deviation (RMSD) between the query features and the hit compound features (Koes et al.,2015). And the compounds screening result was saved in SDF file format for further analysis.

2.2 Molecular docking

Compounds screened from ZINC database with pharmacophore model were docked to BRD2 protein. The screened compounds were added Gasteiger charges and hydrogen atoms and were converted into pdbqt format by Raccoon.py scripts using AutodockTools (Morris et al.,2009) (ADT; Scripps Research Institute, La Jolla,San Diego, USA). The crystal structure of the second bromodomain of human BRD2 bound with a quinazolinone inhibitor (RVX-208) was downloaded from RCSB PDB database (PDB ID: 4MR5, Resolution: 1.63 Å) (https://www.rcsb.org/structure/4mr6). After removing water molecular and co-crystallized ligand, the protein was prepared for docking by adding hydrogen atoms and Gasteiger partial charges using AutodockTools. The prepared protein structure was then converted into pdbqt format. The docking analysis between BRD2 protein and screened compounds was performed by AutoDock Vina (Trott and Olson, 2010) (The Scripps Research Institute, La Jolla, San Diego,USA). A grid box (X = 43.610000, Y = -7.375208, Z = 24.717417, 25*25*25) was generated around the binding sites of protein from co-crystallized ligand (RVX-208) using BIOVIA Discovery Studio Visualizer (Dassault Syst emes, San Diego, CA, USA). The parameters for docking procedures were set to default. The docking pose and docking score were generated from AutoDock Vina.

For each docking analysis, a total of 9 docking conformations are generated, and the conformation with the lowest docking energy is selected as the preferred docking conformation. The lower the binding affinity, indicating that the closer the binding with BRD2, the more stable the docked complex. Interaction between BRD2 and compound was analyzed and visualized using BIOVIA Discovery Studio Visualizer. Before docking procedures, co-crystallized ligand was re-docked to BRD2 protein to validate docking parameters and scheme.

2.3 Molecular dynamics simulation

To evaluate stabilization of the hits docked complexes, top five hits, control compound-BRD2 docked complexes and apo protein (unbound with ligand) were selected to carry out molecular dynamics simulation using Gromacs 2021.4 software (van der Spoel., 2013). The CHARMM36 force field in GROMACS format (charmm36-jul2021 version) (Huang and MacKerell, 2013) was employed for proteins to generate the topology of BRD2, whereas the CHARMM general force field version 4.6 (Vanommeslaeghe et al., 2010) was employed for ligands to generate the topology file using CGenFF server (Vanommeslaeghe and MacKerell, 2012; Vanommeslaeghe et al., 2012) (https://cgenff.umaryland.edu/). Each complex was placed in a dodecahedron box with 1 nm distance to the edges. The explicit water model in the water box is the TIP3P water model (Leontyev and Stuchebrukhov, 2011), and an appropriate amount of chloride ions or sodium ions was added to ensure that the entire simulation system is electrically neutral. The steepest descent energy minimization with 5000 steps was conducted for each system. Next, two steps of equilibration were performed. The first step was a NVT equilibration. In this step, the temperature of the system was gradually heated from 0 K to 310 K using V-rescale thermostat method and the simulation time was 200 ps under the NVT ensemble. The second step was a NPT equilibration. In this step, the pressure of 1.01325105 Pa was controlled using Parrinello-Rahman method and the simulation time was 200 ps under the NPT ensemble. The particle-mesh-Ewald method (Sagui and Darden, 1999) was used to deal with the long-range electrostatic interactions, and van der Waals interactions were evaluated by a cut-off of 10 Å. Finally, the molecular dynamics simulation of 100 ns was carried out under the NPT ensemble. The simulation time step is 2 fs, and the trajectories was saved every 10 ps for subsequent analysis.

2.4 Simulation trajectories analysis

The obtained molecular dynamics simulation trajectories were analyzed to measure root mean square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), hydrogen bond interaction and protein secondary structure using gmx rms, gmx rmsf, gmx gyrate, gmx sasa, gmx hbond and gmx do_dssp module in GROMACS. Principal component analysis and dynamics cross-correlation matrix analysis were performed using bio3d R package (Grant et al.,2021). The free energy landscape analysis was carried out using gmx covar, gmx anaeig and gmx sham module in GROMACS.

2.5 Binding free energy

The binding free energy (ΔG_bind) of each system between ligand and protein was calculated by Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method (Wang et al.,2018) using g_mmpbsa module (Kumari et al.,2014) in GROMACS. The binding free energy is calculated according to the following formula: $$\mathrm{\Delta }{\text{G}}_{\text{bind}}={\text{G}}_{\text{complex}}-\left({\text{G}}_{\text{receptor}}+{\text{G}}_{\text{ligand}}\right)$$

G_complex, G_protein and G_ligand represent the free energy of complex, receptor and ligand, respectively. Each can be expressed separately as following: $$\mathrm{\Delta }{\text{G}}_{\text{bind}}=\mathrm{\Delta }{\text{E}}_{\text{MM}}+\mathrm{\Delta }{\text{G}}_{\text{PB}}+\mathrm{\Delta }{\text{G}}_{\text{SA}}-\text{T}\mathrm{\Delta }\text{S}$$ $$\mathrm{\Delta }{\text{E}}_{\text{MM}}=\mathrm{\Delta }{\text{E}}_{\text{electrostatic}}+\mathrm{\Delta }{\text{E}}_{\text{vdw}}$$

2.6 ADME analysis

The physicochemical properties and ADME (absorption, distribution, metabolism and excretion) properties of the top five hits from docking were analyzed using swissADME web server (Daina et al.,2017) (https://www.swissadme.ch/).

3.1 Construction of the pharmacophore model and pharmacophore-based virtual screening

Pharmacophore feature refers to a certain atom, atomic group, or some active energy characteristics, such as hydrophobic characteristics, in the drug molecules and target structure fusion and plays a key role in drug development (Seidel et al.,2019). A pharmacophore model based on interaction of RVX-208 with the second bromodomain of BRD2 was constructed using ZINCPharmer server and the pharmacophore features were shown in Fig. 1. The pharmacophore hypotheses were chosen for having two hydrogen bond acceptors, one hydrophobic characteristic and one aromatic characteristic. The radius of each hydrogen bond acceptor feature, the hydrophobic feature and aromatic feature were 0.5, 1 and 1.1 Å, respectively. The ZINC database was used to screen hits with the generated pharmacophore model. In total, 10,842 compounds out of 21,723,923 compounds of ZINC database were screened to match the all pharmacophore features using ZINCPharmer server.

Close

Fig. 1

Structure-based pharmacophore model derived from crystal structures of the second bromodomain of human BRD2 bound with a quinazolinone inhibitor (RVX-208) using ZINCPharmer. The hydrogen bond acceptors, the hydrophobic feature and aromatic feature were shown in yellow, green, and purple spheres, respectively. The arrows show the constraint direction.

Export to PPT

3.2 Molecular docking

Before docking procedures, co-crystallized ligand was re-docked to BRD2 protein to validate reliability of docking parameters and scheme. The conformation of re-docked ligand almost overlapped with the original conformation, and root mean square deviation is 0.1 nm (Fig. S1), indicating that the docking parameters and scheme in this study could reliability predict the conformation of protein–ligand docked complex (Hevener et al., 2009). RVX208 (IC₅₀ = 510 nM), co-crystallized ligand of BRD2 protein, was selected as a control compound. In docking analysis, 10,842 compounds screened from pharmacophore model were docked the BRD2 protein and docking affinity was generated.

A total of 102 compounds were screened, which displayed docking affinity to BRD2 higher than the control compound (−8.1 kcal/mol). The top five hits, ZINC20417563, ZINC12322175, ZINC12391479, ZINC13098192, and ZINC20407881, which displayed high affinity against the BRD2 protein with docking score of −10, −9.9, −9.8, −9.7 and −9.7 kcal/mol in comparison of −8.1 kcal/mol for control compound were selected for subsequent analysis. The structure and docking scores of the top five hits were listed in Table 1. We then analyzed the interaction between the five compounds and residues of BRD2 protein using BIOVIA Discovery Studio Visualizer. Interaction details were listed in Table 2. ZINC20417563 molecule had hydrophobic interactions with Phe372, Lys374, Val376, Trp370, Cys425 and Val435 residues, formed two hydrogen bonds with Pro371 and Asn429 residues, and had a Pi-Sulfur interaction attributed to Cys425 residue of BRD2 protein (Fig. 2). In term of ZINC12322175, three hydrogen bonds were yielded with Val376, Asp377 and Asn429 residues, and hydrophobic interactions were formed with Trp370, Pro371, Val376, Leu381, Cys425 and Val435 residues, and a Pi-Sulfur interaction attributed to Trp370 residue of BRD2 protein (Fig. 3). In case of ZINC12391479, the formation of hydrogen bonds interacted with Lys374 and Val376 residues, hydrophobic interactions with Trp370, Pro371, Val376, Leu381, Cys425 and Val435 residues, and two Pi-Sulfur interactions with Trp370 and Cys425 residues were observed (Fig. 4). Fig. 5 showed the interaction diagram of ZINC13098192-BRD2 docked complex. Two residues of BRD2 protein, Trp370 and Pro371, formed two hydrogen bonds interacted with ZINC13098192. Resides Lys374, Val376, Leu381, Leu383, Cys425, His433 and Val435 of BRD2 protein yielded hydrophobic interactions with ZINC13098192. Moreover, a Pi-Sulfur interaction was formed between ZINC13098192 and TRP370 residue of BRD2 protein. As shown in Fig. 6, the BRD2 protein yielded hydrophobic interactions with ZINC20407881 through Lys374, Val376, Phe372, Cys425 and Val435 residues. Moreover, two hydrogen bonds for ZINC20407881 with Pro371 and Asn429 residues and a Pi-Sulfur interaction also formed with Trp370 residue of BRD2 protein were observed. Compared to the RVX208 (exited BRD2 inhibitor) (Picaud et al.,2013) and NSC127133 (a new BRD2 inhibitor identified by Tripathi) (Tripathi et al., 2016), the interactions between the screened compounds in this study and BRD2 include not only hydrogen bonds and hydrophobic interactions but also Pi-Sulfur interactions. Above docking results show that the five compounds can potentially interfere with some of the important residues interactions such as Trp370, Pro371, Leu381 and Cys425 (Picaud et al.,2013). These specific residues are located at the binding interface between substrate and BD2 domain of BRD2. According to docking analysis, the top compounds were occupied in substrate binding pocket of BRD2, and formed hydrogen bond, hydrophobic interaction, and Pi-Sulfur interaction with specific residues located at substrate binding interface. It indicated that the screened compounds can prevent the acetyl-lysine substrate from entering the binding sites to consistent BRD2 activity.

Table 1 The structure and docking score of top five hits.

	Molecular Formula	IUPAC Name	Structure	Affinity (kcal/mol)
ZINC20417563	C27H29N3O3	7-methyl-4-[[4-[[5-methyl-2-(3-methylphenyl)-1,3-oxazol-4-yl]methyl]piperazin-1-yl]methyl]chromen-2-one		−10
ZINC12322175	C27H25N5O2S	N-benzyl-2-(2-benzyl-5-oxo-8,9,10,11-tetrahydro[1]benzothieno[3,2-e][1,2,4]triazolo[1,5-c]pyrimidin-6(5H)-yl)acetamide		−9.9
ZINC12391479	C27H25N5O3S	2-(2-benzyl-5-oxo-8,9,10,11-tetrahydro[1]benzothieno[3,2-e][1,2,4]triazolo[1,5-c]pyrimidin-6(5H)-yl)-N-(3-methoxyphenyl)acetamide		−9.8
ZINC13098192	C27H26FN5OS	[4-[[4-Benzyl-3-(2-fluorophenyl)-5-sulfanylidene-1,2,4-triazol-1-yl]methyl]piperazin-1-yl]-phenylmethanone		−9.7
ZINC20407881	C26H27N3O3	7-methyl-4-({4-[(5-methyl-2-phenyl-1,3-oxazol-4-yl)methyl]piperazin-1-yl}methyl)–2H-chromen-2-one		−9.7

Close

Fig. 2

The 3D interaction (a), docked pose view (b) and 2D interaction (c) profile of ZINC20417563 docked with BRD2 protein. Hydrogen bonds are displayed as green ball and stick, Pi-Sulfur are displayed as gold ball and stick, hydrophobic (pi-pi/pi-alkyl stacking) are displayed as pink ball and stick, and carbon-hydrogen bond are displayed as white ball and stick.

Export to PPT

Close

Fig. 3

The 3D interaction (a), docked pose view (b) and 2D interaction (c) profile of ZINC12322175 docked with BRD2 protein. Hydrogen bonds are displayed as green ball and stick, Pi-Sulfur are displayed as gold ball and stick, hydrophobic (pi-pi/pi-alkyl stacking) are displayed as pink ball and stick, and carbon-hydrogen bond are displayed as white ball and stick.

Export to PPT

Close

Fig. 4

The 3D interaction (a), docked pose view (b) and 2D interaction (c) profile of ZINC12391479 docked with BRD2 protein. Hydrogen bonds are displayed as green ball and stick, Pi-Sulfur are displayed as gold ball and stick, hydrophobic (pi-pi/pi-alkyl stacking) are displayed as pink ball and stick, and carbon-hydrogen bond are displayed as white ball and stick.

Export to PPT

Close

Fig. 5

The 3D interaction (a), docked pose view (b) and 2D interaction (c) profile of ZINC13098192 docked with BRD2 protein. Hydrogen bonds are displayed as green ball and stick, Pi-Sulfur are displayed as gold ball and stick, hydrophobic (pi-pi/pi-alkyl stacking) are displayed as pink ball and stick, and carbon-hydrogen bond are displayed as white ball and stick.

Export to PPT

Close

Fig. 6

The 3D interaction (a), docked pose view (b) and 2D interaction (c) profile of ZINC20407881 docked with BRD2 protein. Hydrogen bonds are displayed as green ball and stick, Pi-Sulfur are displayed as gold ball and stick, hydrophobic (pi-pi/pi-alkyl stacking) are displayed as pink ball and stick, and carbon-hydrogen bond are displayed as white ball and stick.

Export to PPT

3.3 Molecular dynamics simulation

3.3.1

3.3.1 Root mean standard deviation (RMSD)

The root mean standard deviation (RMSD) is a quantitative parameter that can evaluate the stability of protein–ligand system (Stark et al., 2003; Brüschweiler et al., 2003). In order to study the conformational stability of the protein–ligand complex system, the conformation of each frame in the trajectory is compared with the initial conformation, and RMSD of the Cα atom is calculated. As shown in Fig. 7A, the fluctuations in the curves of all complex systems tend to be balanced after 10 ns and remain at ∼ 0.25 nm, indicating that the structures of all systems are in a stable state after a dynamic simulation. The averaged RMSD values of the ZINC13098192, ZINC20417563, ZINC12322175, ZINC12391479, ZINC20407881, control-BRD2 docked complexes and apo BRD2 protein were 0.252 nm, 0.251 nm, 0.247 nm, 0.253 nm, 0.248 nm, 0.253 nm, and 0.246 nm, respectively. The RMSD value of apo BRD2 protein was similar as BRD2 bound with compounds over the simulation, indicating that binding of the compounds had no effect on protein stability.

Close

Fig. 7

Conformational stability analysis of the top five hits, control compound-BRD2 docked complexes and apo protein. (a) the RMSD of the Cα atom, (b) the RMSF profile, (c) the Rg profile, and (d) the SASA profile of the protein of top five hits, control compound-BRD2 docked complexes and apo protein.

Export to PPT

3.4 Hydrogen bond analysis

Hydrogen bonding between proteins and their ligands provides the basic directional and specific interactions for molecular recognition (Korth, 2011). The number of hydrogen bonds between top five hits and control compound and BRD2 protein over simulation time were analyzed using gmx_hbond module in GROMACS. The number of the hydrogen bond formed during protein–ligand complex interaction throughout the 200 ns MD simulation was presented in Fig. 8. ZINC20417563 compound formed up to three hydrogen bonds, ZINC12322175 compound formed up to four hydrogen bonds, ZINC12391479 compound formed up to four hydrogen bonds, ZINC12391479 compound formed up to four hydrogen bonds, ZINC20407881 compound formed up to three hydrogen bonds and control compound formed up to four three bonds, whereas ZINC13098192 compound rarely formed two hydrogen bonds over the 100 ns MD simulation (Fig. 8). These results revealed that the hydrogen bond interaction between all compounds and BRD2 was relatively stable and continuous over simulation time.

Close

Fig. 8

The hydrogen bond analysis between the top five hits, control compound and SGLT2 protein during MD simulation. The hydrogen bond number varies with simulation time.

Export to PPT

Next, hydrogen bond existence maps for six complexes were plotted to explore the consistency of key residues involved in hydrogen bonding (Fig. 9). Fig. 9A showed TRP370, LYS374, SER379, TYR386, ASN429, and HIS433 residues participated in formation of hydrogen bonds interacted with ZINC20417563. Among them, the hydrogen bond between residue ASN433 and atom O3 of ZINC20417563 was seen constantly during the 200 nm simulation. PRO371 and ASN429 residues of BRD2 protein were found to be participated in forming hydrogen bonds with ZINC12322175. And these were observed incessantly throughout the simulation (Fig. 9B). For ZINC12391479, it has been observed to form hydrogen bonds with these residues LYS363, ALA366, TRP370, PRO371, LYS374 and ASP377 of BRD2 protein. Among them, TRP370, LYS374 and ASP377 residues yielded hydrogen bonds with different atoms of ZINC12391479 in which only one hydrogen bond between ASP377 residue and atom N1 of ligand was seen insistently during 0 ∼ 25 ns and 40 ∼ 80 ns simulation time (Fig. 9C). In case of ZINC13098192, BRD2 protein residues TRP370, TYR373, LYS374, ASP377 and HIS433 were found to yield several hydrogen bonds. However, these hydrogen bonds were fleeting and less consistent over simulation time (Fig. 9D). ZINC20407881 formed eight hydrogen bonds with LYS374, TYR386, ASN429 and HIS433 residues in which only one hydrogen bond established between ASN429 residue and atom O3 of ligand was observed constantly throughout the simulation time (Fig. 9E). Similarly, only one hydrogen bond of ASN429 residue of BRD2 protein was found for atom O5 of control compound among several hydrogen bonds formed between BRD2 protein and ligand which was seen consistently over the simulation. (Fig. 9F). These results revealed that the top five hits interacted with the BRD2 with higher stability than control compound. Combining the activity of each system with the hydrogen bond, we believe that the hydrogen bond had a positive effect on the inhibitory activity of the selected hits and BRD2 protein.

Close

Fig. 9

The hydrogen bond existence map for ligands (a) ZINC20417563 (b) ZINC12322175 (c) ZINC12391479 (d) ZINC13098192 (e) ZINC20407881 (f) control compound over simulation.

Export to PPT

3.5 Secondary structure analysis

Protein secondary structure analysis can be utilized to monitor the folding process of protein, the change of protein structure caused by mutation, and the effect of ligand binding on protein structure. The secondary structural elements (a-helix, b-sheet, b-brige, trun, coil, bend, and 5-helix) changes of BRD2 protein bound with hits and control compound during MD simulation were analyzed by dictionary of secondary structure for proteins (DSSP) module (Taghvaei et al.,2022) in GROMACS. DSSP uses hydrogen bonds and special geometric features to identify different secondary structures. The output of DSSP can be displayed as a color map to visually observe the secondary structure changes in each region of the protein during the simulation. The percentage of various secondary structures in the whole simulation time can also be obtained by processing the output data of DSSP. The distribution of the secondary structural elements was seen in Fig. 10, the percentage of various secondary structures was detailed in Table 3, and the secondary structure changes of each region of BRD2 protein were observed in Fig. 11. The secondary structure profile of BRD2 bound with hits and apo BRD2 had similar distribution and changes of the secondary structural elements in most regions. The conserved α-helices and coils with a little of turn and bend elements were mainly formed over majority simulation time in secondary structure profile of BRD2-hit complexes and apo protein. When binding to hits, the percentage of the structure (α-helix + β-sheet + β-bridge + turn) and α-helices were decreased as shown in Fig. 10 and Table 3. The binding of ligands contributed to transforming α-helix to turn and bend in residues 378 ∼ 382 region during the period of 25 ∼ 100 ns simulation as seen in Fig. 11. The binding of the identified compounds caused fluctuation of some residues in BRD2 during simulation between different secondary structural elements. These transformed residues region was located at acetyl-lysine binding pocket at second bromodomain of BRD2. These changes in the secondary structure will affect the conformation of the binding pocket, resulting in the inability of BRD2 to accurately recognize and bind the substrate (acetyl-lysine), thereby reducing the activity of BRD2.

Close

Fig. 10

The percentage of the secondary structure of the protein during the simulation time for (a) ZINC20417563 (b) ZINC12322175 (c) ZINC12391479 (d) ZINC13098192 (e) ZINC20407881 (f) control compound and (g) apo protein.

Export to PPT

Close

Fig. 11

Color maps to for the secondary structure changes of each region of the protein during the simulation time for (a) ZINC20417563 (b) ZINC12322175 (c) ZINC12391479 (d) ZINC13098192 (e) ZINC20407881 (f) control compound and (g) apo protein.

Export to PPT

3.6 Principal component analysis

Studies have shown that although the degree of freedom of the protein is very large, most of its motion is concentrated on a very small number of degrees of freedom. These degrees of freedom form a basic subspace, in which the motion is closely related to the function of the protein. Principal component analysis (PCA) can separate this subspace of proteins, making it possible to understand the functional movement patterns of proteins more clearly (Amadei et al.,1993). Based on the MD simulation trajectory of the protein–ligand complex and apo protein system, the covariance matrix is constructed for Cα atom, and the covariance matrix is then diagonalized into eigenvalues and corresponding eigenvectors. Two eigenvectors with the largest eigenvalue are selected as the projection basis vectors, that is, the direction of functional motion. The simulated trajectory is projected onto the basis vectors, and the two projection coordinates with the largest deviation are stacked to visually see the trend and amplitude of functional motion. The PCA results were visualized by the percentage of variance of eigenvectors graph and essential subspace projection of PC1 vs PC2 graph as shown in Fig. 12. We found that the first three PCs is primarily responsible for the variance seen in the collective motion of BRD2 protein bound with or without ligand. It has been shown that the first three PCs of ZINC13098192, ZINC20417563, ZINC12322175, ZINC12391479, ZINC20407881, control-BRD2 docked complexes account for 51.4 %, 71.6 %, 54.1 %, 51.7 %, 55.8 %, and 56.8 % of the variance, whereas the apo BRD2 protein account for 75.8 % of the variance. It should be observed that the first PC was adequate to account for around 50 % of the overall variance in apo BRD2 protein, whereas contribution of PC2 and PC3 was necessary to account for the identical variance in BRD2-ligand complexes. When analyzing the dominating motions in PCA, the first PC is a vital factor. As seen in the graph (Fig. 12), the first PC of apo BRD2 protein accounts for 53.1 % of the variance, while the ZINC13098192, ZINC20417563, ZINC12322175, ZINC12391479, ZINC20407881, and control-BRD2 complexes account for 23.6 %, 52.2 %, 33.7 %, 26.9 %, 26 %, and 27.4 % of the overall variance, respectively. The apo BRD2 protein had a considerable drop in variance regarding the initial three PCs, which indicated that the unbound protein exhibited significant conformational motions. However, when binding ligand, BRD2 exhibited a reduction in such conformational motions.

Close

Fig. 12

PCA analysis for (a) ZINC20417563-BRD2 complex (b) ZINC12322175-BRD2 complex (c) ZINC12391479-BRD2 complex (d) ZINC13098192-BRD2 complex (e) ZINC20407881-BRD2 complex (f) control compound-BRD2 complex and (g) apo protein. The percentage of variance of eigenvectors graphs (on the right), and essential subspace projection of PC1 vs PC2 graphs (on the left).

Export to PPT

The influence of ligand molecules on the essential dynamics of BRD2 protein was seen in the graphs of essential subspace projection of PC1 vs PC2 (Fig. 12). Protein-ligand complexes had bigger covered area in comparison with the apo BRD2 protein, indicating that a greater portion of the conformational space is sampled, particularly along the first PC. As a result, proteins bound with ligand exhibited more dynamic tertiary structural conformations. Such dynamic tertiary structural conformations are considered to contribute to the formation of stable complexes by ligands in the active pocket of the protein. Furthermore, the region covered by the BRD2 protein bound with hits was also bigger than the protein bound with control compound, showed that this scenario is more flexible than the control case.

3.7 Free energy landscape generated

To depict the energy minimum landscape of control and hit compounds docked BRD2 protein, free energy landscape (FEL) analysis was produced. The Gibbs free energy map of the resulting PC1 and PC2 of eigenvectors was used to investigate the conformational stabilities of the docked complexes as shown in Fig. 13. The values of the Gibbs free energy for ZINC13098192, ZINC20417563, ZINC12322175, ZINC12391479, ZINC20407881, control-BRD2 docked complexes and apo protein ranged from 0 kJ/mol to 12.29 kJ/mol, 14.91 kJ/mol, 12.83 kJ/mol, 11.89 kJ/mol, 11.89 kJ/mol, 15.19 kJ/mol, and 15.44 kJ/mol, respectively. The findings of the global free energy minima demonstrated that the docked complexes displayed more stable states on the folding behavior of BRD2 than apo protein. According to the FEL contour plot (Fig. 13), the binding of selected hits in the binding pocket was more significant and stable than the binding of control compounds. It found that the docked complexes were persistent with significant energy minima, suggesting that the residues involved the docking poses are critical for BRD2 stability and function.

Close

Fig. 13

Free energy landscape for (a) ZINC20417563-BRD2 complex (b) ZINC12322175-BRD2 complex (c) ZINC12391479-BRD2 complex (d) ZINC13098192-BRD2 complex (e) ZINC20407881-BRD2 complex (f) control compound-BRD2 complex and (g) apo protein.

Export to PPT

3.8 Dynamics cross-correlation matrix analysis

The concerted movement of protein domains plays an important role in the binding of proteins to their ligands and other specific biological functions. The atomic motion synergy of domains which can play a specific function in protein system is higher. At present, MD simulation sampling is generally used, the dynamic characteristics of the trajectory are analyzed by statistical methods, and finally the correlation of the system motion is obtained using bio3d R package (Yu and Dalby, 2020). The dynamic cross-correlation matrix (DCCM) was calculated by using the coordinates of Cα atoms from the 100 ns MD trajectories. The DCCM plots showed the difference in correlated motions between BRD2 protein bound with and without compound (Fig. 14). The reduction of correlated motions were observed for BRD2 protein bound with selected five hits and control compound compared to apo protein. The BRD2-ZINCZINC20417563 docked complex displayed a slight decrease in correlated motions, while maximum decrease in correlated motions was seen for BRD2 protein bound with ZINC13098192 compound. The DCCM analysis revealed that the activity of BRD2 protein decreased when binding with compounds.

Close

Fig. 14

The dynamics cross-correlation matrix (DCCM) plots for (a) ZINC20417563-BRD2 complex (b) ZINC12322175-BRD2 complex (c) ZINC12391479-BRD2 complex (d) ZINC13098192-BRD2 complex (e) ZINC20407881-BRD2 complex (f) control compound-BRD2 complex and (g) apo protein. The positive value represents correlated motions in cyan, while the negative values represent the anti-correlated motions in pink.

Export to PPT

3.9 Binding free energy

In order to study the interaction between BRD2 protein and the hits in more detail, binding free energy of each system between hit and protein can be calculated by the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method using the MM-PBSA method in g-mmpbsa tools. We selected 1000 frames uniformly from the 100 ns trajectory of simulation to calculate its binding free energy and the listed the results in Table 4. The order of binding free energiey of six systems is ZINC20417563 (−12.988 kcal/mol) < ZINC20407881 (−9.176 kcal/mol) < ZINC12322175 (−8.47 kcal/mol) < ZINC12391479 (−7.172 kcal/mol) < ZINC13098192 (−7.02 kcal/mol) < control (−6.257 kcal/mol). It indicated that control interacted with the BRD2 with lower stability than selected hits. From the results of binding free energy calculation, it can be seen that the van der Waals force (ΔE_vdw), electrostatic interaction (ΔE_ele) and non-polar solvation free energy (ΔE_nonpolar) are favorable for the binding of protein receptors and ligands, while the polar solvation free energy (ΔE_polar) plays a major adverse effect on the binding. In short, van der Waals force, electrostatic interaction and nonpolar interaction provided the main driving forces for binding these compounds to BRD2.

To further analyze the amino acid residues that play an important role in the binding process, we decomposed the binding free energy between the protein and ligand into each amino acid residue. Most of the key residues in six systems are similar, suggesting that the binding patterns of the six small compounds are not significantly different (Fig. 15). Among the six complexes, Pro371, Val376, Leu381, Leu383, Ile663 have important contributions.

Close

Fig. 15

Residue binding energy contribution of (a) ZINC20417563-BRD2 complex (b) ZINC12322175-BRD2 complex (c) ZINC12391479-BRD2 complex (d) ZINC13098192-BRD2 complex (e) ZINC20407881-BRD2 complex and (f) control compound-BRD2 complex.

Export to PPT

3.10 The ADME properties

The ADME properties refer to the absorption, distribution, metabolism and excretion of drugs in the human body, involving the pharmacokinetic properties of drugs in the body, and are key factors in measuring whether a compound can be a drug (Daoud et al.,2021). The physicochemical properties and ADME properties of the top five hits from docking were assessed using swissADME web server as shown in Table 5. The results showed that pharmacokinetic and ADMET properties of the five selected compounds were in an acceptable range, suggesting these compounds may be potential drug candidates.