Identification of a potential DNA methyltransferase (DNMT) inhibitor

Abstract

DNA methyltransferases (DNMTs) play an important role in the epigenetic regulation of gene expression through the methylation of DNA. Since hypermethylation and consequent suppression of tumor suppressor genes are associated with cancer development and progression, DNA hypomethylating agents (HMAs) such as DNMT inhibitors have been proposed for cancer therapy. Two nucleoside analogues approved for the treatment of hematological cancers, decitabine and azacytidine, have poor pharmacokinetic properties, and hence there is a critical need for identifying novel HMAs. Virtual screening of a library of ∼40,000 compounds from the ZINC database, followed by molecular docking of 4,000 compounds having potential druggable properties with DNMT1, DNMT3A and DNMT3B were performed. One unique inhibitor (ZINC167686681) was identified that successfully passed through the Lipinski Rule of 5, geometry constraints as well as ADME/Tox filters and having strong binding energy for DNMTs. Further, molecular dynamics simulations of the docked complexes showed detailed structural features critical for its binding with the DNMTs and the stability of their interaction. Our study identified a compound with potential druggable properties and predicted to bind and inhibit DNMTs. Further investigations involving cellular and animal models of ZINC167686681 will help in potentially taking it into clinical trials for the treatment of cancers.

Communicated by Ramaswamy H. Sarma

KEYWORD:

• DNA methyltransferases
• DNMT inhibitors
• drug discovery
• DNA hypomethylating agents
• virtual screening
• cancer therapy

1. Introduction

DNA methylation is a process of the addition of methyl group at the fifth position of the pyrimidine ring of cytosine residues to generate 5-methylcytosine (5-mC), which commonly takes place at CpG sites in the genome (Kumar et al., 2018). In normal development, it is a heritable epigenetic process that regulates gene function (Moore et al., 2013). Though there are overlaps in the functions, DNA methyltransferase 1 (DNMT1) is an epigenetic regulator enzyme that has a primary role in maintaining methylation (Lyko, 2018). The epigenetic alteration of DNA cytosine residues by methylation plays an important role in controlling gene expression without changing the original DNA sequence (Gao & Das, 2014). During DNA methylation, the MTase domain transfers the adenosyl-l-methionine (SAM) methyl group (–CH3) to the 5′-cytosine position in the DNA sequence, resulting in 5-methylcytosine (5mC) and S-adenosyl-l-homocysteine (SAH) (Moore et al., 2013). DNMT3A and DNMT3B are de novo methyltransferases that are in charge of establishing DNA methylation patterns during embryogenesis and establishing genomic imprints throughout germ cell development (Zhang & Xu, 2017). DNMT3A primarily methylates a subset of genes and sequences in the late stage of embryonic development, particularly after birth, whereas DNMT3B alters a wide region of genomic sequences in early embryos (Zhang & Xu, 2017). More functional details of DNMTs can be found in this review (Lyko, 2018).

DNA methylation is dysregulated in human cancer, resulting in localized hypermethylation and subsequent transcriptional suppression of CpG island-associated gene promoters (Anwar & Lehmann, 2014; Jithesh et al., 2013). DNMT1 is overexpressed in many forms of malignancies, such as breast cancer and acute myeloid leukemia (AML) (Wong, 2021; Wong et al., 2019), which is linked to cellular transformation, whereas lower expression of DNMT1 appears to have a protective effect. Overexpression of DNMT3A and DNMT3B were identified in cancer tissues including triple negative breast cancer (TNBC) and AML (Jinawath et al., 2005; Roll et al., 2008). Furthermore, DNMT3A has been associated with transcriptional activation of DNMTs in many cancers (Jinawath et al., 2005).

Targeting DNMTs is a promising strategy for cancer treatment due to the reversibility of these epigenetic alterations. A few such DNMT inhibitors, either nucleoside analogues or non-nucleoside compounds have been identified, examples of which are shown in Figure 1, and reviewed in detail elsewhere (Lopez et al., 2022). Decitabine (DAC) and azacytidine (AZA) are conventional DNA hypomethylating agents (HMAs) that have demonstrated therapeutic efficacy in the management of hematologic malignancies (Jones & Taylor, 1980). These nucleoside analogues are integrated into DNA where they covalently bind DNMT1, −3A, and –3B, causing dose-limiting toxicities and DNA damage (Agrawal et al., 2018). However, these drugs have poor pharmacokinetic characteristics and stability, which limit target engagement in peripheral tissues (Stresemann & Lyko, 2008). A DNA methylation inhibitor with improved pharmacological characteristics and reduced toxicity is required since these combined qualities may restrict the clinical potential of current DNMT inhibitors. Thus, research into finding better HMAs is currently active, with a recent report of the first-in-class reversible DNMT1-specific HMA (GSK3685032), which is in the preclinical stage (Pappalardi et al., 2021). We hypothesized that screening a large number of compounds and analyzing the interactions with the DNMTs in silico will reveal potential DNMT inhibitors that may be taken forward for further characterization. In this study, we screened 40,000 compounds from the ZINC database (https://zinc.docking.org/) followed by docking with DNMTs and molecular dynamics simulations of selected compounds, revealing the potential catalytic site binding characteristics of some of the compounds.

Figure 1. Examples of DNMT inhibitors. Chemical structure of some of the nucleoside analogues (left) and non-nucleoside analogues (right) are shown here.

2. Methods

2.1. Preparation of the target proteins

The 3D structures of DNMT1, DNMT3A and DNMT3B were retrieved from the PDB (https://www.rcsb.org/), with IDs 4WXX (Zhang et al., 2015), 6W8B (Zhang et al., 2018), and 6U8V (Gao et al., 2020) respectively. Chimera 1.6.1 (https://www.rbvi.ucsf.edu/chimera) was employed for the refinement of the crystal structure of DNMTs through the energy minimization process by applying an amber force field while the geometry optimization was done by using VEGA ZZ (https://www.ddl.unimi.it). Molecular Operating Environment (MOE) was also utilized for the preparation of the structures before docking (Scholz et al., 2015).

2.2. Preparation of the ligand

We obtained the library of 40,000 compounds from the ZINC database (https://zinc.docking.org/) (Irwin & Shoichet, 2005; Sterling & Irwin, 2015). All compounds were downloaded in .sdf format. MONA 2.13 tool was employed to remove redundancy among these compounds. Open Babel was utilized for the conversion of the .sdf format to .pdb (Yoshikawa & Hutchison, 2019). The MOE software was used for optimization of all downloaded compounds by adding partial charges through the protonate3D module. For energy minimization, the MMFF94X force fields were applied for each compound and then added to the MOE ligands database for the purpose of docking (Scholz et al., 2015).

2.3. Virtual screening

In virtual screening, all ligands with poor drug-like properties were eliminated, and the remaining hit molecules were passed through a series of filters for a better selection of drug-like compounds. The Osiris Property Explorer Applet and Osiris Data Warrior were employed to explore the Lipinski filter and geometry constraints of the screened compounds. From this enrichment process, we filtered 4,000 compounds with well-designed geometric constraints, LogP value of 1–5, molecular weight 180–450 kDa, five or less rotatable bonds, <5 hydrogen bonds, and <10 hydrogen bond acceptors (HBAs). Moreover, SCIfinder was used to check whether our filtered hits were previously mentioned in the literature, and then search results were compiled and all those hits with >70% similarity index were evaluated.

2.4. Molecular docking

The catalytic residues of DNMT proteins were identified based on available reports (Horton et al., 2022; Zhang et al., 2015). Site finder tools were also used to find the potential active residues of DNMTs, and then electrostatic surface maps were created around these residues to define docking sites. The MOE tools (Vilar et al., 2008) were used to perform the docking of DNMTs with 4,000 compounds. A triangular algorithm was used to find the 10 best poses of docking molecules (Vilar et al., 2008). We applied London dg function and force fields refinement algorithms to create 10 best poses per molecule and then minimized them. These top 10 molecules were selected based on their binding energies for further interaction analysis with active residues of DNMTs.

2.5. Interaction analysis and ADME/tox properties

LigX tool of MOE was used to prepare the 2D plots of receptor-ligand interactions focused on electrostatic/non-electrostatic interactions, hydrogen bonding, and hydrophobic interactions. For complexes, these interactions assist in validating the binding pattern of the compounds within the active site of the protein (Vilar et al., 2008). For further visualization and interaction analysis, Pymol was used. The 3D structures of protein-ligand complexes were created through Pymol.

After the interaction analysis, the top two best interacting compounds were selected for ADME/Tox (Absorption, Distribution, Metabolism, Excretion, and Toxicity) analysis, which is important to find the possible fate of compounds after administration into the human body. Pk-CSM (Pires et al., 2015) was employed to check the ADME/Tox for the selected compounds and then SwissADME (Daina et al., 2017) and ADMETopt (Yang et al., 2018) were used to crosscheck the pharmacokinetics properties. Finally, the top best compound was selected for MD simulation studies passing the ADME/Tox filter as well as ‘Lipinski’s Rule of Five’ from the molecular docking studies against DNMTs.

2.6. Molecular dynamics simulations

Molecular dynamics (MD) simulations of the holo- and apo DNMT structures were performed using the NAMD 2.13 software (Phillips et al., 2005) and CHARMM36 force field (Huang et al., 2017). Topology and parameter files for the simulations were generated using CHARMM-GUI online server (Jo et al., 2008). The biomolecular simulation systems were generated by dissolving the 3D structural models in explicit solvent using TIP3P cubic water box (Jorgensen et al., 1983) with 10 Å minimum distance between edge of box and any of the atoms in the model. Charges were neutralized by adding 0.15 M NaCl to the solvated system. The biomolecular simulation systems consisted of about 272544, 128529 and 135781 atoms for DNMT1, DNMT3A, and DNMT3B, respectively. Preceding the production run, the biomolecular simulation systems were subjected to energy minimization and thermal equilibration with Periodic Boundary Conditions applied as described previously (Ahmed et al., 2021; Altamash et al., 2021; Geethakumari et al., 2022). Following equilibration, 100 ns production simulation runs were initiated. Timestep of integration was set up at 2 fs and trajectory frames were saved every 10,000 steps. A 12 Å cut-off with 10 Å switching distance was chosen to handle short-range non-bonded interactions, while long-range non-bonded electrostatic interactions were handled using Particle-mesh scheme at 1 Å PME grid spacing (Essmann et al., 1995; Feller et al., 1995; Steinbach & Brooks, 1994).

Analysis of MD trajectories was performed using the available tools in VMD (Humphrey et al., 1996). Backbone root-mean-square deviation (RMSD) measurements were performed using ‘RMSD-trajectory Tool’ in VMD. Calculating root-mean-square fluctuation (RMSF) was performed for residues’ Cα atoms. Radius of gyration (Rog), solvent accessible surface area (SASA), and centre-of-mass distance between protein and ligand were all performed using VMD. Hydrogen bond (H-bond) analysis was performed at a cut-off distance and A-D-H angle of 3.5 Å and 20°, respectively, using the ‘Hydrogen Bonds’ analysis plugin in VMD. Trajectory movies representing the 100 ns simulation runs were generated using VMD ‘Movie Maker’ extension and compiled using Fiji distribution of ImageJ software (Schindelin et al., 2012).

Ligand-binding free energy was estimated through the linear interaction energy (LIE) method [Aqvist, 1994 #5][Hansson, 1998 #6] using the CaFE 1.0 tool [Liu, 2016 #7] in VMD [Humphrey et al., 1996 #8]. Dynamic cross-correlation (DCC) analysis of Cα atoms of the three DNMTs in the absence or in the presence of the ligand was performed using the Bio3D package in R software (Feller et al., 1995; Steinbach & Brooks, 1994). Each trajectory was aligned using the Cα atoms prior to the analysis. The outcome of DCC analysis was represented as heatmap that range from +1 (high correlation, red) to 0 (no correlation, white) to −1 (high anti-correlation, blue).

3. Results

3.1. Virtual screening and molecular docking studies

The current study was conducted to perform a high-throughput integrative computational analysis to identify the most promising drug candidate against DNMTs. Figure 2 shows the overall workflow employed and the criteria used for selection of compounds and the results at each stage in this study.

Figure 2. Overall workflow showing the criteria used for screening and the outcomes at various stages.

A dataset of 40,000 drug-like compounds were retrieved from the ZINC database (Irwin & Shoichet, 2005; Sterling & Irwin, 2015), and virtual screening was performed. The top 4,000 compounds were shortlisted based on geometry constraints and Lipinski Rule of 5 for further studies. To investigate the interactions of lead compounds with the active residues of DNMTs, these compounds were subjected to molecular docking analysis against the catalytic pocket of DNMTs. This analysis resulted in two best-docked compounds (ZINC167686681 and ZINC201939708) with strong binding energy (higher than −10 k cal/mol) towards the active site pocket of DNMTs (Table 1).

Table 1. Binding energies of the top two shortlisted compounds.

Compounds ID	Binding energies (kcal/mol)
	DNMT1	DNMT3A	DNMT3B
ZINC167686681	−13.7,591	−16.3,026	−14.2,510
ZINC201939708	−11.3,841	−12.8,998	−12.8,384

Of all the compounds screened, ZINC167686681 was found to be on top of the list exhibiting best binding energy and potential inhibitory activity against targets DNMT1, DNMT3A, and DNMT3B (Table 1) and hence described here in detail. Molecular docking revealed that the lead inhibitor was bound with the active site of the DNMTs (Figure 3). Comparing the calculated bond lengths between ZINC167686681 and binding residues from the targeted receptors (DNMT1, DNMT3A, and DNMT3B) binding site, it was identified that ligand binding was facilitated by hydrogen bond formation (with a cut-off distance of less than 3.5 Å). Interaction in the ZINC167686681-DNMT1 complex had the maximum number of the hydrogen bond propensities. Both the carbonyl groups of an oxane-3,4-dione moiety of ZINC167686681 were found to be engaged by three hydrogen bonds with C1226 and Q1227 residues at distances 2.2, 3.4 and 2.3 Å, respectively. Similarly, the hydroxyl groups of the benzene-1,3-diol region were connected by five hydrogen bonds with the key residues N1267, V1268, R1310 and R1312 at a distance measured 2.8, 3.3, 2.6, 2.7 and 2.6 Å, respectively. Furthermore, the interaction of ZINC167686681 bound with DNMT3A involved five hydrogen bonds in which the oxane-3,4-dione moiety was engaged by R891 at close vicinity of carbonyl with a bond length of 2.4 Å each. Interestingly, the key residues T645, L888, and S892 located in the close vicinity of the hydroxyl groups of benzene-1,3-diol region was engaged by three hydrogen bonds at a distance of 2.8, 3.1 and 2.2 Å, respectively. Moreover, the binding interaction mode of ZINC167686681 in complex with DNMT3B was through two hydrogen bonds with L829 and C651 residues (distances: 2.9 and 2.8 Å, respectively) facilitated with hydroxyl of benzene-1,3-diol and oxygen of methoxycyclohexane. Regarding the protein-ligand complex formed with ZINC167686681, the ligand was engaged with partially negatively charged groups (carbonyl and hydroxyl groups) at oxane-3,4-dione and benzene-1,3-diol moieties. Thus, the formation of hydrogen bonds with the lead inhibitor was identified in the molecular docking, with hydrogen bonds formed most strongly with DNMT1 amongst the 3 DNMTs studied here.

Figure 3. Interaction of DNMTs with the selected compound. Docked complexes of (A) DNMT3A, (B) DNMT3B and (C) DNMT1 with the lead compound ZINC167686681 are shown.

3.2. Physicochemical and pharmacokinetic properties

To better understand the clinical use of the two ligands, we employed pkCSM, ADMET-SAR and SwissADME to characterize the physicochemical and pharmacokinetic properties. ZINC167686681 had a favorable pharmacokinetic profile. The molecule has a weight of 300.266 g/mol, LogP close to 2 and molecular refractivity index of 76.24 (Table 2). Moreover, its intestinal absorption is high and is a substrate of P-glycoprotein. Furthermore, it cannot cross the central nervous system (CNS) barrier (Table 3), and no hepato- or Ames- toxicities were observed.

Table 2. Molecular properties of the final screening hit (ZINC167686681).

Molecular properties	Compound
Structure
Smiles	COc1ccc([C@@H]2Oc3cc(O)cc(O)c3C(=O)C2 = O)cc1
Molecular formula	C16H12O6
LogP	1.992
Molecular weight	300.266 g/mol
Rotatable bonds	2
HBA	6
HBD	2
Molar refractivity	76.24
Water solubility	Soluble
Surface Area	125.81
TPSA	93.06 A^2

HBA: hydrogen bond acceptor; HBD: hydrogen bond donor; TPSA: topological polar surface area.

Table 3. Pharmacokinetic/toxicity profiles and drug-likeness of the final screening hit (ZINC167686681).

Absorption		Toxicity
ESOL class solubility	Soluble	Ames toxicity	No
Skin permeability	−2.752	Hepatoxicity	No
Ali class solubility	Soluble	Skin sensitization	No
Silicos-IT class solubility	Soluble	Drug likeness
Caco 2 permeability	0.57	Lipinski violations	0
Intestinal absorption (human)	High	Egan violations	0
Gastrointestinal absorption	High	Brenk alerts	0
P-glycoprotein substrate	Yes	Bioavailability Score	0.55
P-glycoprotein I inhibitor	No	Veber violations	0
P-glycoprotein II inhibitor	No	Ghose violations	0
Distribution		Muegge violations	0
BBB permeant	No	PAINS alerts	0
CNS permeant	No	Drug likeness	Excellent
Metabolism		Medicinal chemistry	Excellent
CYP1A2 inhibitor	No	Synthetic accessibility	3.2
CYP2C19 inhibitor	No	Leadlikeness violations	0
CYP2C9 inhibitor	No
CYP2D6 inhibitor	No
CYP3A4 inhibitor	Yes
Excretion
Renal OCT2 substrate	No

ESOL: estimated solubility; BBB: blood-brain barrier; CNS: central nervous system; PAINS: Pan-Assay INterference compoundS.

The compound with the second-best binding energies (Table 1), ZINC201939708, was found to bind to the active site of DNMTs in our molecular docking studies. Based on the measured bond lengths (with a cut-off distance of less than 3.5 Å) between ZINC201939708 and the active site residues in the receptor (DNMT1, DNMT3A, and DNMT3B) binding sites, binding of the compound was found to be facilitated by hydrogen bond interactions (Figure A1 in Appendix). For DNMT1, we found that residues F1229, K1275, K1242 and Q573 formed hydrogen bonds with the compound and had a docking score of −11.3841 kcal/mol. For DNMT3A, the residues R790, R792, V758 and W893 formed hydrogen bonds with a docking score of −12.8998 kcal/mol (Figure A1 in Appendix). Similarly, DNMT3B showed good binding interactions with ZINC201939708 with a docking score of −12.8384 kcal/mol (Table 1). However, ZINC201939708 had a less favorable pharmacokinetic profile (Table A1 in Appendix). The molecule has a weight of 397.42 g/mol and molecular refractivity index of 97.87. Moreover, its intestinal absorption is very low and is a substrate of P-glycoprotein efflux transporter. Furthermore, it can cross the central nervous system (CNS) barrier, though no hepato- or Ames- toxicities were observed. ZINC201939708 also violated some drug-likeness rules like Ghose, Veber, Muegge, and Egan. Hence, we did not further explore ZINC201939708, but took forward the compound with better overall druggable qualities, ZINC167686681, for further analysis using molecular dynamics simulations to understand the stability of binding interactions with the DNMTs.

3.3. Molecular dynamic simulation analysis

After ADME/Tox analysis, we performed molecular dynamics (MD) simulation to study the stability of the ligand-protein complex formed by ZINC167686681 binding to the three DNMTs (Figure 4(A)). All atom, explicit solvent MD simulations were performed for a total of 100 ns using NAMD for each ligand-DNMT complex (holo). Additionally, we performed the simulations of the three proteins alone (apo) for the purpose of comparison. Root-mean-square deviation (RMSD) analysis of Cα of the proteins revealed more dynamic changes in both the apo and holo forms of DNMT1 compared to DNMT3A and DNMT3B, while DNMT3A showed a higher RMSD value compared to DNMT3B (Figure 4(B)). All structures reached an equilibrium RMSD value before 100 ns, suggesting that 100 ns is sufficient for analyzing the trajectories.

Figure 4. Stable interaction of the inhibitor with DNMTs. (A) Composite image of snapshots representing 100 ns of MD simulation, taken 10 ns apart. (B–E) Graph showing Cα root mean square deviation (RMSD) (B), Cα root mean square fluctuation (RMSF) (C), radius of gyration (Rog) (D), solvent accessible surface area (SASA) (E) analyses of apo- and holo-DNMTs over the 100 ns MD simulations. (F–H) Graph showing ligand SASA (F), protein-ligand center-of-mass distance (G) and number of protein-ligand H-bonds (H) in holo-DNMTs over the 100 ns MD simulations.

We then analyzed the root-mean-square fluctuation (RMSF) of Cα atoms, which revealed a region of high dynamic fluctuation in all DNMTs (aa 1450 − 1500 in DNMT1) with an enhanced fluctuation in DNMT1 in this region of the protein (Figure 4(C)). We also evaluated the overall shape of the proteins in the absence and in the presence of the ligand by measuring the radius of gyration (Rog) and solvent accessible surface area (SASA). This analysis revealed a relatively higher Rog and SASA for DNMT1 compared to DNMT3A and DNMT3B, with the latter two showing similar Rog and protein SASA profiles (Figure 4(D,E)). More importantly, there was no observable difference in RMSD, RMSF, Rog, and protein SASA between the apo and holo structure for the most part of the simulation, suggesting that ligand binding did not cause any discernable change in the protein dynamics. Calculation of ligand SASA in the DNMT-ligand complexes showed that the ligand was most buried in the DNMT3B binding pocket (Figure 4(F)). We note that while there was an initial surge in the ligand SASA values, the value returned to basal and was maintained throughout the simulation time in the case of the DNMT3B complex. Ligand SASA in DNMT1 complex showed a brief surge between ∼35 and ∼55 ns and then returned to basal value. On the other hand, ligand SASA values increased dramatically in DNMT3A at about 87 ns and reached the maximum value of ∼500 Ų at ∼93 ns (Figure 4(F)). Calculation of the center-of-mass distances between the ligand and the proteins showed that the ligand leaves the binding site of DNMT3A but returns to interact with the protein at the surface (Figure 4(G)). Importantly, the distance was maintained throughout the simulation run in the case of DNMT3B complex while the DNMT1 complex showed some increase in the middle of the run. Finally, we analyzed the hydrogen bonding between the protein and the ligand in the three complexes which revealed a generally higher number of H-bonds in the DNMT1 and DNMT3A complexes compared to DNMT3B (Figure 4(H)), which is consistent with the docking results. Importantly, we found a sustained H-bond between the ligand and the Glu1266 side chain in the case of DNMT1 with an occupancy of 60.88% (Table A2 in Appendix). Such a highly sustained H-bond (with the Glu1266 residue) was not observed in the case of either DNMT3A or 3B although residues Thr645 in DNMT3A and Asp582 in DNMT3B showed significantly sustained H-bond with the ligand with occupancies of 21.16 and 21.12% (Table A2 in Appendix). Overall, it appears that the ligand, ZINC167686681, can interact with the DNMTs in structurally distinct ways.

To understand the structural dynamics of DNMTs and any alteration in the presence of ligands, we performed dynamic cross-correlation (DCC) analysis of the MD simulation trajectories. DCC has been used to determine dynamically correlated motions in proteins through the analysis of individual amino acid residue pairs in the protein structure over a simulation trajectory. For this, we first aligned all six trajectories using Cα atoms and determined DCC using the Bio3D package in R software (Grant et al., 2006; Skjærven et al., 2014). Our results indicate higher global dynamic correlated and anticorrelated motions in DNMT1 compared to DNMT3a or DNMT3b (Figure 5). Moreover, residue positions between 1,500 and 1,600 seem to have higher dynamically correlated motions in the holo DNMT1 compared to the apo structure. While more global increase in this type of motion is observed in the holo DNMT3a complex compared to the apo. On the other hand, there seem to be a rather global decrease in the dynamically correlated motions in holo DNMT3b compared to apo structure. Finally, we determined the free energy change of ligand binding to the DNMTs using the LIE method (Figure 6) and found the values to be −15.27 ± 2.48, −13.24 ± 5.45 and −16.52 ± 1.45 kcal/mol for DNMT1, DNMT3a and DNMT3b, respectively. Overall, these data provide additional supportive evidence that the inhibitor interacts differently with the different DNMT subtypes.

Figure 5. Altered structural dynamics of DNMTs upon ligand binding. Plots showing dynamic cross correlation of Cα atoms in the apo (upper panel) and the holo (lower panel) DNMT1 (left panel), DNMT3a (middle panel), and DNMT3b (right panel) obtained from MD simulations.

Figure 6. Free energy change of ligand binding. Graphs showing free energy change (ΔG) over the course of 100 ns MD simulations (with an interval of 10 ns) (left panel) and average (ΔG) (right panel) of ligand binding to DNMT1, DNMT3a and DNMT3b.

4. Discussion

The present work finds a novel inhibitor that can potentially inhibit the aberrant function of DNMT1, DNMT3A, and DNMT3B. Molecular docking and molecular dynamics are commonly utilized during the drug discovery process to investigate ligand-receptor interactions (Sheikh et al., 2020). Here, we used virtual screening, molecular docking as well as ADMET analysis, and identified the best inhibitor for DNMTs from a large number of compounds. Moreover, in order to understand the interaction of the identified inhibitor, we performed all-atom, explicit solvent molecular dynamics (MD) simulations using the docked complexes. Additionally, we also performed simulations of the proteins without the ligand (apo) for the purpose of comparison. Generally, the inhibitor appears to bind to all three proteins with the exception of DNMT3A where it starts to unbind after around 90 ns but moves back towards the binding pocket, although the rebinding may not be as tight as it was in the beginning since the ligand continues to show higher solvent accessible surface area. Additionally, it was interesting to observe only a small increase in the protein RMSD in each case.

Two nucleoside analog DNMT inhibitors have been approved by FDA: azacytidine (5-azacytidine) and decitabine (5-aza-20-deoxycytidine) for the treatment of acute myeloid leukaemia (AML), myelodysplastic syndromes (MDS), and chronic myelomonocytic leukaemia (CMML) (Esteller, 2008). Pappalardi et al., reported a novel non-nucleoside DNMT1 inhibitor (GSK3685032) that showed better results in in vitro and animal studies for leukaemia (Pappalardi et al., 2021). Previously, it was reported that Val1580, Asn1578, Met1169 and Trp893, Arg891, Val665 are key residues in the binding site of DNMT1 and DNMT3A for the selectivity of inhibitors (Xie et al., 2019). However, our result is slightly different from this observation, as we observe only residue Arg891 of DNMT3A bind with our inhibitor and other key residues are different from previous findings.

Absorption, distribution, metabolism and excretion are important for efficacy and drug safety, identification of which is an important step in clinical trials (Reddy & Zhang, 2013). We compared the pharmacokinetic properties of our inhibitor with 31 previously reported, identified, and repurposed inhibitors. We identified that our inhibitor showed excellent properties as compared to all other reported inhibitors, including GSK3685032 and decitabine (Table A3 in Appendix). For example, GSK3685032 and decitabine may cause hepatotoxicity (Daina et al., 2017; Pires et al., 2015). The following factors indicate the potential for drug absorption: compound water solubility, cell permeability in the CaCO₂ cell line, skin permeability, and human intestinal absorption, as well as whether the molecule is a P-glycoprotein substrate or inhibitor (Daina et al., 2017). The results obtained in our study showed that ZINC167686681 is soluble in water with high intestinal absorption and has lower CaCO₂ permeability (0.57) as compared to all other compounds (Table A3 in Appendix) whereas, bazedoxifene, a recently reported novel drug for the treatment of TNBC (Huang et al., 2017) is poorly soluble in water and CaCO₂ permeability is 1.068. The pharmacokinetic properties showed that ZINC167686681 follows the Lipinski’s Rule of 5 (Jo et al. 2007) having molecular weight <500 g/mol, hydrogen bond donor <5, and hydrogen bond acceptors <10. Moreover, drug-induced hepatoxicity is a very critical factor in drug designing as more than 1000 drugs induce liver injury and toxicity (Kuna et al., 2018). Our inhibitor was predicted not to induce hepatoxicity and hence can be considered a safe drug. Based on all these findings, we believe that ZINC167686681 could be used as a potential DNMT inhibitor, however, further in vitro and in vivo assays need to be performed to substantiate this observation. Lack of experimental validation is indeed an important limitation of our study.

5. Conclusions

In conclusion, we discovered a new potential DNMT inhibitor (ZINC167686681) having excellent pharmacokinetic properties that demonstrated key interaction with DNMTs. Further preclinical and clinical studies will reveal the potential of this inhibitor in treatment against cancers.

Author contributions

ZJ, WSA performed experiments, analysed data, and contributed to writing. PVJ and KHB conceptualized, implemented, supervised, analysed, and wrote the manuscript. All authors read and approved the final manuscript.

Abbreviations
DNMT	=	DNA methyltransferase
MD	=	molecular dynamics
RMSD	=	root mean square deviation
RMSF	=	root mean square fluctuation
Rog	=	radius of gyration
SASA	=	solvent-accessible surface area
DCC	=	dynamic cross-correlation
DAC	=	decitabine
AZA	=	azacytidine
ESOL	=	estimated solubility
BBB	=	blood-brain barrier
CNS	=	central nervous system
PAINS	=	Pan-Assay INterference compoundS

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Publicly available datasets and tools were employed for this study.

Additional information

Funding

This work was supported by the College of Health & Life Sciences, Hamad Bin Khalifa University. Some of the computational research work reported in the manuscript were performed using high-performance computer resources and services provided by the Research Computing group in Texas A&M University at Qatar. Research Computing is funded by the Qatar Foundation for Education, Science and Community Development (http://www.qf.org.qa). Open Access funding provided by the Qatar National Library.

Appendix

Figure A1. Interaction analysis of DNMT1, DNMT3A and DNMT3B with ZINC201939708. (A) Presents the docked complex of DNMT1 with ZINC201939708. Left side of the figure shows the docked complex and the domains of DNMT1 (CXXC ; BAH1; BAH2; MTase; RFTS). The interacting residues are shown in the inset picture. (B) Shows the docked complexes of DNMT3A and DNMT3B with ZINC201939708. Superimposed image of the docked complex shows the highly similar structure of DNMT3A and DNMT3B. The inset figure shows the interacting residues of DNMT3A and DNMT3B with ZINC201939708.

Table A1. Physicochemical and pharmacokinetic properties of ZINC201939708.

Physicochemical properties
Smiles	C(=O) N[C@@H]2[C@@H] (O)C(CCCc1ccccc1) C(C(=O) O) O[C@H]2[C@@H] (O)[C@H] (O)CO
Structure
Formula	C19H27NO8
Molecular weight	397.42 g/mol
Num. heavy atoms	28
Num. arom. heavy atoms	6
Fraction Csp3	0.58
Num. rotatable bonds	10
Molar refractivity	97.87
TPSA	156.55 Ų
Lipophilicity
Log Po/w (iLOGP)	1.39
Log Po/w (XLOGP3)	0.07
Log Po/w (WLOGP)	−1.33
Log Po/w (MLOGP)	−0.79
Log Po/w (SILICOS-IT)	−0.04
Consensus Log Po/w	−0.14
Druglikeness
Lipinski	Yes
Ghose	No; 1 violation: WLOGP< −0.4
Veber	No; 1 violation: TPSA > 140
Egan	No; 1 violation: TPSA > 131.6
Muegge	No; 2 violations: TPSA > 150, H-don > 5
Bioavailability score	0.11
Medicinal chemistry
Leadlikeness	No; 2 violations: MW > 350, Rotors > 7
Synthetic accessibility	4.62
Pharmacokinetics
Absorption
Water solubility	Soluble
Caco2 permeability	−0.365
Intestinal absorption	24.614 (very low)
Skin permeability	−2.735
P-glycoprotein substrate	Yes
P-glycoprotein I inhibitor	No
P-glycoprotein II inhibitor	No
Distribution
VDss (human)	−1.252
Fraction unbound (human)	0.48
BBB permeability	−1.135
CNS permeability	−4.187
Metabolism
CYP2D6 substrate	No
CYP3A4 substrate	No
CYP1A2 inhibitior	No
CYP2C19 inhibitior	No
CYP2C9 inhibitior	No
CYP2D6 inhibitior	No
CYP3A4 inhibitior	No
Excretion
Total clearance	1.424
Renal OCT2 substrate	No
Toxicity
AMES toxicity	No
Max. tolerated dose (human)	0.913
hERG I inhibitor	No
hERG II inhibitor	No
Oral rat acute toxicity (LD50)	2.125
Oral rat chronic toxicity (LOAEL)	4.836
Hepatotoxicity	No
Skin sensitisation	No
T.pyriformis toxicity	0.285
Minnow toxicity	4.778

Table A2. Occupancy of hydrogen bonds between the protein and the ligand in the complexes.

	DNMT1			DNMT3a			DNMT3b
#	Donor	Acceptor	Occupancy (%)	Donor	Acceptor	Occupancy (%)	Donor	Acceptor	Occupancy (%)
1	Ligand	GLU1266-Side	60.88	Ligand	THR645-Side	21.16	Ligand	ASP582-Side	21.12
2	Ligand	PRO1224-Main	26.55	TRP893-Main	Ligand	12.18	ARG832-Side	Ligand	13.15
3	Ligand	GLU1168-Side	8.18	ARG891-Side	Ligand	9.78	Ligand	GLU697-Side	4.38
4	ASN1267-Main	Ligand	6.99	Ligand	TRP893-Side	4.19	ASN652-Side	Ligand	3.59
5	GLN1227-Side	Ligand	5.59	Ligand	ARG891-Main	3.99	SER833-Main	Ligand	2.39
6	CYS1226-Main	Ligand	4.39	LEU888-Main	Ligand	3.59	Ligand	THR586-Side	2.19
7	Ligand	GLU1266-Main	4.39	THR645-Side	Ligand	1.20	TRP834-Main	Ligand	1.39
8	CYS1226-Side	Ligand	2.40	SER892-Main	Ligand	1.20	Ligand	SER649-Main	1.00
9	PRO1225-Main	Ligand	2.20	Ligand	ASP641-Side	1.20	Ligand	LEU580-Main	1.00
10	Ligand	ASN1578-Main	1.80	ASN711-Side	Ligand	1.00	Ligand	TRP834-Side	0.80
11	ARG1310-Side	Ligand	1.80	PRO709-Side	Ligand	1.00	SER649-Main	Ligand	0.80
12	Ligand	GLY1223-Main	1.60	Ligand	GLY707-Main	1.00	TRP834-Side	Ligand	0.80
13	ASN1578-Side	Ligand	1.40	GLY642-Main	Ligand	0.80	PRO650-Side	Ligand	0.60
14	ASN1267-Side	Ligand	1.40	SER892-Side	Ligand	0.80	ARG733-Side	Ligand	0.60
15	MET1169-Main	Ligand	0.80	Ligand	GLU664-Side	0.80	Ligand	ARG832-Main	0.60
16	VAL1268-Side	Ligand	0.60	LEU888-Side	Ligand	0.60	THR776-Side	Ligand	0.40
17	VAL1268-Main	Ligand	0.60	SER708-Side	Ligand	0.60	PRO650-Main	Ligand	0.20
18	ARG1574-Side	Ligand	0.40	SER669-Side	Ligand	0.60	CYS651-Main	Ligand	0.20
19	GLY1223-Main	Ligand	0.20	Ligand	PHE640-Main	0.40	Ligand	LEU829-Main	0.20
20	ASN1578-Main	Ligand	0.20	TRP893-Side	Ligand	0.40	ASN698-Main	Ligand	0.20
21				Ligand	SER892-Side	0.40	GLY647-Main	Ligand	0.20
22				Ligand	TYR874-Side	0.40	Ligand	PRO650-Main	0.20
23				ILE643-Side	Ligand	0.20	Ligand	ASN652-Side	0.20
24				Ligand	ARG891-Side	0.20	GLY648-Main	Ligand	0.20
25				PRO709-Main	Ligand	0.20
26				ARG887-Side	Ligand	0.20
27				Ligand	SER708-Side	0.20
28				GLN846-Side	Ligand	0.20

Table A3. Pharmacokinetic/toxicity profiles and drug-likeness of the lead compound compared with several medications including DNMT inhibitors.

Ligands	ZINC167686681	GSK3685032	Decitabine	α -yohimbine	Propranolol	Atenolol	Metoprolol	Bazedoxifene	Flubendazole	Niclosamide	Osthole	Risedronate Sodium	Thioridazine	Fluphenazine	Trifluoperazine	Bicalutamide	Enzalutamide	Abiraterone acetate	Orteronel	Seviteronel	Enobosarm	Zoledronic acid	L-NMMA	Chloroquine	Benztropine mesylate	ATRA	Vitamin D3	Pyrvinium pamoate	Salinomycin	Olaparib	Talazoparib	Sacituzumab govitecan
Water solubility	−2.995	−1.238	−2.04	−3.106	−2.882	−1.708	−1.708	−3.873	−3.194	−3.971	−3.603	−0.855	−5.771	−4.336	−4.837	−5.009	−5.755	−6.132	−2.883	−3.149	−4.584	−2.138	Soluble	−4.249	−3.12	−4.924	−6.92	−2.892	−4.053	−4.014	−3.2	−2.892
ESOL class solubility	Soluble	Soluble	Soluble	Moderately soluble	Soluble	Soluble	Soluble	Poorly soluble	Soluble	Moderately soluble	Soluble	Soluble	Moderately soluble	Moderately soluble	Moderately soluble	Soluble	Moderately soluble	Moderately soluble	Soluble	Moderately soluble	Soluble	Soluble	Soluble	Moderately soluble	Soluble	Moderately soluble	Poorly soluble	Insoluble	Poorly soluble	Soluble	Moderately soluble	Moderately soluble
Silico-IT class solubility	Soluble	Soluble	Soluble	Soluble	Soluble	Soluble	Soluble	Poorly soluble	Moderately soluble	Poorly soluble	Soluble	Soluble	Poorly soluble	Poorly soluble	Poorly soluble	Poorly soluble	Poorly soluble	Poorly soluble	Moderately soluble	Moderately soluble	Moderately soluble	Soluble	Soluble	Poorly soluble	Moderately soluble	Soluble	Poorly soluble	Poorly soluble	Poorly soluble	Poorly soluble	Poorly soluble	Insoluble
Caco2 permeability	0.57	1.767	0.092	0.975	1.629	0.632	0.632	1.068	0.885	−3.971	1.644	0.451	1.013	0.949	1.009	1.361	0.719	1.27	0.993	1.242	1.319	−0.796	−0.661	1.624	0.791	1.61	1.223	−1.016	−0.233	1.08	0.59	−1.299
Intestinal absorption (human)	High	High	Low	High	High	Low	Low	High	High	High	High	Very Low	High	High	High	High	High	High	High	High	High	Very Low	Very Low	High	High	High	High	High	Low	High	High	Very Low
Skin permeability	−2.752	−2.105	−2.748	−2.767	−2.783	−3.022	−3.022	−2.735	−2.735	−3.971	−2.009	−2.739	−2.578	−2.746	1.009	−2.794	−2.846	−3.061	−2.737	−2.735	−2.837	−2.735	−2.735	−2.679	−2.735	−2.703	−2.762	−2.735	−2.735	−2.736	−2.735	−2.735
P-glycoprotein substrate	No	No	No	Yes	Yes	Yes	Yes	Yes	No	Yes	No	Yes	Yes	Yes	Yes	Yes	No	No	Yes	Yes	Yes	No	No	Yes	Yes	No	No	Yes	Yes	Yes	Yes	Yes
BBS permeability	No	Yes	No	Yes	Yes	No	No	No	No	No	Yes	No	Yes	Yes	Yes	No	No	No	Yes	No	No	No	No	Yes	No	Yes	No	No	No	No	No	No
CNS permeability	No	No	No	No	No	No	No	No	No	No	Yes	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No
GI absorption	High	High	Low	High	High	Low	High	High	High	High	High	Low	High	High	High	Low	High	High	High	High	High	Low	Low	High	High	High	Low	Low	Low	High	High	Low
CYP1A2 inhibitor	No	Yes	No	No	Yes	No	No	Yes	Yes	Yes	Yes	No	Yes	Yes	Yes	No	No	No	Yes	Yes	No	No	No	No	No	No	No	No	No	No	Yes	No
CYP2C19 inhibitor	No	No	No	No	No	No	No	Yes	No	Yes	Yes	No	No	No	Yes	Yes	Yes	No	Yes	Yes	Yes	No	No	No	No	Yes	No	No	No	Yes	No	No
CYP2C9 inhibitor	No	No	No	No	Yes	No	No	No	No	Yes	Yes	No	No	No	No	No	Yes	No	No	Yes	No	No	No	No	No	Yes	Yes	No	No	Yes	No	No
CYP2D6 inhibitor	No	No	No	Yes	Yes	No	No	Yes	No	No	No	No	Yes	Yes	Yes	No	No	No	No	No	No	No	No	Yes	No	No	No	No	No	No	No	No
CYP3A4 inhibitor	Yes	No	No	Yes	Yes	No	No	Yes	No	Yes	No	No	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	Yes	No	No	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	Yes
Renal OCT2 substrate	No	No	No	No	Yes	No	No	Yes	No	No	No	No	No	No	No	No	No	No	Yes	No	No	No	No	Yes	No	No	No	No	No	No	No	No
AMES toxicity	No	No	No	No	No	No	No	No	Yes	Yes	No	No	No	No	No	No	No	No	Yes	No	No	Yes	No	Yes	No	No	No	Yes	No	No	No	No
Hepatoxicity	No	Yes	Yes	No	No	Yes	Yes	No	Yes	No	No	Yes	Yes	No	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	No	Yes	Yes	Yes	No	No	No	Yes	Yes	Yes
Skin sensitisation	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	No	Yes	No	No	No	No	No	No
Lipinsk violations	0	0	0	0	0	0	0	0	0	0	0	0	1 violation	0	1 violation	0	0	1 violation	0	0	0	0	0	0	0	1 Violation	1 Violation	3 Violations	2 violations	0	0	3 violations
Ghose violatons	0	0	1 Violation	0	0	0	0	2 Violations	0	0	0	0	0	0	0	0	0	1 Violation	0	1 Violations	0	1 Violation	1 Violation	0	0	1 Violation	2 Violations	4 Violations	4 violations	0	0	3 violations
Veber violations	0	0	0	0	0	0	0	0	0	0	0	1 Violation	0	0	0	0	0	0	0	0	0	1 Violation	1 Violation	0	0	0	0	2 Violations	2 violations	0	0	2 violation
Egan violations	0	0	0	0	0	0	0	1 Violation	0	0	0	1 Violation	0	0	0	0	0	1 Violation	0	0	0	1 Violation	1 Violation	0	0	0	1 Violation	2 Violations	2 violations	0	0	1 violation
Muegge violations	0	1 Violation	0	0	0	0	0	1 Violation	0	1 Violation	0	2 Violation	1 Violation	0	1 Violation	0	0	1 Violation	0	0	0	2 Violations	2 Violations	0	0	1 Violation	2 Violations	3 Violations	4 violations	0	0	7 violation
Bioavailability Score	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.85	0.55	0.56	0.11	0.55	0.55	0.17
Druglikeness	Excellent	Best	Best	Best	Best	Best	Best	Intermediate	Best	Best	Best	Intermediate	Best	Best	Intermediate	Best	Best	Best	Best	Best	Best	Intermediate	Poor	Best	Best	Poor	Poor	Poor	Poor	Best	Best	Poor
Medicinal chemistry	Excellent	Best	Best	Poor	Best	Best	Best	Intermediate	Best	Best	Poor	Intermediate	Poor	Poor	Poor	Best	Poor	Poor	Best	Poor	Best	Best	Poor	Poor	Poor	Poor	Poor	Poor	Poor	Poor	Poor	Poor
Synthetic accessibility	3.2	1.62	3.46	4.21	2.43	2.51	2.51	3.48	2.3	2.01	3.13	3.51	4.06	3.58	3.47	3.53	3.17	5.23	3.01	3.4	3	3.14	2.86	2.76	4.79	4.1	6.02	7.9	9.11	3.15	3.82	10
Leadlikeness violations	0	1 Violation	1 Violation	1 Violation	0	1 Violation	1 Violation	2 Violations	0	1 Violation	2 Violations	0	2 Violations	2 Violations	2 Violations	1 Violation	2 Violations	2 Violations	0	2 Violations	1 Violation	0	1 Violation	2 Violations	1 Violation	1 Violation	2 Violations	3 Violations	3 violations	1 violation	1 violation	2 violations

Table A4. Comparison of binding energies and interaction analysis of DNMT1, DNMT3A and DNMT3B with ZINC201939708 and ZINC167686681.

Properties		Compounds
		ZINC201939708	ZINC167686681
Binding energies	DNMT1	−11.3841	−13.7591
	DNMT3A	−12.8998	−16.3026
	DNMT3B	−12.8384	−14.2510
Interacting residues	DNMT1	K1242, Q573, K127, R1238	C1226, Q1227, N1267, V1268, R1310, R1312
	DNMT3A	C710, R891, G890, R790	R891, T645, L888, S892
	DNMT3B	V606, D582, F581, G648, R832	L829, C651