Next Article in Journal
A Multivariate Analysis-Driven Workflow to Tackle Uncertainties in Miniaturized NIR Data
Next Article in Special Issue
Mechanism of Mutation-Induced Effects on the Catalytic Function of TEV Protease: A Molecular Dynamics Study
Previous Article in Journal
Synthesis, Characterization, and Reactivity Studies of New Cyclam-Based Y(III) Complexes
Previous Article in Special Issue
Exploring the Potential of Plant Bioactive Compounds against Male Infertility: An In Silico and In Vivo Study
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Molecules
Volume 28
Issue 24
10.3390/molecules28247997
molecules-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Insight into the Inhibitory Mechanism of Embryonic Ectoderm Development Subunit by Triazolopyrimidine Derivatives as Inhibitors through Molecular Dynamics Simulation

by
Jianan Ju
1,
Hao Zhang
1,
Shanshan Guan
2,3,
Chang Liu
1,
Juan Du
1,
Xiaoli Shen
1 and
Song Wang
1,*
1
Institute of Theoretical Chemistry, College of Chemistry, Jilin University, 2 Liutiao Road, Changchun 130023, China
2
College of Biology and Food Engineering, Jilin Engineering Normal University, Changchun 130052, China
3
Key Laboratory of Molecular Nutrition at Universities of Jilin Province, Changchun 130052, China
*
Author to whom correspondence should be addressed.
Molecules 2023, 28(24), 7997; https://doi.org/10.3390/molecules28247997
Submission received: 16 November 2023 / Revised: 2 December 2023 / Accepted: 5 December 2023 / Published: 7 December 2023
(This article belongs to the Special Issue Molecular Dynamics Simulations of Biomacromolecules)
Browse Figures
Versions Notes

Abstract

:
Inhibition of the Embryonic Ectoderm Development (EED) subunit in Polycomb Repressive Complex 2 (PRC2) can inhibit tumor growth. In this paper, we selected six experimentally designed EED competitive Inhibitors of the triazolopyrimidine derivatives class. We investigated the difference in the binding mode of the natural substrate to the Inhibitors and the effects of differences in the parent nuclei, heads, and tails of the Inhibitors on the inhibitory capacity. The results showed that the binding free energy of this class of Inhibitors was close to or lower compared to the natural substrate, providing an energetic basis for competitive inhibition. For the Inhibitors, the presence of a strong negatively charged group at the 6-position of the parent nucleus or the 8′-position of the head would make the hydrogen atom on the head imino group prone to flip, resulting in the vertical movement of the parent nucleus, which significantly decreased the inhibitory ability. When the 6-position of the parent nucleus was a nonpolar group, the parent nucleus would move horizontally, slightly decreasing the inhibitory ability. When the 8′-position of the head was methylene, it formed an intramolecular hydrophobic interaction with the benzene ring on the tail, resulting in a significant increase in inhibition ability.

Graphical Abstract

1. Introduction

Polycomb Repressive Complex 2 (PRC2) is a conserved multi-protein, repressive chromatin complex essential for developing and maintaining eukaryotic organisms. PRC2 mediates mono-methylation, di-methylation, and tri-methylation of lysine 27 on histone H3 (H3K27me1/2/3), of which H3K27me3 is a repressor mark associated with chromatin compression and site-specific gene silencing [1,2,3,4]. As a critical epigenetic modulator, PRC2 plays a role in multiple biological processes [5], such as cell proliferation, stem cell self-renewal, differentiation, DNA repair, and epithelial-to-mesenchymal transition [6]. However, PRC2 overexpression and activating mutations in PRC2 individual subunits can cause various cancers, including prostate cancer [7], lung cancer [8,9], breast cancer [10], liver cancer [11], gastric cancer [12], ovarian cancer [13], and so on [14]. Therefore, the inhibition of PRC2 activity can inhibit the growth of certain cancers, making it an attractive target for cancer therapy [15,16].
The core subunits of PRC2 are composed of Enhancer of Zeste Homolog 2 (EZH2) or its highly related homolog EZH1, Embryonic Ectoderm Development (EED), Suppressor of Zeste 12 (SUZ12), and Retinoblastoma Binding Protein4/7 (RBBP4/7), shown in Figure 1 [17]. One of them, EZH1/2, has histone methyltransferase activity and is responsible for transferring methyl from the natural methyl donor cofactor S-adenosylmethionine (SAM) to the ε-amino group of the acceptor H3K27. However, it must be catalytically active in the presence of EED and SUZ12. Once EZH1/2 trimethylates H3K27 on a specific nucleosome, the resulting H3K27me3 product binds to EED. This binding then triggers EZH1/2 to denature, making the next catalytic step easier. This process is called the denatured activation of PRC2 [18]. Multiple groups have developed inhibitors targeting EZH2, the primary catalytic subunit in PRC2. Some of these Inhibitors are dual Inhibitors that target both EZH1 and EZH2. This class of Inhibitors competes with SAM and has demonstrated antitumor activity in preclinical and clinical models [19,20,21]. However, new research indicates that when treated with SAM-competitive EZH2 inhibitors, different lymphoma cancer cell lines can develop new mutations in both the wild-type and mutant EZH2. These mutations disrupt the binding of EZH2 inhibitors to PRC2, resulting in PRC2-acquired resistance to these catalytic inhibitors [22,23,24]. These occurrences make current inhibitors that target EZH2 inadequate to inhibit its cancer-causing activity. Based on the catalytic process, the binding of EED to H3K27me3 is crucial for PRC2 [25]. Therefore, targeting the EED subunit could be a promising new approach to inhibiting PRC2 and may offer a solution to the challenges associated with EZH2 inhibitors [26].
EED is a protein containing the WD40 repeat domain, and its C-terminal (81–441 residues) has seven WD40 repeat groups. This motif consists of approximately 40 residues and forms a four-stranded antiparallel β-folded sheet [27]. The repeat seven WD40 domains fold into a typical seven-bladed β-propeller structure. Each surface of the EED has a central pocket with a smaller plane on the top surface and a larger plane on the bottom surface. The co-crystal structure of EED bound to a trimethylated pentapeptide (Ala1-Arg2-Lys3me3-Ser4-Ala5) truncated from H3K27me3 is shown in Figure 2a [25]. This Complex shows that EED has an “aromatic cage” consisting of four aromatic amino acids (Phe97, Tyr148, Trp364, and Tyr365) on its top surface. The ARLme3SA side chain is inserted into this “aromatic cage” and remains to stabilize the combination [28]. Furthermore, the bottom surface of EED can bind to the N-terminal EBD structural domain (residues 39–68) of EZH2 in Figure 2b [29,30].
After the trimethylation of lysine 27 on histone H3 by EZH2, by binding its H3K27me3 product via EED, PRC2 can spread to neighboring nucleosomes to deposit H3K27me3, initiating a positive feedback loop mechanism that efficiently propagates the repressive H3K27me3 histone mark on chromatin. H3K27me3, as a gene silencing marker, can cause various cancers, so inhibiting the interaction between EED and H3K27me3 can effectively inhibit tumor production [5]. Andreas Lingel et al. had proven that compounds that occupied the central pocket on the top surface of EED, which bounds to H3K27me3, could inhibit the EED-H3K27me3 interaction [31]. At present, inhibitors targeting EED mainly include EED226, EEDi-5285, EEDi-5273, and MAK683. EED226 achieved complete tumor regression at a dose of 300 mg/kg in the KARPAS300 xenograft model [32]. In addition, the literature has shown that EED226 could inhibit excessive histone trivalent methylation in liver cancer cells, thereby relieving the inhibition of the mediator of cell death (Bim) and cyclin-dependent kinase inhibitor 1A (p21)’s expression and achieving the effect of inhibiting liver cancer cell proliferation [33]. The Inhibitor EEDi-5285 exhibited excellent pharmacokinetics (PK) characteristics in mice, achieving complete and persistent tumor regression in the KARPAS422 xenograft model at an oral dose of 50 mg/kg [34]. EEDi-5285 was 100 times more potent in binding to EED than the Inhibitor EED226 and 300 times more potent in inhibiting cell growth than the Inhibitor EED226. This suggested that inhibitors with a strong ability to bind to EED had a greater capacity for growth regulation at the cellular level. The Inhibitor EEDi-5273 exhibited good PK and absorption, metabolism, distribution, and excretion (ADME) characteristics, and could achieve complete and persistent tumor regression at 50 mg/kg in the KARPAS422 xenograft model without any signs of toxicity [35]. Huang et al. progressively optimized the Inhibitor EED226 to obtain a more potent and less toxic Inhibitor, MAK683, which was able to achieve complete tumor regression in nude mice at 10 mg/kg in the Karpas 422 xenograft model and was well tolerated at relevant plasma exposure levels [30]. Their optimization process also resulted in the synthesis of a series of triazolopyrimidine derivatives with the same parent nucleus as EED226, and these compounds also possessed competitive inhibitory abilities of the EED protein. However, the mechanism of inhibition of this class of inhibitors at the molecular level was not known. Therefore, in the present work, we chose six of these representative derivatives as subjects for the study of this class of inhibitors. The six Inhibitors have the structure shown in Figure 3 (the atomic numbers of each part are also labeled in the figure), in which the red part is called the head, the blue part is called the parent nucleus, and the green part is called the tail. The different substitutions of the three parts could result in a more than 100-fold differences in inhibition ability between the Inhibitors. Furthermore, Huang et al. utilized a molecular docking method to examine the interaction between the EED protein and Inhibitors 4 and 5 on a molecular scale [30].
Huang et al. also utilized a molecular docking method to examine the interaction between the EED protein and Inhibitors 4 and 5 on a molecular scale [30]. However, the protein was rigid in the docking calculations and could only be studied qualitatively, which did not give binding energies. Additionally, this approach did not clarify why various inhibitors had different abilities to inhibit the EED protein. To address the above problems, six representative Inhibitors of the triazolopyrimidine derivatives, as shown in Figure 3, divided into three comparison groups; among them, Inhibitor 1, Inhibitor 2, and Inhibitor 3 differed only in the substituent of the 6-position on the parent nucleus as the first group; Inhibitor 4 and Inhibitor 5 differed only in the substitution position of the oxygen atom in the head as the second group; and Inhibitor 4 and Inhibitor 6 differed only in the tail as the third group. Using molecular dynamics simulations and binding free energy calculations, we could examine the binding modes and energies of PRC2 with six Inhibitors. Additionally, by analyzing the Inhibitors within each group, we could explore how different substituent groups in the head, parent nucleus, and tail affected their inhibitory ability. The findings of this study can provide valuable guidance for creating new inhibitors.

2. Results

2.1. Molecular Docking Results

Molecular docking results showed that the six Inhibitors were all docked at specific binding sites in aromatic cages of the EED protein and had similar postures (Figure 4a,b). The affinity of Inhibitors 1–6 was −10.2, −10.4, −9.8, −10.7, −10.4, −11.0 kJ·mol−1. All the Inhibitors formed hydrogen bonds to Asn194, Lys211, and Tyr365 in the protein, and Inhibitor 2 also included additional hydrogen bonds to Arg367 (Figure 4c). This result is generally accorded with previous docking results in the literature, differing only in the absence of any mention of the hydrogen bond formed by Arg367 with Inhibitor 2 [30].

2.2. Stability Analysis of the Simulation System

After 200 ns of molecular dynamics simulations for each of the six EED-inhibitor complexes as well as for the EED-pentapeptide complex (EED-ARLme3SA), we have selected one set of parallel simulations from two sets of each system to be discussed (the results obtained from copies were not identical to the master samples, but the main conclusions were basically the same, so we did not discuss them in detail). The Root Mean Square Deviation (RMSD) during the simulation of the seven systems with the Cα atom in the initial conformation as a reference is shown in Figure 5. As can be seen from the figure, the RMSD of each system tends to stabilize after an initial rapid increase and reaches equilibrium after 100 ns, and the fluctuation of the RMSD after equilibrium is less than 0.1 nm.

2.3. Free Energy Landscape and Sampling

The simulation trajectories were projected onto the first two eigenvectors, PC1 and PC2, to create a free energy landscape (FEL), which reflected the changes in conformational energies in the simulated trajectories. The FEL of the six complex systems and EED-ARLme3SA are shown in Figure 6, where the darker blue color represents the lower energy of the conformation. We sampled from the clusters with the lowest free energy and region where the RMSD was as smooth as possible, as the primary goal of the later conformational analysis. Based on the two conditions, seven samples of the simulation systems for the corresponding simulation time were 152.73 ns, 157.30 ns, 116.99 ns, 101.53 ns, 111.16 ns, 136.20 ns, and 62.06 ns. The 3D conformations of the seven samples can be found in Supplementary Materials’ complexes 1~6.pdb and EED-ARLme3SA.pdb. The 3D conformations of all copies are also placed in Supplementary Materials.

2.4. Sample Conformation Analysis

Comparing the sample conformation with the initial conformation obtained by docking, in Complexes 1, 4, and 6, the binding positions of Inhibitors differed little from the initial conformation for Inhibitors. In Complexes 3 and 5, the Inhibitors moved out of the binding pocket vertically, as shown in the enlarged box in Figure 7, making the part inserted into the binding pocket shorter, especially since the vertical movement of Inhibitor 3 was more pronounced than Inhibitor 5. This vertical movement caused the parent nucleus of the Inhibitor to separate from the aromatic cage to some extent and made the tail trimethylamine group more into the solvent. In Complex 2, the Inhibitor underwent horizontal movement perpendicular to the direction of motion in Complexes 3 and 5 (Figure 7). Unlike the vertical movement, this horizontal movement did not separate the parent nucleus of the Inhibitor from the aromatic cage.
A particular phenomenon noticed when making conformational comparisons was that the orientation of the hydrogen atom on the head imino group differed in Complexes 3 and 5 than in the other Complexes. Figure 8a shows the 1-hydrogen atoms on the head imino group of Complexes 3 and 5 are on the opposite side of the C3-N2 bond with the 4-nitrogen atom on the triazole ring, while the other Complexes are on the same side. To ensure this flip was not an isolated phenomenon in the sample, we counted the H1-N2-C3-N4 dihedral angles throughout the simulation. Figure 8b shows that the dihedral angle in the other four Complexes is stable at around 0 degrees. In contrast, Complexes 3 and 5 had a certain percentage of flips, incredibly Complex 3 after 50 ns when almost the entire conformation left the vicinity of 0 degrees, and a part of it had even reached about −200 degrees. This phenomenon was consistent with what we observed in our sampled conformations.
The inversion of the hydrogen atom on the head imino group of Inhibitor 3 might result from the more negative charge of the unique cyano group at the 6-position of its parent nucleus, which attracted the positively charged hydrogen atom to the cyano group side (Figure 8a). Comparison of energy of the Inhibitor conformations before and after the hydrogen atom on the head imino group flip using Gaussian 09 calculations could give a quantitative indication of the strength of this intramolecular interaction. Because the head of Inhibitor 1 was the same as that of Inhibitor 3 and the hydrogen atom did not flip, we removed the head of Inhibitor 3. We attached the head of Inhibitor 1 to Inhibitor 3, naming it Inhibitor 3′, as shown in Figure 9a. Then, Gaussian 09 optimized Inhibitors 3 and 3′ at the B3LYP/6-31G(d) level. To prevent the optimized conformations from changing too much, we fixed the atoms with relatively significant conformational changes to keep our calculated conformations the same as those obtained from the simulations. Calculations showed that the energy of Inhibitor 3 was 40.07 kJ·mol−1 lower than 3′, indicating this intramolecular interaction made this conformation less energetic when the hydrogen flipped to the cyano group side.
Unlike several other Inhibitors, the negatively charged oxygen atom on the head of Inhibitor 5 was closer to the positively charged hydrogen atom on the head imino group, thus attracting the hydrogen atom to turn to the side nearer to it and forming intramolecular interactions (Figure 8a). Inhibitors 4 and 5 had the same number of atoms but different oxygen atom positions. We removed the head of Inhibitor 5, attached the head of Inhibitor 4 to Inhibitor 5, and exchanged the oxygen atom to the same side as Inhibitor 5, naming it Inhibitor 5′, as shown in Figure 9b. Then, Gaussian 09 optimized Inhibitors 5 and 5′ at the B3LYP/6-31G(d) level. Calculations showed that the energy of Inhibitor 5 was 25.22 kJ·mol−1 lower than that of 5′, indicating this intramolecular interaction made the conformation less energetic when the hydrogen atom flipped to the oxygen atom side. The higher energy of the conformational transition of Inhibitor 3 compared to 5 also indicated why there was a more significant flipping ratio of the hydrogen atom on the head imino group of Inhibitor 3 than 5.
The inversion of the hydrogen atom on the head imino group could result in a more considerable distance between the parent nucleus and the head benzene ring, and the head of this class of Inhibitors was relatively fixed in its protein-binding position so that the part of the parent nucleus moved toward the outside of the protein. The previously observed vertical movement of the parent nucleus in Inhibitors 3 and 5 was partly related.

2.5. Hydrogen Bond Analysis

The hydrogen bonds between the protein and Inhibitor are often an essential component of affinity. In Table 1, we counted the occurrence rates of hydrogen bonding for some vital residues throughout the simulations. In this case, the statistics for the occurrence of hydrogen bonding for important residues are summed over all types of hydrogen bonds formed by that residue with the inhibitor, and there may be more than one hydrogen bond present in some snapshots, which is why the occurrence of hydrogen bonding for Asn194 exceeds 100% in Complexes 4 and 6.
As can be shown in Table 1, there is a significant reduction in the hydrogen bonding occurrence of Asn194 in Complexes 3 and 5, which may be related to the inversion of the hydrogen atom on the head imino group in Complexes 3 and 5. The reduction in hydrogen bonding might be one of the reasons why the inhibition of these two Inhibitors was lower than that of other Inhibitors. The occurrence rate of hydrogen bonding for other vital residues will be discussed with the binding free energy in Section 2.6 below.

2.6. Binding Free Energy Analysis

Binding free energy can determine quantitatively the strength of the interaction between Inhibitor and protein. The six Inhibitors studied in this paper are all competitive, and their inhibition abilities are closely related to their affinity with protein. We used the MM/PBSA method [36] to calculate the binding free energy of six Inhibitors and the EED protein. Various binding free energy contributions between the EED protein and Inhibitors of the six complex systems are shown in Table 2, and the order of their total binding free energy contributions ΔGbinding is consistent with the sequence of Inhibitor capabilities measured by experiments. In all complex systems, the Van der Waals Force contribution was much more significant than the electrostatic contribution, indicating that the primary binding of protein to Inhibitors relied on Van der Waals forces.
In addition, Table 2 shows that the Van der Waals Force and electrostatic contributions of Complex 3 and Complex 5 are significantly smaller than those of the other Complexes. The decrease in their Van der Waals Force contribution might be due to a change in the Inhibitor binding site, which detached its parent nucleus from the aromatic cage. At the same time, the decrease in electrostatic interactions might be due to the inversion of the hydrogen atom on the head imino group, which reduced the occurrence of hydrogen bonding with Asn194 or the electrostatic interactions with other charged residues.
To identify critical residues in the EED protein that interacted with Inhibitors, we also calculated the decomposition contribution of each residue to the total binding free energy in each of the six systems. Table 3 shows the EMM values of residues with more significant contributions.
Table 3 illustrates that the three amino acids Tyr365, Tyr148, and Phe97, which formed the aromatic cage, played a significant role in the EMM of each Complex. These amino acids were the main reason for the inhibitory ability of this class of competitive inhibitors. Conformational observations showed that the aromatic rings of Tyr365 and Tyr148 in all Complexes sandwiched the aromatic ring of the Inhibitor parent nucleus to form a stable sandwich structure (Figure 10a–c). There was a large π-π stacking interaction between these aromatic rings. The EMM contributions of these two residues in Complexes 3 and 5 (the two weakest Inhibitors) were significantly smaller than in the other four Complexes because the parent nucleus of Inhibitors partially moved away from the aromatic cage, making the π-π stacking interaction between the parent nucleus and Tyr365 and Tyr148 weaker (Figure 10a,c). The EMM contributions of Tyr365 and Tyr148 in Complexes 4 and 6 (the two most potent Inhibitors) were significantly higher than in the other four Complexes, possibly because the extra methylene group on the dihydrofuran ring at the head of the Inhibitors 4 and 6 had more hydrophobic interactions with these two residues than in the other Complexes. Phe97 formed mainly the π-π stacking interaction with the aromatic ring of the Inhibitor tail or hydrophobic interaction with the parent nucleus. Its EMM contributions were also minimized in Complexes 3 and 5, as these two Complexes moved outward and disrupted the π-π stacking interaction between Phe97 and the Inhibitor tail (Figure 10d). However, Phe97 had the more considerable EMM contribution in Complexes 1 and 6, probably because the benzene rings at the tail in Complexes 1 and 6 were closer to this residue and formed stronger interactions (Figure 10e,f). The last residue constituting the aromatic cage, Trp364, had a weak EMM contribution in all the Complexes, suggesting that although it was closer to the Inhibitor in position, the spacing of Tyr365 in the middle did not allow it to form more robust interactions with the Inhibitors.
Lys211, Glu238, and Asp237 were located on the same side of the Inhibitors. The charged side chains of these three residues were very close, where positively charged Lys211 formed a salt bridge with negatively charged Glu238, and Asp237 was next to Glu238, also negatively charged but slightly farther away from the Inhibitors than Glu238 (Figure 11). The group that primarily interacted with Lys211 was the more negatively charged 2-nitrogen atom on the parent nucleus of the Inhibitors (the 1-nitrogen atom also had a weaker interaction). In Complex 2, due to the horizontal movement of the Inhibitor parent nucleus in the direction of Lys211, its 2-nitrogen atom on the Inhibitor was closest to the positively charged center nitrogen atom of Lys211 (0.43 nm), and its EMM contribution was the largest (Figure 11a). In Complexes 4 and 5, the distance between the 2-nitrogen atom on the Inhibitor and the positively charged center nitrogen atom of Lys211 was not much different (0.50 nm and 0.53 nm), and these two nitrogen atoms could form a certain percentage of hydrogen bonds, whose EMM contribution was close. This distance was relatively far in Complexes 1, 3, and 6 (0.58 nm, 0.55 nm, and 0.56 nm) and could not form hydrogen bonds, so the EMM contributions were similar and smaller.
Glu238 and Asp237, as opposed to Lys211, had negative side charges. The greater the electrostatic attraction exhibited by Lys211, the greater the repulsive force indicated by Glu238 and Asp237, which was more pronounced for Glu238 due to its proximity to the Inhibitor. However, Glu238 in Complexes 3 and 5 exhibited additional repulsive forces, which might be related to the different orientations of the hydrogen atom on the head imino group. The shift of the hydrogen on the imine with electrostatic attraction to Glu238 to the other side led to a decrease in the electrostatic attraction of the imine with Glu238 and an increase in the repulsive force (Figure 11b,c). In addition, in Complex 6, Lys211 and Glu238 exhibited additional electrostatic attraction compared with other Complexes, which might be because the Inhibitor tail was a pyridine ring with more uneven charge distribution, and the charged atoms on the ring could have additional interactions with Lys211 and Glu238 (Figure 11d).
The positively charged guanidine group of Arg414 was located on one side of the pyrimidine ring of the parent nucleus, which was close to the 6-position cyano group on the parent nucleus and the oxygen atom on the head. It could form an electrostatic interaction with the charged groups in these two parts (Figure 11). Among them, the 6-position cyano group in Complex 3 was the closest to this residue and carried an enormous negative charge on the nitrogen atom of the cyano group (0.53e). Therefore, the EMM contribution of Arg414 in Complex 3 was the largest (Figure 11b). Because Arg414 was in the outer layer of the EED protein, its more interactions with the 6-position cyano group caused Inhibitor 3 to move outward overall by its attraction, so the vertical movement was more significant in Complex 3 than in Complex 5. In contrast to Complex 3, in Complex 2 with no polar nitrogen atom at the 6-position (the charge on the C-H group was only 0.11e), Arg414 had a minor interaction with the parent nucleus and a little EMM contribution, resulting in Inhibitor 2 being more prone to horizontal movement due to the lack of fixation of this interaction (Figure 11a). The remaining four complexes were all nitrogen atoms at the 6-position, with the EMM contribution of Arg414 intermediate between Complexes 2 and 3.
Asn194 was also the residue with a more considerable EMM contribution to the binding free energy. Previous hydrogen bonding analysis showed that this residue formed more significant occurrence rates of hydrogen bonding with the hydrogen atom on the head imino group of the Inhibitors (Figure 12a). The hydrogen bonding analysis also showed a considerable decrease in the Asn194 hydrogen bonding occurrence rates in Complexes 3 and 5, related to the flipping of the hydrogen atom on the head imino group (Figure 12b). Consequently, the EMM contribution of Asn194 was smaller in these two Complexes than in the other four Complexes.
Arg367, Asp310, and Met366 were all located near the heads of the inhibitors, and the first two could form a variety of interactions with charged groups in the inhibitor due to their charged side chains. Among them, the electrostatic interactions were strongest with the oxygen atoms of the head that are near them (Arg367 for electrostatic attraction and Asp310 for repulsive forces). In Met366, on the other hand, it was the -NH- in the peptide plane that formed electrostatic attraction with the head oxygen atom. In Complex 2, due to the different positions of the head oxygen atom (as shown in Figure 12), all three residues in this complex exhibit different EMM contributions than in the other complexes.
Asn146 was close to the tails of the Inhibitors. The tail of Inhibitor 6 was not a benzene ring but a pyridine ring (Figure 12d). In the simulation process, the pyridine ring was prone to rotation. When the nitrogen atom with a relatively concentrated negative charge was transferred to the side of Asn146, a certain proportion (17.61%) of a hydrogen bond would be formed with the amino group on Asn146. Therefore, its EMM contribution was more significant. The aromatic rings at the tail of the remaining five Complexes were all benzene rings with relatively dispersed charges. However, they all carried polar trimethylamine groups, which could also interact with Asn146. The tail trimethylamine of Inhibitor 1 was significantly closer to Asn146 than the others (Figure 12c), so its EMM contribution was also more substantial.
Asp430 was located between the parent nucleus 6-position and the tail of the inhibitor. Its negatively charged side chain could exhibit repulsion with the negatively charged group at the 6-position of the parent nucleus. In addition, in Complex 6, this residue also exhibited repulsion with the nitrogen atom on the tailed aryl ring of the inhibitor (as shown in Figure 12f), whereas in the other five complexes, it would exhibit electrostatic attraction with the tailed trimethylamine group. Thus, inhibitor 6, which had a nitrogen atom on the tail aromatic ring of the inhibitor and no trimethylamine group, had the greatest repulsive force with Asp430 (as shown by the EMM contribution having the largest positive value). Inhibitor 3, which had the trimethylamine group closest to Asp430, exhibited the greatest electrostatic attraction with Asp430.

3. Discussion

3.1. Comparison of Inhibitors and Natural Substrate Binding Modes

Comparing the binding modes of inhibitors and natural substrate to EED protein provide more insight into the main sources of inhibitor competitiveness. We performed simulations of EED-ARLme3SA and calculated its binding free energy. The results showed that the binding posture of the ARLme3SA on EED was quite different from that of the Inhibitors of the triazolopyrimidine derivatives class. This might be due to the large number of residues in the EED protein with charged side chains (44 negatively charged glutamate/aspartate and positively charged lysine/arginine out of a total of 362 residues). Of these, the natural substrate binding was more negatively charged on the top surface and more positively charged on the side of the protein and the bottom of the center pocket. However, the natural substrate had a positive charge on both Arg2 and Lys3me3, which repelled the positively charged residues at the bottom of the central pocket, and it did not interact strongly with the aromatic cage due to the relatively small number of nonpolar side chains on the natural substrate. As a result, the Lysme3 side chain in the natural substrate was inserted more shallowly into the central pocket compared to the Inhibitors of the triazolopyrimidine derivatives class (as shown in Figure 13) and instead tended to form ionic interactions with the region of negative charge concentration on the top surface.
The results of binding free energy calculation could prove the above conclusions quantitatively. The Van der Waals interaction between the protein and the natural substrate was only −13.95 kJ·mol−1, which was much smaller than the electrostatic interaction between the two (−162.93 kJ·mol−1), which was the opposite of the results for the six Inhibitors. The decomposition contribution of each residue indicated that the residues that formed the strongest interactions with the natural substrate were those with charged side chains, such as Glu238 (−38.65 kJ·mol−1), Asp310 (−30.72 kJ·mol−1) and Asp237 (−27.74 kJ·mol−1) that exhibited electrostatic attraction, and Lys211 (33.57 kJ·mol−1), Arg367 (29.35 kJ·mol−1) and Arg414 (28.89 kJ·mol−1) that exhibited repulsive forces. The first three of these residues formed a negatively charged pocket on the top surface and were the primary region for binding to the positively charged Arg2 and Lys3me3 in the natural substrate. Arg367 was located at the bottom of the central pocket and was the primary residue that resisted the penetration of Lys3me3 into the central pocket. In contrast, Tyr365, Trp364, and Tyr148 among the aromatic cage residues formed the largest hydrophobic interactions with the natural substrate, but all of them had smaller EMM contributions of −2.8142 kJ·mol−1, −2.2539 kJ·mol−1, and −1.6934 kJ·mol−1, respectively, with Tyr365 and Tyr148 forming hydrophobic interactions mainly with Lys3me3 in the natural substrate and Trp364 mainly formed hydrophobic interactions with the hydrophobic portion on the Arg2 side chain of the natural substrate as well as the side chain methyl group of Ala5 (as shown in Figure 13).
Although ARLme3SA formed ionic interactions with proteins with a strength greater than or equal to that of hydrogen bonding and hydrophobic interactions, the simultaneous presence of many positively charged residues in EED protein caused their repulsive forces with ARLme3SA to offset a large portion of the electrostatic attraction. Therefore, the ΔGbinding of ARLme3SA with EED was only −139.14 kJ/mol, which was only lower than the weaker Inhibitors 2 and 5 among the six inhibitors, which provided an opportunity for the development of inhibitors targeting EED.

3.2. Differences in the Inhibitory Capacity of the Six Inhibitors

The six Inhibitors studied in our study could be grouped and compared to determine the effect of different parts of Inhibitors on inhibitory ability.
(i) The head and tail of Inhibitors 1,2 and 3 were identical, only the substituents on the 6-position of the parent nucleus were different, and their inhibitory capacities were more than 100 times different. Among them, Inhibitor 3 had the weakest inhibitory ability (IC50 = 13,000 nM), and the difference between Inhibitors 1 and 2 was only three times (IC50 = 110 nM and 320 nM). In the previous conformation comparison, in Inhibitor 3, the hydrogen atom on the head imino group was prone to flip due to the stronger electrostatic attraction of the cyano group on the 6-position of the parent nucleus. This flip resulted in a decrease in the hydrogen bond formed between it and Asn194, resulting in a considerable weakening of the binding of Inhibitor 3 to Asn194. In addition, the parent nucleus of Inhibitor 3 underwent a maximal vertical movement, which partially removed the parent nucleus and the trailing benzene ring from the optimal binding position with the aromatic cage residues (Phe97, Tyr148, and Tyr365), resulting in a weakening of the binding of the Inhibitor 3 to the aromatic cage residues. There was a significant reduction in the EMM contribution of Phe97. The above results indicate that in the design of this kind of drug, it is necessary to avoid adding strong negatively charged groups to the 6-position of the parent nucleus to prevent the reversal of the hydrogen atom on the head imino group.
The conformation comparison also showed that the parent nucleus of Inhibitor 2 moved horizontally compared to Inhibitor 1 due to the lack of fixation by Arg414. Unlike vertical movement, this horizontal movement did not cause the Inhibitor parent nucleus to leave the aromatic cage, so the binding of Inhibitor 2 to the aromatic cage residues was somewhat reduced but not as significant as in Inhibitor 3. At the same time, the horizontal movement resulted in a significant decrease in the distance between the 2-nitrogen atom on the parent nucleus triazole ring and the charged residues such as Lys211, Asp237, and Glu238, and the overall manifestation was an increase in the binding of Inhibitor 2 to the EED protein. The above interactions cancel each other, and the combined result made the binding free energy of Inhibitor 2 and the EED protein slightly smaller than that of Inhibitor 1, and the inhibition ability is slightly decreased.
(ii) The parent nucleus and tail of Inhibitors 1,4 and 5 were identical; only the head was different. Inhibitors 4 (IC50 = 20 nM) and 5 (IC50 = 2600 nM) formed a dihydrofuran ring on their head, but the location of the oxygen atom in the ring was different, making their inhibition ability different by 130 times. Previous results had shown that the hydrogen atom on the head imino in Inhibitor 5 was susceptible to flipping due to the intramolecular electrostatic attraction of the head oxygen atom. Like the previous Inhibitor 3, the flipping of the hydrogen atom resulted in a weakening of the binding of Inhibitor 5 to Asn194, as well as a vertical movement of the parent nucleus and tail of Inhibitor 5, and finally a decrease in the binding of Inhibitor 5 to the aromatic cage residues Tyr148 and Phe97. In addition, the difference in the position of the oxygen atom in the head also resulted in a reduced electrostatic attraction between Inhibitor 5 and Arg367 and Met366, which also reduced the affinity between the EED protein and Inhibitor 5.
Inhibitor 1 (IC50 = 110 nM) and Inhibitor 4 (IC50 = 20 nM) had essentially the same head oxygen atom position but had an open-loop structure in Inhibitor 1. A comparison of the conformations revealed that the head and parent nucleus positions of Inhibitors 1 and 4 overlap almost precisely, but there was a significant difference in the tails. This difference might be due to intramolecular interactions. That is, the position of the tailed benzene ring in Inhibitor 4 was closer to the dihydrofuran ring in the head, forming specific intramolecular hydrophobic interactions. To quantify the stability of the two tail orientations, we attached the benzodihydrofuran ring on the head of Inhibitor 3 to the head of Inhibitor 1. However, the tail of Inhibitor 1 remained unchanged, naming it Inhibitor 1′. We immobilized atoms with significant conformational changes and optimized Inhibitors 1′ and 3 on B3LYP/6-31G(d) using Gaussian 09. The results showed that the energy of Inhibitor 3 was 13.50 kJ·mol−1 lower than Inhibitor 1′, suggesting that intramolecular interactions made the conformation of Inhibitor 3 more stable.
The change in the tail conformation of the Inhibitor allowed the tail benzene ring to be fixed in a position that facilitated stable hydrophobic interaction with Tyr365, Trp364, and Tyr148 of the EED protein. As a result, each of these residues bonded more strongly in Inhibitor 3 than in Inhibitor 1. The above results suggest that in drug design, such Inhibitors should not carry a strong negatively charged group at the 8′-position of the head but preferably with a hydrophobic group to maintain a more favorable tail conformation.
(iii) The heads and parent nuclei of Inhibitors 4 (IC50 = 20 nM) and 6 (IC50 = 24 nM) were identical, but their tails differed. The difference in the binding positions of Inhibitor 4 and Inhibitor 6 in the EED protein was negligible. The EMM contributions showed that the binding of the two Inhibitors with the same parent nucleus to the aromatic cage residues differed little. The different Inhibitor tails mainly affected residues Phe97, Asn146, and Asp430, which were closer to it, as well as residues Arg414, Lys211, Glu238, and Asp237, which were slightly further away but with a charged side chain that could interact with the tail in a long-range electrostatic force. However, these residues positively and negatively affected the binding free energy between the EED and Inhibitor, essentially canceling each other. As a result, the difference in the total binding free energy of the two Complexes was slight, and the experimental results also showed a small difference in the inhibitory capacity of the two Inhibitors.

4. Materials and Methods

4.1. Acquisition of the EED Protein and Inhibitor Structures

The X-ray crystal diffraction structure (PDBID:7QK4, resolution: 1.6 A) [30] of the EED-MAK683 Complex was obtained from the RCSB database, and the solvent molecules and the original inhibitor MAK683 in the crystal were removed to get the EED protein. Complete the missing residues in the middle of the EED protein using the model loops panel in Chimera [37]. All structures of Inhibitors were constructed from Gaussview 6.0 and optimized at B3LYP/6-31G(d) level using The Gaussian 09 software package to obtain stable initial structure and Mulliken charge [38].

4.2. Molecular Docking

We used AutoDock Vina [39,40,41] as molecular docking software to obtain EED protein Complexes with the six Inhibitors we studied. Since these are all competitive Inhibitors [30], we selected the binding pocket of the natural substrate as the docking area, the aromatic cage area mentioned in the introduction. We used the built-in grid of AutoDock Vina to generate a grid box with 24 × 24 × 24 grid points in the x, y, and z directions with a grid spacing of 0.1 nm, and the center coordinates of the grid in the x, y, and z directions were 8.291, 17.934, and 17.401. Then, each of the six Inhibitors was molecular docked with the EED protein in the grid. AutoDock Vina outputted the top 20 site poses for each Inhibitor and selected the conformations that scored high and matched the binding modes reported in the literature as the initial structures for the next step of molecular dynamics simulations. The Complexes obtained from the six docking groups were named Complexes 1–6.

4.3. Molecular Dynamics Simulation

Taking the structures of the Complexes obtained by molecular docking as the initial conformations, we performed molecular dynamics simulations of all the complex systems using the Gromacs2018.8 program [42,43,44] As a reference, the same simulation was performed for the EED-ARLme3SA complex with natural substrates (PDBID:3IIW). The EED Protein selected the AMBER14SB-parmbsc1.ff force field [45]. The antechamber program in the AMBER16 [46] software package generated the GAFF force field for the Inhibitors [47]. Then, the TIP3P water molecule model [48] was used to solvate the entire complex systems. The shortest distance between the protein edge and the cubic solvent box boundary was 1 nm. The total number of water molecules added after the solvation of different complex systems was about 22,600. To maintain the electrical neutrality of the system, we added two sodium ions to neutralize the negative charge in the protein. The steepest descent method minimized the system’s energy by 1000 steps to ensure no atomic collisions or unreasonable spatial structure. Then the balance of two stages was carried out. The first stage was the NVT balance of 500 ps at 300 K temperature. In the second stage, the NPT balance of 500 ps was performed at the same temperature, and the pressure was set at 1atm. The final product simulations were carried out separately for 200 ns. The integration step was set to 2 fs, and the long-range electrostatic cut-off radius was 1.2 nm. The coordinates and energy information were saved every 10 ps during the simulations. Two sets of parallel simulations were performed for each system to rule out randomization of results.

4.4. Binding Free Energy Calculation

We used the Molecular Mechanics/Poisson–Boltzmann surface area (MM/PBSA) method to calculate the binding free energy based on the simulation results [36,49]. A point was selected every 20 frames from the equilibrium trajectory simulated by molecular dynamics, and 1000 frames of conformation were chosen with a duration of 20 ns. The binding free energy of each protein–ligand complex system was calculated using the g_mmpbsa tool in the Gromacs2018.8 software package, and the contribution of each residue to the decomposition of free energy was also estimated. The binding free energy of protein-ligand complexes is defined as follows [50]:
ΔGbind = ⟨Gcomplex⟩ − ⟨Gprotein⟩ − ⟨Gligand⟩
Gcomplex, Gprotein, and Gligand represent the free energy of the protein–ligand complex, protein, and ligand. The individual entity free energies of each of them can be further described as the following:
Gx = ⟨EMM⟩ − TS + ⟨Gsolv⟩
where x stands for complex, protein, and ligand, the free energy G is divided into two terms, namely molecular mechanical term and solvation energy. <EMM> represents the potential energy of molecules in a vacuum, <Gsolv> represents the free energy of solvation, and TS refers to the contribution of entropy to the free energy in a vacuum. In the Gromacs program, the TS term is approximately zero.
The vacuum molecular mechanical potential energy (EMM) is the sum of the energy of Ebonded and Enonbonded interactions, as shown in the following equation:
EMM = Ebonded + Enonbonded
Ebonded and Enonbonded can also be expressed as:
Ebonded = Ebond + Eangle + Etorsion
Enonbonded = Eelec + EvdW
where Ebond, Eangle, and Etorsion are represented as bonding interactions of bond length, bond angle, and dihedral angle, while Eelec and EvdW are represented as non-bonding interactions of the Van der Waals Force and the Coulomb Force and are modeled using the Coulomb Force potential function and Lennard–Jones potential function respectively.
The free energy of solvation (Gsolv) is composed of the polar free energy of solvation (Gpolar) and non-polar free energy of solvation (Gnonpolar), which can be expressed as the following formula:
Gsolv = Gpolar + Gnonpolar

5. Conclusions

In this work, molecular dynamics simulations and binding free energy calculations were performed on six triazolopyrimidine inhibitors and natural substrates in complexes with the EED subunit in PRC2. We mainly concluded the following: on the one hand, we compared the inhibitor to a positively charged natural substrate due to the high number of positively charged residues in the central pocket of the EED, making the binding of the natural substrate to the EED unstable and providing an energetic basis for competition for the inhibitor. On the other hand, we found that when the 6-position of the parent nucleus or the 8′-position of the head benzodihydrofuran had a strong negatively charged group, it caused the hydrogen atom on the head imino group to be susceptible to flipping and the parent nucleus moving vertically, resulting in a significant decrease in inhibitory ability. Furthermore, when the 6-position of the parent nucleus was a nonpolar group, its reduced electrostatic interaction with Arg414 resulted in the parent nucleus moving horizontally, which decreased the inhibitory ability slightly. We also discovered that when the head formed a benzodihydrofuran ring and the 8′-position was methylene, the inhibitory ability was significantly improved. Finally, we found that the differences in the tail of the inhibitor had little impact on its inhibitory ability.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules28247997/s1.

Author Contributions

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

Funding

This research was funded by Jilin Provincial Science and Technology Development Plan, grant number, 20220508114RC; Jilin University Graduate Innovation Research Program Project, grant number, 2022171, 2023CX184; Jilin Provincial Association for Science and Technology Domestic and Foreign Academic Exchange Project, grant number, KJLT202312; Jilin Province Higher Education Research Project, grant number, JGJX2022D31; Jilin Provincial Social Science Foundation, grant number, 2022Q3 and Jilin University Graduate Teaching Reform Project, grant number, 2022JGZ032.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article and Supplementary Materials.

Acknowledgments

We thank Yurou Zhang, Xinyue Zhang, Qinqin Li and Xuerui Zhao for their instruction.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Glancy, E.; Ciferri, C.; Bracken, A.P. Structural basis for PRC2 engagement with chromatin. Curr. Opin. Struct. Biol. 2021, 67, 135–144. [Google Scholar] [CrossRef]
  2. Conway, E.; Healy, E.; Bracken, A.P. PRC2 mediated H3K27 methylations in cellular identity and cancer. Curr. Opin. Cell Biol. 2015, 37, 42–48. [Google Scholar] [CrossRef] [PubMed]
  3. Hojfeldt, J.W.; Laugesen, A.; Willumsen, B.M.; Damhofer, H.; Hedehus, L.; Tvardovskiy, A.; Mohammad, F.; Jensen, O.N.; Helin, K. Accurate H3K27 methylation can be established de novo by SUZ12-directed PRC2. Nat. Struct. Mol. Biol. 2018, 25, 225–232. [Google Scholar] [CrossRef] [PubMed]
  4. Lavarone, E.; Barbieri, C.M.; Pasini, D. Dissecting the role of H3K27 acetylation and methylation in PRC2 mediated control of cellular identity. Nat. Commun. 2019, 10, 1679. [Google Scholar] [CrossRef]
  5. Zhao, Y.; Guan, Y.; Zhao, F.; Yu, T.; Zhang, S.; Zhang, Y.; Duan, Y.; Zhou, X. Recent strategies targeting Embryonic Ectoderm Development (EED) for cancer therapy: Allosteric inhibitors, PPI inhibitors, and PROTACs. Eur. J. Med. Chem. 2022, 231, 114–144. [Google Scholar] [CrossRef] [PubMed]
  6. Jürg, M.; Verrijzer, P. Biochemical mechanisms of gene regulation by polycomb group protein complexes. Curr. Opin. Genet. Dev. 2009, 19, 150–158. [Google Scholar]
  7. Ku, S.; Spencer, R.; Wang, Y.; Mu, P.; Seshadri, M.; Goodrich, Z.W.; Goodrich, M.M.; Labbé, D.P.; Gomez, E.C.; Wang, J.; et al. Rb1 and Trp53 cooperate to suppress prostate cancer lineage plasticity, metastasis, and antiandrogen resistance. Science 2017, 355, 78–83. [Google Scholar] [CrossRef]
  8. Sato, T.; Kaneda, A.; Tsuji, S.; Isagawa, T.; Yamamoto, S.; Fujita, T.; Yamanaka, R.; Tanaka, Y.; Nukiwa, T.; Marquez, V.E.; et al. PRC2 overexpression and PRC2-target gene repression relating to poorer prognosis in small cell lung cancer. Sci. Rep. 2013, 3, 1911. [Google Scholar] [CrossRef]
  9. Zhang, H.; Qi, J.; Reyes, J.M.; Li, L.; Rao, P.K.; Li, F.; Lin, C.Y.; Perry, J.A.; Lawlor, M.A.; Federation, A.; et al. Oncogenic Deregulation of EZH2 as an Opportunity for Targeted Therapy in Lung Cancer. Cancer Discov. 2016, 6, 1006–1021. [Google Scholar] [CrossRef]
  10. Puppe, J.; Opdam, M.; Schouten, P.C.; Jóźwiak, K.; Lips, E.; Severson, T.; Van de Ven, M.; Brambillasca, C.; Bouwman, P.; Van Tellingen, O.; et al. EZH2 Is Overexpressed in BRCA1-like Breast Tumors and Predictive for Sensitivity to High-Dose Platinum-Based Chemotherapy. Clin. Cancer Res. 2019, 25, 4351–4362. [Google Scholar] [CrossRef]
  11. Gao, S.B.; Zheng, Q.F.; Xu, B.; Pan, C.B.; Li, K.L.; Zhao, Y.; Zheng, Q.L.; Lin, X.; Xue, L.X.; Jin, G.H. EZH2 represses target genes through H3K27-dependent and H3K27-independent mechanisms in hepatocellular carcinoma. Mol. Cancer Res. 2014, 12, 1388–1397. [Google Scholar] [CrossRef]
  12. Gan, L.; Xu, M.; Hua, R.; Tan, C.; Zhang, J.; Gong, Y.; Wu, Z.; Weng, W.; Sheng, W.; Guo, W. The polycomb group protein EZH2 induces epithelial–mesenchymal transition and pluripotent phenotype of gastric cancer cells by binding to PTEN promoter. J. Hematol. Oncol. 2018, 11, 9. [Google Scholar] [CrossRef]
  13. Jones, B.A.; Varambally, S.; Arend, R.C. Histone Methyltransferase EZH2: A Therapeutic Target for Ovarian Cancer. Mol. Cancer Ther. 2018, 17, 591–602. [Google Scholar] [CrossRef] [PubMed]
  14. Wassef, M.; Margueron, R. The Multiple Facets of PRC2 Alterations in Cancers. J. Mol. Biol. 2017, 429, 1978–1993. [Google Scholar] [CrossRef] [PubMed]
  15. Højfeldt, L.A.; Westergaard, J.; Helin, K. Molecular Mechanisms Directing PRC2 Recruitment and H3K27 Methylation. Mol. Cell 2019, 74, 8–18. [Google Scholar]
  16. Dockerill, M.; Gregson, C.; O’Donovan, D.H. Targeting PRC2 for the treatment of cancer: An updated patent review. Expert Opin. Ther. Pat. 2020, 31, 119–135. [Google Scholar] [CrossRef]
  17. Kasinath, V.; Faini, M.; Poepsel, S.; Reif, D.; Feng, X.A.; Stjepanovic, G.; Aebersold, R.; Nogales, E. Structures of human PRC2 with its cofactors AEBP2 and JARID2. Science 2018, 359, 940–944. [Google Scholar] [CrossRef]
  18. Van Mierlo, G.; Veenstra, G.J.C.; Vermeulen, M.; Marks, H. The Complexity of PRC2 Subcomplexes. Trends Cell Biol. 2019, 29, 660–671. [Google Scholar] [CrossRef]
  19. McCabe, M.T.; Ott, H.M.; Ganji, G.; Korenchuk, S.; Thompson, C.; Van Aller, G.S.; Liu, Y.; Graves, A.P.; Della Pietra, A., 3rd; Diaz, E.; et al. EZH2 inhibition as a therapeutic strategy for lymphoma with EZH2-activating mutations. Nature 2012, 492, 108–112. [Google Scholar] [CrossRef]
  20. Knutson, S.K.; Kawano, S.; Minoshima, Y.; Warholic, N.M.; Huang, K.C.; Xiao, Y.; Kadowaki, T.; Uesugi, M.; Kuznetsov, G.; Kumar, N.; et al. Selective inhibition of EZH2 by EPZ-6438 leads to potent antitumor activity in EZH2-mutant non-Hodgkin lymphoma. Mol. Cancer Ther. 2014, 13, 842–854. [Google Scholar] [CrossRef]
  21. Italiano, A.; Soria, J.C.; Toulmonde, M.; Michot, J.M.; Lucchesi, C.; Varga, A.; Coindre, J.M.; Blakemore, S.J.; Clawson, A.; Suttle, B.; et al. Tazemetostat, an EZH2 inhibitor, in relapsed or refractory B-cell non-Hodgkin lymphoma and advanced solid tumors: A first-in-human, open-label, phase 1 study. Lancet Oncol. 2018, 19, 649–659. [Google Scholar] [CrossRef]
  22. Gibaja, V.; Shen, F.; Harari, J.; Korn, J.; Ruddy, D.; Saenz-Vash, V.; Zhai, H.; Rejtar, T.; Paris, C.G.; Yu, Z.; et al. Development of secondary mutations in wild-type and mutant EZH2 alleles cooperates to confer resistance to EZH2 inhibitors. Oncogene 2016, 35, 558–566. [Google Scholar] [CrossRef]
  23. Eich, M.L.; Athar, M.; Ferguson, J.E., 3rd; Varambally, S. EZH2-Targeted Therapies in Cancer: Hype or a Reality. Cancer Res. 2020, 80, 5449–5458. [Google Scholar] [CrossRef]
  24. Baker, T.; Nerle, S.; Pritchard, J.; Zhao, B.; Rivera, V.M.; Garner, A.; Gonzalvez, F. Acquisition of a single EZH2 D1 domain mutation confers acquired resistance to EZH2-targeted inhibitors. Oncotarget 2015, 6, 32646–32655. [Google Scholar] [CrossRef]
  25. Margueron, R.; Justin, N.; Ohno, K.; Sharpe, M.L.; Son, J.; Drury, W.J., 3rd; Voigt, P.; Martin, S.R.; Taylor, W.R.; De Marco, V.; et al. Role of the polycomb protein EED in the propagation of repressive histone marks. Nature 2009, 461, 762–767. [Google Scholar] [CrossRef] [PubMed]
  26. Qi, W.; Zhao, K.; Gu, J.; Huang, Y.; Wang, Y.; Zhang, H.; Zhang, M.; Zhang, J.; Yu, Z.; Li, L.; et al. An allosteric PRC2 inhibitor targeting the H3K27me3 binding pocket of EED. Nat. Chem. Biol. 2017, 13, 381–388. [Google Scholar] [CrossRef]
  27. Jain, B.P.; Pandey, S. WD40 Repeat Proteins: Signalling Scaffold with Diverse Functions. Protein J. 2018, 37, 391–406. [Google Scholar] [CrossRef] [PubMed]
  28. Sanulli, S.; Justin, N.; Teissandier, A.; Ancelin, K.; Portoso, M.; Caron, M.; Michaud, A.; Lombard, B.; Da Rocha, S.T.; Offer, J.; et al. Jarid2 Methylation via the PRC2 Complex Regulates H3K27me3 Deposition during Cell Differentiation. Mol. Cell 2015, 57, 769–783. [Google Scholar] [CrossRef] [PubMed]
  29. Han, Z.; Xing, X.; Hu, M.; Zhang, Y.; Liu, P.; Chai, J. Structural basis of EZH2 recognition by EED. Structure 2007, 15, 1306–1315. [Google Scholar] [CrossRef]
  30. Huang, Y.; Sendzik, M.; Zhang, J.; Gao, Z.; Sun, Y.; Wang, L.; Gu, J.; Zhao, K.; Yu, Z.; Zhang, L.; et al. Discovery of the Clinical Candidate MAK683: An EED-Directed, Allosteric, and Selective PRC2 Inhibitor for the Treatment of Advanced Malignancies. J. Med. Chem. 2022, 65, 5317–5333. [Google Scholar] [CrossRef] [PubMed]
  31. Lingel, A.; Sendzik, M.; Huang, Y.; Shultz, M.D.; Cantwell, J.; Dillon, M.P.; Fu, X.; Fuller, J.; Gabriel, T.; Gu, J.; et al. Structure-Guided Design of EED Binders Allosterically Inhibiting the Epigenetic Polycomb Repressive Complex 2 (PRC2) Methyltransferase. J. Med. Chem. 2017, 60, 415–427. [Google Scholar] [CrossRef]
  32. Huang, Y.; Zhang, J.; Yu, Z.; Zhang, H.; Wang, Y.; Lingel, A.; Qi, W.; Gu, J.; Zhao, K.; Shultz, M.D.; et al. Discovery of First-in-Class, Potent, and Orally Bioavailable Embryonic Ectoderm Development (EED) Inhibitor with Robust Anticancer Efficacy. J. Med. Chem. 2017, 60, 2215–2226. [Google Scholar] [CrossRef] [PubMed]
  33. Li, G.; Tang, Z.; Fang, H.; Shi, W. Inhibition of hepatocellular carcinoma cell proliferation by embryonic ectodermal development protein inhibitor and its molecular mechanism. Chin. J. Exp. Surg. 2021, 38, 648–650. [Google Scholar]
  34. Rej, R.K.; Wang, C.; Lu, J.; Wang, M.; Petrunak, E.; Zawacki, K.P.; McEachern, D.; Fernandez-Salas, E.; Yang, C.Y.; Wang, L.; et al. EEDi-5285: An Exceptionally Potent, Efficacious, and Orally Active Small-Molecule Inhibitor of Embryonic Ectoderm Development. J. Med. Chem. 2020, 63, 7252–7267. [Google Scholar] [CrossRef]
  35. Rej, R.K.; Wang, C.; Lu, J.; Wang, M.; Petrunak, E.; Zawacki, K.P.; McEachern, D.; Yang, C.Y.; Wang, L.; Li, R.; et al. Discovery of EEDi-5273 as an Exceptionally Potent and Orally Efficacious EED Inhibitor Capable of Achieving Complete and Persistent Tumor Regression. J. Med. Chem. 2021, 64, 14540–14556. [Google Scholar] [CrossRef]
  36. Homeyer, N.; Gohlke, H. Free Energy Calculations by the Molecular Mechanics Poisson-Boltzmann Surface Area Method. Mol. Inf. 2012, 31, 114–122. [Google Scholar] [CrossRef] [PubMed]
  37. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera—A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605–1612. [Google Scholar] [CrossRef] [PubMed]
  38. Tirado-Rives, J.; Jorgensen, W.L. Performance of B3LYP Density Functional Methods for a Large Set of Organic Molecules. J. Chem. Theory Comput. 2008, 4, 297–306. [Google Scholar] [CrossRef]
  39. Trott, O.; Olson, A.J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2010, 31, 455–461. [Google Scholar] [CrossRef]
  40. Tanchuk, V.Y.; Tanin, V.O.; Vovk, A.I.; Poda, G. A New, Improved Hybrid Scoring Function for Molecular Docking and Scoring Based on AutoDock and AutoDock Vina. Chem. Biol. Drug Des. 2016, 87, 618–625. [Google Scholar] [CrossRef]
  41. Wang, Z.; Sun, H.; Yao, X.; Li, D.; Xu, L.; Li, Y.; Tian, S.; Hou, T. Comprehensive evaluation of ten docking programs on a diverse set of protein-ligand complexes: The prediction accuracy of sampling power and scoring power. Phys. Chem. Chem. Phys. 2016, 18, 12964–12975. [Google Scholar] [CrossRef] [PubMed]
  42. Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A.E.; Berendsen, H.J. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701–1718. [Google Scholar] [CrossRef]
  43. Kutzner, C.; Páll, S.; Fechner, M.; Esztermann, A.; De Groot, B.L.; Grubmüller, H. Best bang for your buck: GPU nodes for GROMACS biomolecular simulations. J. Comput. Chem. 2015, 36, 1990–2008. [Google Scholar] [CrossRef] [PubMed]
  44. Abraham, M.J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.C.; Hess, B.; Lindahl, E. High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19–25. [Google Scholar] [CrossRef]
  45. Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. [Google Scholar] [CrossRef]
  46. Dickson, C.J.; Rosso, L.; Betz, R.M.; Walker, R.C.; Gould, I.R. GAFFlipid: A General Amber Force Field for the accurate molecular dynamics simulation of phospholipid. Soft Matter 2012, 8, 9617–9627. [Google Scholar] [CrossRef]
  47. Case, D.A.; Cheatham, T.E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K.M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R.J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668–1688. [Google Scholar] [CrossRef]
  48. Mark, P.; Nilsson, L.J. Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. Phys. Chem. A 2001, 105, 9954–9960. [Google Scholar] [CrossRef]
  49. Kumari, R.; Kumar, R. g_mmpbsa—A GROMACS tool for high-throughput MM-PBSA calculations. J. Chem. Inf. Model. 2014, 54, 1951–1962. [Google Scholar] [CrossRef]
  50. Valdes-Tresanco, M.S.; Valdes-Tresanco, M.E.; Valiente, P.A.; Moreno, E. gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS. Theory Comput. 2021, 17, 6281–6291. [Google Scholar] [CrossRef]
Molecules 28 07997 g001
Figure 1. Conformation of the core subunit in PRC2 protein (PDBID:6c24).
Figure 1. Conformation of the core subunit in PRC2 protein (PDBID:6c24).
Molecules 28 07997 g001
Molecules 28 07997 g002
Figure 2. Conformations of the EED protein. (a) The EED top surface contains H3K27me3 (PDBID:3IIW), where the four light blue residues represent aromatic cage residues and the pink portion is the natural substrate (b) and the EED bottom surface contains the N-terminal structural domain of the EZH2 subunit (PDBID:7QK4).
Figure 2. Conformations of the EED protein. (a) The EED top surface contains H3K27me3 (PDBID:3IIW), where the four light blue residues represent aromatic cage residues and the pink portion is the natural substrate (b) and the EED bottom surface contains the N-terminal structural domain of the EZH2 subunit (PDBID:7QK4).
Molecules 28 07997 g002
Molecules 28 07997 g003
Figure 3. The structures of the six Inhibitors include the ring’s atom numbering and corresponding IC50 value.
Figure 3. The structures of the six Inhibitors include the ring’s atom numbering and corresponding IC50 value.
Molecules 28 07997 g003
Molecules 28 07997 g004
Figure 4. The binding position of inhibitors and the formation of hydrogen bonds in EED. (a,b) The binding positions of the Inhibitors in the EED, and the circular magnification frame detail the comparison of the poses of the six Inhibitors, where the atoms marked in dark red are oxygen atoms. (c) The square magnification frame shows the hydrogen bonding formed between the Inhibitors and the surrounding residues after docking (in the case of inhibitor 2, which had the highest number of hydrogen bonds), where the red dotted lines represent hydrogen bonds and the portion labelled green is the residues that form hydrogen bonds with Inhibitor 2.
Figure 4. The binding position of inhibitors and the formation of hydrogen bonds in EED. (a,b) The binding positions of the Inhibitors in the EED, and the circular magnification frame detail the comparison of the poses of the six Inhibitors, where the atoms marked in dark red are oxygen atoms. (c) The square magnification frame shows the hydrogen bonding formed between the Inhibitors and the surrounding residues after docking (in the case of inhibitor 2, which had the highest number of hydrogen bonds), where the red dotted lines represent hydrogen bonds and the portion labelled green is the residues that form hydrogen bonds with Inhibitor 2.
Molecules 28 07997 g004
Molecules 28 07997 g005
Figure 5. RMSD diagram of the Cα atomic skeleton of six complex systems and EED-ARLme3SA with simulated time.
Figure 5. RMSD diagram of the Cα atomic skeleton of six complex systems and EED-ARLme3SA with simulated time.
Molecules 28 07997 g005
Molecules 28 07997 g006
Figure 6. Free energy landscape and sampling of seven complex systems.
Figure 6. Free energy landscape and sampling of seven complex systems.
Molecules 28 07997 g006
Molecules 28 07997 g007
Figure 7. Complexes overlap the diagram of six samples, where the atoms marked in dark red are oxygen atoms. The square magnification frame is the overlap of six Inhibitors, and the red rectangular frame shows the location of the parent nucleus in Inhibitors 1, 4, and 6.
Figure 7. Complexes overlap the diagram of six samples, where the atoms marked in dark red are oxygen atoms. The square magnification frame is the overlap of six Inhibitors, and the red rectangular frame shows the location of the parent nucleus in Inhibitors 1, 4, and 6.
Molecules 28 07997 g007
Molecules 28 07997 g008
Figure 8. Intramolecular interactions diagrams and dihedral angle diagrams. (a) The hydrogen atom on the head imino group of the Inhibitors reversal and intramolecular interactions, where green represents the Inhibitor 1, yellow represents the Inhibitor 3, light blue represents the Inhibition 5 and red dashed lines represent intramolecular interactions. (b) The change of dihedral angle H1-N2-C3-N4 of the six complex systems with the simulation time during the whole simulation process.
Figure 8. Intramolecular interactions diagrams and dihedral angle diagrams. (a) The hydrogen atom on the head imino group of the Inhibitors reversal and intramolecular interactions, where green represents the Inhibitor 1, yellow represents the Inhibitor 3, light blue represents the Inhibition 5 and red dashed lines represent intramolecular interactions. (b) The change of dihedral angle H1-N2-C3-N4 of the six complex systems with the simulation time during the whole simulation process.
Molecules 28 07997 g008
Molecules 28 07997 g009
Figure 9. (a) Schematic calculation of Gaussian optimization for Inhibitor 3–Inhibitor 3′. (b) Schematic calculation of Gaussian optimization for Inhibitor 5–Inhibitor 5′.
Figure 9. (a) Schematic calculation of Gaussian optimization for Inhibitor 3–Inhibitor 3′. (b) Schematic calculation of Gaussian optimization for Inhibitor 5–Inhibitor 5′.
Molecules 28 07997 g009
Molecules 28 07997 g010
Figure 10. Interaction between aromatic cage residues and Inhibitors, where subfigures (a–c) are residues Tyr148, Trp364, Tyr365 interacting with the Inhibitor and subfigures (d–f) are residues Phe97 interacting with the inhibitor. The yellow dashed lines represent π-π stacking interaction, and the grey dashed lines represent hydrophobic interactions.
Figure 10. Interaction between aromatic cage residues and Inhibitors, where subfigures (a–c) are residues Tyr148, Trp364, Tyr365 interacting with the Inhibitor and subfigures (d–f) are residues Phe97 interacting with the inhibitor. The yellow dashed lines represent π-π stacking interaction, and the grey dashed lines represent hydrophobic interactions.
Molecules 28 07997 g010
Molecules 28 07997 g011
Figure 11. Locations of the four residues interacting with the parent nucleus of the Inhibitors. The green dashed lines represent salt bridges, the blue dashed lines represent the electrostatic attraction, the pink dashed lines represent the electrostatic repulsion and the distances in the figure represent average distances during the simulation.
Figure 11. Locations of the four residues interacting with the parent nucleus of the Inhibitors. The green dashed lines represent salt bridges, the blue dashed lines represent the electrostatic attraction, the pink dashed lines represent the electrostatic repulsion and the distances in the figure represent average distances during the simulation.
Molecules 28 07997 g011
Molecules 28 07997 g012
Figure 12. Locations of the six residues interacting with the head and tail of the Inhibitors. The red dashed lines represent the hydrogen bond, the yellow dashed lines represent π-π stacking interaction, the gray dashed lines represent hydrophobic interactions, the blue dashed lines represent the electrostatic attraction, the pink dashed lines represent the electrostatic repulsion and the distances in the figures represent average distances during the simulation.
Figure 12. Locations of the six residues interacting with the head and tail of the Inhibitors. The red dashed lines represent the hydrogen bond, the yellow dashed lines represent π-π stacking interaction, the gray dashed lines represent hydrophobic interactions, the blue dashed lines represent the electrostatic attraction, the pink dashed lines represent the electrostatic repulsion and the distances in the figures represent average distances during the simulation.
Molecules 28 07997 g012
Molecules 28 07997 g013
Figure 13. Comparison of ARLme3SA and Inhibitor binding poses and interaction of ARLme3SA with its surrounding residues. The green portion is ARLme3SA and the gold portion is the Inhibitor 4. The light blue residues represent the negatively charged region, the pink residues represent the positively charged region, and the yellow residues represent the aromatic cage.
Figure 13. Comparison of ARLme3SA and Inhibitor binding poses and interaction of ARLme3SA with its surrounding residues. The green portion is ARLme3SA and the gold portion is the Inhibitor 4. The light blue residues represent the negatively charged region, the pink residues represent the positively charged region, and the yellow residues represent the aromatic cage.
Molecules 28 07997 g013
Table 1. Occurrence of vital residues in protein forming hydrogen bonds with Inhibitors.
Table 1. Occurrence of vital residues in protein forming hydrogen bonds with Inhibitors.
ComplexBonding ResidueOccurrence (%)
Complex 1Asn19485.34
ARG3678.47
ASN1461.33
Complex 2ASN19498.24
ARG3678.34
LYS2116.06
Complex 3ASN19437.28
ARG41419.18
ARG3674.14
Complex 4ASN194108.51
ARG3676.65
LYS2113.70
Complex 5ASN19463.60
ARG4144.57
LYS2114.10
Complex 6ASN194104.13
ASN14617.61
ARG3676.12
Table 2. Binding free energy of six compounds calculated by MM/PBSA (kJ·mol−1).
Table 2. Binding free energy of six compounds calculated by MM/PBSA (kJ·mol−1).
Energy
Contribution		Complex 1Complex 2Complex 3Complex 4Complex 5Complex 6
ΔEvdW−209.73−192.71−150.48−227.44−176.55−219.95
ΔEelec−76.93−79.11−45.68−83.26−37.68−93.09
ΔGPB151.89145.65100.42154.29111.65159.39
ΔGSA−17.66−16.22−12.07−16.87−13.74−17.34
ΔGpolar a74.9666.5454.7471.0373.9766.30
ΔGnonpolar b−227.39−208.93−162.55−244.31−190.29−237.29
ΔGbinding c−152.43−142.39−107.81−173.28−116.32−170.99
a ΔGpolar = ΔGelec + ΔGPB, b ΔGnonpolar = ΔGvdW + ΔGSA, c ΔGbind = ΔGpolar + ΔGnonpolar.
Table 3. EMM contributions of significant residues (kJ·mol−1).
Table 3. EMM contributions of significant residues (kJ·mol−1).
ResidueComplex 1Complex 2Complex 3Complex 4Complex 5Complex 6
Tyr365−22.04−25.58−21.86−27.71−22.18−27.76
Tyr148−16.80−14.11−13.19−17.36−12.77−18.74
Phe97−16.03−12.24−8.53−14.50−10.94−17.06
Trp364−1.38−2.21−1.26−3.67−2.02−1.74
Arg367−8.80−4.83−7.61−9.99−0.99−7.65
Asn194−20.71−20.27−6.99−22.94−11.15−23.50
Lys211−5.56−21.37−5.24−11.64−12.93−5.90
Glu238−3.294.28−0.84−0.925.41−5.84
Asp2370.532.650.101.021.37−0.20
Arg414−3.53−1.91−6.86−2.80−5.63−5.02
Asp310−3.54−4.21−3.34−3.34−10.85−4.37
Met366−3.61−4.02−2.35−5.83−1.05−5.59
Asn146−4.68−2.39−2.84−1.71−1.55−4.15
Asp430−1.25−1.680.20−3.95−0.891.50
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.

Share and Cite

MDPI and ACS Style

Ju, J.; Zhang, H.; Guan, S.; Liu, C.; Du, J.; Shen, X.; Wang, S. Insight into the Inhibitory Mechanism of Embryonic Ectoderm Development Subunit by Triazolopyrimidine Derivatives as Inhibitors through Molecular Dynamics Simulation. Molecules 2023, 28, 7997. https://doi.org/10.3390/molecules28247997

AMA Style

Ju J, Zhang H, Guan S, Liu C, Du J, Shen X, Wang S. Insight into the Inhibitory Mechanism of Embryonic Ectoderm Development Subunit by Triazolopyrimidine Derivatives as Inhibitors through Molecular Dynamics Simulation. Molecules. 2023; 28(24):7997. https://doi.org/10.3390/molecules28247997

Chicago/Turabian Style

Ju, Jianan, Hao Zhang, Shanshan Guan, Chang Liu, Juan Du, Xiaoli Shen, and Song Wang. 2023. "Insight into the Inhibitory Mechanism of Embryonic Ectoderm Development Subunit by Triazolopyrimidine Derivatives as Inhibitors through Molecular Dynamics Simulation" Molecules 28, no. 24: 7997. https://doi.org/10.3390/molecules28247997

APA Style

Ju, J., Zhang, H., Guan, S., Liu, C., Du, J., Shen, X., & Wang, S. (2023). Insight into the Inhibitory Mechanism of Embryonic Ectoderm Development Subunit by Triazolopyrimidine Derivatives as Inhibitors through Molecular Dynamics Simulation. Molecules, 28(24), 7997. https://doi.org/10.3390/molecules28247997

Article Metrics

Back to TopTop