Docking and 3D-QSAR Studies on the Imidazo[1,5-c]pyrimidine Derivative as EED Inhibitors①

CHEN Xiao-Zhong LI Guang-Ping SHEN Yan HU Yong WANG Juan WANG Yuan-Qiang LIN Zhi-Hua

(Department of Pharmacy and Bioengineering, Chongqing University of Technology, Chongqing 400054, China)

ABSTRACT Embryonic ectoderm development (EED) has become a novel target for cancer treatment. In this study, a series of EED inhibitors was subjected to a three-dimensional quantitative structure-activity relationship(3D-QSAR) and molecular docking. Accordingly, this is the first of such 3D-QSAR studies in a series of EED inhibitors displaying anti-cancer pharmacological profiles. The CoMFA (q2 = 0.792, r2 = 0.994, r2pred = 0.74) and CoMSIA (q2 = 0.873, r2 = 0.994, r2pred = 0.81) models demonstrated good robustness and predictive ability.Moreover, molecular docking suggested that cation-?,?-??stacking and hydrogen bonding interactions were the main factors affecting the activity of these inhibitors. Five new small molecules were designed based on the CoMFA and CoMSIA contour maps. These molecules were then submitted to further ADME studies, in which the ADME properties of the five designed molecules were found to be within a reasonable range. In view of the corresponding findings, this study may provide theoretical guidance for the rational design of novel EED inhibitors.

Keywords: embryonic ectoderm development, CoMFA, CoMSIA, molecular docking;

1 INTRODUCTION

The polycomb repressive complex 2 (PRC2) consists of three proteins: EZH2 (an enhancer of zeste homolog 2), EED(embryonic ectoderm development), and SUZ12 (suppressor of zeste 12). PRC2 plays an essential role in regulating gene expression, and dysregulation of PRC2 is observed in multiple human cancers, such as lymphoblastic leukemia and malignant peripheral nerve sheath tumors[1]. Therefore,inhibiting PRC2 activity has become a strategy in the treatment of cancer. In this regard, one approach is developing EZH2 inhibitors that target the S-adenosyl-Lmethionine (SAM) binding site of EZH2[2]. Tazemetosta, an EZH2 inhibitor, has been approved by the FDA in January 2020. However, acquired drug resistance in EZH2 mutations is of concern[3]. Consequently, small molecule targeting binding sites of embryonic ectoderm development (EED)offers a novel therapeutic strategy in inhibiting PRC2 activity[4]. Novartis researchers have reported that EED226(Fig. 1), an EED inhibitor, binds to the histone-binding pocket of EED with triazolo[4,3-c]pyrimidine scaffold.Following a high-throughput screen against EED, A-395, a potent EED inhibitor was identified by AbbVie researchers[5].Using EED226 as a starting point, Rohan et al identified EEDi-5285, which binds to EED with anIC50 value of 0.2 nmol[6]. Promising EDD inhibitors such as EED226, A-395,and EEDi-5285 have antitumor activity in EZH2 inhibitor-resistant models. Fig. 1 illustrates various representative EED inhibitors such as EED226, A-395, and BR-001 in preclinical studies.

The 3D-QSAR study consisted of a comparative molecular field analysis (CoMFA)[7]and comparative molecular similarity indices analysis (CoMSIA)[8]. CoMFA and CoMSIA are two highly cited methods used to investigate the structure-activity relationships of a series of compounds in drug design. Molecular docking is a methodology that characterizes the binding modes between a receptor and its ligands.The combination of 3D-QSAR and molecular docking may demonstrate the detailed pharmacophoric features of the compounds and put forward guidelines for rational new compounds. This combination is effective for the discovery of potential novel drugs.

There are no existing reports on the 3D-QSAR studies on EED inhibitors. Consequently, in this study, CoMFA,CoMSIA, and molecular docking are used to study a set of novel EED inhibitors with imidazo[1,5-c]pyrimidine scaffold.According to the CoMFA and CoMSIA contour maps, five new molecules with high predictive activity were designed.In addition, an ADME study was carried out on the designed compounds. Accordingly, this study may be applied to the discovery of novel EED inhibitors.

Fig. 1. Representative EED inhibitors and EZH2 inhibitors

2 EXPERIMENTAL

2. 1 Dataset

A series of 34 EED inhibitors which have been recently described[6]was adopted for this research. Among these compounds, theIC50value of the most active compound 22 was 0.2 nmol, while theIC50value of the least active compound 10 was 207 nmol. All inhibitors were built using the SYBYL.2.1 software package. The structures of the inhibitors and their pIC50(–logIC50) values are listed in Table 1. A proper minimization of the molecule is vital for obtaining accurate 3D-QSAR models. All compounds were minimized using the Tripos force field[9]and Powell methods.In addition, the maximum iteration was set to 1000, with the other parameters as default values. Gasteiger-Hückel charges[10]were then assigned to each molecule, in which the training and test sets had 27 and 7 compounds, respectively.Test set compounds were marked with *, as shown in Table 1.

2. 2 Molecular alignment

All inhibitors were aligned with the universal chemical structure and compound 22 with the highest activity as the template molecule[11]. The superimposed common skeleton was imidazo[1,5-c]pyrimidine, and the 3Dstructure of the aligned compounds is shown in Fig. 2.

Table 1. Molecular Structures of Compounds and Their Observed pIC50 and Predicted pIC50 from CoMFA and CoMSIA Models

?

?

?

?

Fig. 2. Alignments of 35 imidazo[1,5-c]pyrimidine EED inhibitors

2. 3 Development CoMFA and CoMSIA models

CoMFA and CoMSIA models were built using the SYBYL.2.1.1 QSAR module. Lennard-Jones and the Coulomb potentials were then used to calculate the CoMFA steric and electrostatic fields, respectively[12]. Ansp3carbon atom probe with a radius of 1.53 ? and a charge of +1.0 was used to calculate the steric and electrostatic energies between the probe and molecules[13]. CoMSIA descriptors have increased hydrophobic field, hydrogen bond-acceptor field,and hydrogen bond-donor field compared to CoMFA[14].

The CoMFA and CoMSIA models were built using the partial least squares (PLS) method, which generates a linear relationship between molecular fields with its activity[15].CoMFA and CoMSIA descriptors and EED inhibitors pIC50 served as independent and dependent variables in thePLSmethod, respectively.LOO(Leave-One-Out) cross-validation was used to obtain the optimal number of components(ONC)[16]and the highest cross-validated coefficientq2. The optimal number of components used to perform the final PLS analysis. The model's internal predictive capability was evaluated byq2andr2, and its external predictive capability was evaluated byr2pred. The values ofq2,r2, andr2predwere determined by equations (1)～(3), respectively.

2. 4 Molecular docking

Selecting a docking program and modifying its parameters are critical in the docking results. Surflex has higher accuracy compared to other docking programs such as DOCK, FlexX,ICM, and PhDOCK[17]. Therefore, a docking study was completed in the surflex module of the SYBYL-X2.1.1 software. To obtain accurate docking results, GeomX mode was used to perform docking. The X-ray crystal structure of EED complexed with compound 22 (PDB code: 6W7F) that was resolved at 2.2 ? was downloaded from the protein data bank (https://www.rcsb.org). The water molecules were then deleted, and the original ligands were extracted from the crystal structure, after which the protomol generated base on the original ligands. The Sanitize protocol was selected as the ligand preparation method. The protein was allowed to move,and 50 conformations per ligands were produced in the process of molecular docking. The pose of ligands with the highest total scores served as the docking results.

2. 5 AMDE prediction

It is costly to predict the ADME (absorption, distribution,metabolism, and excretion) properties of compounds. A widely used method is to evaluate the ADME of a compound with computer tools prior to synthesis. Thus, the online free web tool SwissADME[18]was used to predict the ADME of newly designed EED inhibitors.

3 RESULTS AND DISCUSSION

3. 1 CoMFA and CoMSIA statistical analyses

The statistical parameters of CoMFA model are shown in Table 2. The steric and electrostatic fields contributed 53.5%and 46.5%, respectively. The cross-validated coefficientq2,no cross-validated coefficientr2and predictive correlation coefficientr2predwere 0.792, 0.994 and 0.740, respectively.The model exhibited excellent robustness and predictive ability based on the general evaluation criteria:q2> 0.5 andr2pred> 0.6[19]. Moreover, to make a satisfactory 3D-QSAR model,r2–q2should not be more than 0.3[20]. Here, the value ofr2–q2in the CoMFA model was 0.204, suggesting the model is acceptable. The predicted and experimental pIC50of the 34 EED inhibitors are listed in Table 1, and the scatter plots of the predicted and experimental pIC50values for the CoMFA model are depicted in Fig. 3.

To obtain a better CoMSIA model, the five fields were freely combined to generate 33 different CoMSIA models.The first parameter that evaluates the statistical robustness of a QSAR model is the value ofq2. The results of the top five highestq2models are listed in Table 3. Here, CoMSIA-SHAand CoMSIA-SHDAwere observed to have higherq2values than that of CoMSIA-SEHDA. In a satisfactory CoMFA or CoMSIA model, theONCshould be less than one-third of the total number of compounds studied. Thus, CoMSIA-SEHDA(ONC＜ 9) was chosen as the best CoMSIA model for further analysis. The statistical parameters of the CoMSIA-SEHDAmodel are shown in Table 2, and the statistical parameters of the CoMSIA-SEHDAmodel were observed to be better than the CoMFA model. The differences of statistical parameters of the two models can be attributed to the energy functions used to calculate the field values. In CoMFA, the energy function is very sensitive with the changes in position. In CoMSIA, these fields are calculated using much smoother potentials which are not as steep as the Lennard-Jones and Coulombic functions and have a finite value even at the atomic positions. The contributions of steric field, electrostatic field, hydrophobic field, hydrogen-bond acceptor field and hydrogen-bond donor field are 7.3%,20.6%, 30.1%, 41.1% and 0.9% in the CoMSIA-SEHDAmodel, respectively, indicating that the hydrogen-bond acceptor group has an essential influence on affinity. The pIC50values predicted by the CoMSIA-SEHDAmodel are given in Table 1. Additionally, the corresponding scatter plots of predicted and experimental pIC50value for the CoMSIA analyses are shown in Fig. 3.

Fig. 3. Plots of experimental pIC50 versus predicted pIC50 for the training and test set compounds by (a) CoMFA and (b) CoMSIA models

Table 2. Statistical Parameters of CoMFA and CoMSIA Analyses

Table 3. Different Field Combinations of the Top Five Highest q2 in CoMSIA Analysis

3. 2 CoMFA contour maps

The CoMFA contour map is shown in Fig. 4, in which the green region (80% contribution) represents the favorable steric region, while the yellow region (20% contribution)represents the unfavorable steric region. A large green contour was found to be located in the R1 position,suggesting that a bulky group was favorable in this region.For example, compound 2 has a bulky ethyl ester group in the R1 position, which has a higher affinity than compound 1 (R1= H) in binding EED. The blue region (80% contribution)refers to the area where the electron-withdrawing group was favorable, while the red region (20% contribution) represents the area where electron-donating group was favorable.Additionally, a large blue region is present in the R1 position,signifying that electron-withdrawing groups in this region would increase activity. This may explain why the activity of compound 15 (pIC50= 8.22) is higher than that of compound 12 (pIC50= 7.80). Moreover, a red region was observed near cyclopropyl in the R3 position, suggesting that electronwithdrawing groups here decreased activity. This may serve as a possible reason as to why the activity of compound 25 contained a strong electron-withdrawing group in the R3 position 4 times less than that of compound 22. Furthermore,blue and red contours were present in the R2 position,demonstrating that both the electron-withdrawing and electron-donating groups in this area uniformly affect its activity.

Fig. 4. Contour map of CoMFA model: green displays sterically favored regions, yellow is sterically disfavored regions.Blue contours indicate the area where electropositive substituents are favorable, and red region represents the area where electronegative substituents are favorable. Compound 22 with the highest activity is displayed as a reference

3. 3 CoMSIA contour maps

The CoMSIA contour maps are shown in Fig. 5, where compound 22 is found to be superimposed for reference. The steric and electrostatic contour maps of CoMSIA (Fig. 5d)were equivalent to that of the CoMFA model. Fig. 5a displays the hydrophobic contour map of the CoMSIA model. Here,the yellow region (80% contribution) and white region (20%contribution) signify the areas where hydrophobic and hydrophilic groups are favored and unfavored, respectively.One white region was observed to be located in the R1 position, demonstrating that compound 2 with an ethyl ester group in the R1 position has greater activity than compound 1.The other white regions were found near the pyridine ring of the R3 position, revealing that the hydrophobic substituent,like methyl, on the pyridine ring acted to decrease the activity,which may exhibit the following affinity order: compound 19< 25 and compound 27 < 24.

The CoMSIA contour map of the hydrogen-bond acceptor was plotted, as shown in Fig. 5b. The hydrogen-bond acceptor group in magenta (80% contribution) was found to enhance the activity. In contrast, the red areas (20%contribution) represent regions where hydrogen-bond acceptor groups were unfavorable. The large magenta color contour was observed to be in the R1 position, implying that the methylsulfonyl group in the R1 position may have hydrogen-bond interactions with EED. Moreover, a large red contour was located in the R2 position, suggesting that the hydrogen-bond accepter on this part decreased activity. For example, compound 9 (pIC50= 5.66) was less potent than compound 11 (pIC50= 7.35).

Fig. 5c shows the hydrogen-bond donor interaction in the CoMSIA model. The purple contours represent regions where a hydrogen-bond donor was observed to be unfavorable.Additionally, a large purple contour was found in the R1 position, suggesting that the hydrogen-bond donor groups in the R1 position were responsible for the decrease in activity.

Fig. 5. CoMSIA contour map with the combination of compound 22: (a) Hydrophobic contour maps. Yellow and white are the favorable and unfavorable regions. (b) Hydrogen-bond acceptor contour maps. Magenta and red represent the favorable and unfavorable regions.(c) Hydrogen-bond donor contour maps. Cyan and purple are the favorable and unfavorable regions. (d) Steric and electrostatic contour maps. Green displays the sterically favored regions, yellow is the sterically disfavored regions; Blue represents the electropositive favored regions, and red shows the electronegative favored regions

3. 4 Molecular docking

To verify the docking protocol, the conformation of compound 28 extracted from the X-ray structure was compared with its re-docked conformation. The two poses were superimposed, in which their RMSD was found to be 0.887 ?, indicating that the protocol is acceptable. To better explain the binding modes between EED and its inhibitors,molecular docking was performed for the highest active compound 22 as well as for the lower active compounds 1 and 10. The docking total scores of compounds 28, 7, and 16 were 13.34, 10.77, and 9.84, respectively. The three compounds were located at the active site of EED, and the key interactions between these three compounds and EED are given in Fig. 6. Compound 22 formed hydrogen bonding interactions with amino acid residues ASN194 and LYS211,consistent with that of compounds 10 and 1 (Fig. 6).Moreover, the fluorine atom of compound 22 had a halogen interaction with EED, whereas compounds 10 and 1 did not exhibit similar interactions due to the absence of halogen groups. Compound 28 oxygen atom of the methylsulfonyl substituent, as well as the compound 16 oxygen atom of the ethyl ester, which acted as a hydrogen-bond acceptor, formed a hydrogen bond with LYS 211. However, compound 7 did not exhibit the same interaction as a hydrogen-bond acceptor was absent in the R1 position, which was consistent with the hydrogen-bond acceptor contour map of CoMSIA.Compound 22 has cation-πinteractions with ARG367 andπ-πstacking interactions with TYR148 and TYR365, similar to compounds 10 and 1. Notably, the methylsulfonyl group of compound 22 reduced imidazo[1,5-c]pyrimidine core electron cloud density and possessed more interactions with EED than compound 10. The corresponding findings validate the CoMFA model in that the electron-withdrawing group acts to reduce activity in the R1 position.

Fig. 6. Molecular docking interactions of inhibitors with EED. (a) Compound 22 with EED.(b) Compound 10 with EED. (c) Compound 1 with EED.

3. 5 Design of new derivatives

The CoMFA and CoMSIA contour maps (Fig. 7) provide useful information in designing novel EED inhibitors with imidazo[1,5-c]pyrimidine core. Five compounds were designed, which adopted CoMFA and CoMSIA models to predict their pIC50values. The structure and pIC50values of the designed molecules are displayed in Table 4. These newly designed compounds showed predicted pIC50values close to compound 22. Due to more consideration of the CoMISA model, the pIC50predicted by the CoMISA model for the designed molecule is higher than that by the CoMFA model.To further validate the designed molecules, compound D1 was selected with the highest predicted pIC50values for docking, which showed equivalent interactions with compound 22 and a higher docking score than compound 22.

Fig. 7. 3D-QSAR information obtained from CoMFA and CoMISA contour maps

Table 4. Novel Designed Compounds and Their Predicted pIC50 by the 3D-QSAR Model

D2 9.108 9.785 D3 D4 9.229 9.999 9.293 9.783 D5 9.084 9.823

3. 6 ADME prediction

The ADME properties of the drug candidate were found to be closely related to its therapeutic efficacy. To predict their pharmacokinetics and drug-likeness, the newly designed compounds were submitted to SwissADME, and the results are shown in Table 5. The LogP and LogS were used to evaluate the molecule lipophilicity and solubility, respectively,and the values fell within reasonable ranges, demonstrating that the designed molecules possessed good absorbency and solubility. All of the designed compounds had high human gastrointestinal absorption (HIA) and did not have bloodbrain barrier (BBB) permeation. Moreover, all of the designed molecules had the probability to be CYP3A4 inhibitors, which are CYP enzymes that lower the metabolism of drugs. Furthermore, these molecules satisfied Lipinski's rule, which evaluates the drug-likeness of compounds.

Table 5. ADME Prediction for Novel Designed Compounds

4 CONCLUSION

In the present study, 3D-QSAR and docking were adopted in a series of imidazo[1,5-c]pyrimidine derivatives as EED inhibitors. The CoMFA and CoMISA models were shown to possess significant statistical parameters (CoMFA,q2= 0.792,r2= 0.994,r2pred= 0.740. CoMSIA,q2= 0.873,r2= 0.994,r2pred= 0.810) and demonstrated the relationship between the molecular features and activity of these inhibitors. The docking study revealed the mode of interactions between these inhibitors and EED, confirming the CoMFA and CoMSIA results. According to the CoMFA and CoMSIA contour maps, five small molecules were designed, and additional ADME predictions were carried out for the designed molecules. The ADME prediction results demonstrated that the designed molecules had the potential to serve as anticancer drugs. Accordingly, this study may provide theoretical guidance for the rational design of potential EED inhibitors by adopting a novel strategy in discovering anti-cancer agents through the exploration of EED in silico.