PDNA methyltransferase 1 (DNMT1) has emerged as a potential epigenetic target for development of novel anticancer drugs. In the present investigation, we elaborate systematic virtual screening of large collection of natural products as DNMT1 inhibitors. The virtual screening was executed by implementing several docking programs, such as Glide_XP, Autodock, CDOCKER and LigandFit. Screening of 5120 commercially available FDA-approved drugs (mostly phytochemicals) from ZINC database ranked ZINC38517271, ZINC70455592 and ZINC31983781 as top 3 hits with respective Glide score as high as -14.33, -14.18 and -13.08 kcal/mol. In order to investigate the inhibitors that can endure mutations along the active site residues of DNMT1, mutational analysis at Cys1191 residue, followed by docking and molecular dynamic (MD) simulation were carried out. The inhibitors exhibited minimum binding efficiency on substitution to C1191I while binding affinity remained constant for C1191K mutation. The identified novel ligands were subjected to MD simulation, free energy calculation and energy decomposition per residue. The free energy calculation was based on molecular mechanics- Poisson−Boltzmann surface area (MMPBSA) computation. The free binding energy was high for interaction of the novel inhibitors for both wild type (WT) and C1191K mutant DNMT1, while drastic decrease in binding affinity was seen on non-polar substitution of C1191I. The residues, which contribute to higher binding energies, are Phe1145, Glu1168, Met1169, Cys1191/Lys 1191 and Arg1574. Finally the WT/ mutant DNMT1 -inhibitor complexes were analyzed for interaction with hemimethylated DNA. From the results it is evident that interaction of novel inhibitors with both WT and mutant (C1191I, C1191K) DNMT1 decreases the binding affinity for the hemimethylated DNA.
Keywords: DNA methyltransferases 1; ZINC database; Docking; Mutation; Molecular Dynamics Simulations
DNMT1- DNA Methyltransferase 1;
SAH- S-AdenosylHomocysteine (SAH);
PCNA- Proliferating Cell Nuclear Antigen (PCNA) domain;
NLS- Nuclear Localization Sequence;
MD- Molecular Dynamics;
MMPBSA- Molecular Mechanics-Poisson Boltzmann Surface Area;
MMGBSA- Molecular Mechanics-Generalized Born Surface Area;
WT- Wild Type;
RMSD- Root Mean Square Deviation;
RMSF- Root Mean Square Fluctuation;
PMF- Potential Mean Force;
PLP1- Piecewise Linear Potential;
PDB- Protein Data Bank;
CHARMm- Chemistry at Harvard Molecular Mechanics.
Epigenetic molecular marks characterized by hypermethylation of tumor suppressor genes and hypomethylation of oncogenes, play crucial role in the cancer etiology. These marks are reversible in contrast to genetic mutations, thus they are attractive targets for therapeutic interventions of cancer. Pharmacologic inhibition directly affects the changes in gene expression and can be used in novel cancer therapy. Among the compounds that inhibit epigenetic processes, the most extensively studied are the inhibitors against DNA methylation .
In vertebrates, DNA methylation is associated with the covalent addition of methyl group to 5’-cytosine of CpG dinucleotide clustered in CpG islands at the promoter region of genes. These CpG islands are susceptible to methylation during the normal developmental processes, mainly during embryonic development , stem cell differentiation [3,4], genomic- imprinting , X-chromosome inactivation  and suppression of transposable elements . However, aberrant DNA methylation characterised by hypermethylation of tumor suppressor genes results in gene silencing, which in turn leads to cancer. It has been anticipated that tumor suppressor genes harbor an array of several hundreds of hypermethylated CpG islands mainly on the promoter region . These aberrations have been reported in large number of genes involved in cancer development suggesting that the inhibitors against DNA methylation process has potential for therapeutic intervention of cancer.
Cellular DNA methylation is established and maintained by DNA methyltransferases (DNMTs) which catalyse the transfer of methyl group at 5’cytosine in presence of co-factor SAdenosyl- L-methionine (SAM or AdoMet). There are several isoenzymes of DNMTs known till date, of which DNMT1 play significant role in maintenance of genomic integrity during the process of growth and development. It is responsible for restoring post replicative hemi-methylation to full DNA methylation [9,10]. The functional property of this enzyme is characteristic feature of its structural integrity. The crystal structure of enzyme constitutes 1620 amino acids in complex with S-adenosylhomocysteine (SAH or AdoHcy) extending from N-terminal domain to C-terminal domains . These domain sandwiches 10 characteristic motifs, of which C-terminal domain constituting I, IV, VI, VIII, IX and X motifs are highly conserved . The conserved motifs are responsible for catalytic activity of the enzyme such that motif I and X fold together to form AdoMet binding site while motif IV containing proline-cysteine dipepetide donate thiolate group and Motif IX maintains integrity of target recognition domain. The N-terminal or the regulatory domain constitute proliferating cell nuclear antigen (PCNA) domain, nuclear localization sequence (NLS domain), CXXC-zinc-binding domain and two tandemly arranged bromo-adjacent homology (BAH) domain [12-14]. The regulatory and the catalytic domains act in co-ordination to introduce methylation across the hemi-methylated DNA. The maintenance of DNA methylation brought about by DNMT1 has significant role in adult somatic perturbation cells however perturbation in methylation pattern results in diseases such as imprinting disorders and cancer [15,16]. DNMT1 has emerged as a potential anticancer target. On targeting the catalytic domain of DNMT1 with suitable inhibitors will reduce it activity and subsequently benefit in resorting tumor suppresser gene expression .
A great progress has been made in developing pharmacological inhibitors and some of them have passed through the clinical trials. 5-azacytidine (Vidaza, trade mark), a nucleoside inhibitor arrests DNMT activity by directly incorporating into DNA forming covalent protein-DNA adduct [18,19]. Similarly 5-aza-2-deoxycytidine (decitabine), analogous to 5-azacytidine inhibits DNMT activity and possess antitumor activity in myeloid malignancies [20,21]. However, due to substantial toxicity of these drugs at higher doses there is limitation for their clinical potentials. Beside these nucleoside inhibitors, there is another group of inhibitors i.e. nonnucleoside inhibitors. These inhibitors directly bind to the active site pocket of DNMT1, thus reducing its activity. Some of the non-nucleoside inhibitors known till date are EGCG , curcumin , hydralazine , procaine , procainamide  etc. Although these inhibitors are comparatively less toxic but they have their own limitations as they are less specific in action. Therefore, there is requirement to identify and optimize novel inhibitors which can successfully inhibit DNMT` activity with minimal side effects and higher specificity in action.
It has been reported that mutation of Met 1235→Ala along the active site pocket of mouse DNMT1 (mDNMT1) results in major decrease in catalytic activity for both CpG DNA substrates and hemi-methylated CpG islands. Similarly mutation of Lys1537→Ala influences the substrate recognition potential of mDNMT1. In addition, triple mutation along Arg1237 → Ala, Ser981 → Ala, Tyr983 → Ala, Lys985 → Ala residues exhibit decrease in catalytic activities of the enzyme . Similarly mutation of active site residue Met1235→Ser in human DNMT1 reduced the enzyme activity by 10 fold for 30-mer hemi-methylated substrate. However, the affinity of enzyme for unmodified 30-mer substrate dropped only to 2.2 fold on mutation . Thus DNMT1 enzyme methylates DNA possessively with higher preference and accuracy to hemimethylated target site than unmethylated DNA . As mentioned above, most of the known DNMT inhibitors lack potency, physiochemical and pharmacological properties required to be successful in clinical settings. Our present study aims in identification of novel inhibitors against DNMT1 fulfilling the above challenges. We also intend to investigate how does the mutation along the active site residues influence the catalytic efficiency and binding of the identified novel DNMT1 inhibitors. Thus in quest of novel inhibitors, structure based virtual screening was carried to identify active chemical scaffold leading to the formulation of novel compound. Molecular docking studies play significant role in virtual screening procedure wherein the compounds from large database are screened against the receptor. The compounds from ZINC database library were screened at SAH binding pocket of DNMT1. Herein we employed various docking algorithms to identify compounds of higher efficacy. It was followed by the molecular dynamics (MD) simulation analysis for refinement of final structure and predicting the stability of complexes. The dynamic of the complexes is determined using the time-dependent evolution of the system during the simulation. The non-covalent interaction
constituting hydrogen bonding, van-der Waals and electrostatic occupancies were monitored throughout the docking and simulations. Further ahead MD simulation was merged with molecular mechanics-Poisson, Boltzmann surface area (MM-PBSA, MMGBSA) binding energy (-ΔG) evaluation. The MMPBSA, MMGBSA were carried out to study the dynamic behavior of wild and mutant DNMT1-inhibitor complexes. The change of binding energy (-ΔG) was decomposed on per residue contribution and non-covalent interactions constituting hydrogen bonding, van-der Waals and electrostatic occupancies were monitored throughout the docking and simulations. Several studies done so far based on docking and molecular dynamic simulations suggest that these methods are reliable and have accuracy for determining protein- ligand interactions and calculation of binding free energy [30-32]. Finally the affinity of wild and mutant DNMT1- ligand complexes for substrate (hemi-methylated DNA) was analyzed using Hex-program.
Taken together, the molecular docking and simulations carried for wild and mutant DNMT1-inhibitors complexes constitute multidisciplinary and combinatorial efforts to further develop better inhibitor against DNMT1.
Material and methods
DNMT1 Crystal Structure Retrieval
The X-ray crystallographic structure of wild type (WT) human DNMT1 (PDB id: 3PTA) co-crystallized with SAH at resolutions of 3.6 Å was retrieved from RCSB Protein Data Bank . The hetero-atoms including SAH, zinc ion and water molecules (other than active site) were edited using chimera software. The double stranded DNA attached was clipped off. It followed the addition of polar hydrogen and Gasteiger charges. Further ahead mutant form of DNMT1 for active site residue was constructed using Swiss-Pdb Viewer . In this mutant model, Cys1191 residue was replaced by the remaining 19 amino acids broadly classified into polar, non-polar and electrically charged groups. Both wild and mutant protein structures were optimized using GROMOS 96.1 (43A2) force field within GROMACS 4.5 software package . The protein structure was enclosed in cubic box of 4.0 Å edge lengths. The box was solvated with SPC water model (47,670 molecules) and charged ions, mainly17 chloride ions. The solvated box containing protein and charged ions were subjected to 5000 steps steepest descent minimization 30+ for maximum force of <1000.0 kJ/mol. The optimised wild and mutant DNMT1 structures were subjected to molecular docking analysis.
Structure Based Virtual Screening
Structure based virtual screening has become integral part of drug discovery program, where the molecules are screened against three-dimensional structure of target or reference molecule. Nearly 5120 natural compounds were compiled from ZINC database and screened against the active site pocket of DNMT1. These compounds were prepared using “Prepare ligands” protocol at Discovery Studio 2.5. The preparation of ligand involved removal of duplicate structure, generation of tautomer, isomers, change of ionization state and generating 3D structure.
Molecular docking studies with various algorithms were employed to carry out the screening of these phytochemicals in order to affirm the binding accuracy of the inhibitors. Each algorithm constitutes alternative way of scoring and treating ligand flexibility keeping protein structure in rigid state in order to reduce conformation search space. The various algorithms used for handling ligand flexibility constituted “Lamarckian Genetic Algorithm (LGA) (Autodock)  , monte-carlo conformational search (LigandFit) , descriptor matching (Glide_XP)  and molecular dynamics simulated annealing (CDOCKER) . Results were evaluated on the basis of binding energy (-ΔG) parameter. The more negative the binding, better is the enzyme-inhibitor complex.
AutoDock 3.0.5 is freely available software availed from Scripps Research Institute. Here, the inhibitors were treated as flexible ligands by modifying their rotatable torsions while the protein template was considered to be a rigid receptor. The compounds from the ZINC database were screened to the active site pocket of optimised structure of WT and mutant DNMT1. Grid maps were prepared using auto grid to fix the active site of protein having specific co-ordinates and dimension of 40 × 40 × 40 and a resolution of 0.375 Å. Docking parameters were set as: number of individuals in the population (set to 150), maximum number of energy evaluations (set to 2500000), maximum number of generations (set to 27000), and number of hybrid GA-LS runs (set to 100).
The protein structures (WT and mutant DNMT1) were prepared using the protein preparation module of Schrödinger software. The co-crystallized water molecules were removed. All the selected ligands from ZINC database were assigned an appropriate bond order using the LigPrep 2.4.107 script and converted to .mae format (Maestro, Schrödinger, Inc.) and optimization was carried out by means of the OPLS_2005 force field. The Protein ligand docking studies were performed using Maestro 9.1.107. Parameters having default values were selected and docking was carried out using Glide Extra Precision (XP Glide), version 4.5.19. After the complete preparation of protein and ligand for docking, receptor-grid files were generated. Here van-der Waal radii were scaled of receptor atoms by 1Å with partial atomic charge of 0.25 for running the grid generation module.
CDOCKER is an in-house docking protocol of “Accelyrs Discovery studio”. For the initial stage MD, a soft-core potential was used. Each of the structures from the MD run were then located and fully minimized. The solutions were then clustered according to position and conformation and ranked by energy. CHARMm charges were used for the protein structure, i.e., the param19, toph19 parameter set38 using only polar hydrogen. CDOCKER only allows for flexible ligand treatment. Here the docking model constitutes receptor in its rigid state and static protein conformation of biding site is described using 1.0 Å grid. For every point grid, interaction energy of 20 types of probe atoms was calculated. The three dimensional grid was calculated such that radius of 8Å extend in all directions from any atom in the ligand. Subsequent to simulated annealing conformational search of the flexible ligand, the grid was removed and minimization of all atoms of protein-ligand was performed by fixing the coordinates of the protein using the standard all atom potential function with a distance dependent dielectric (RDIE). This interaction energy was taken as the score for the final ligand pose.
LigandFit is another docking programme of Accelyrs Discovery Studio. It is based on protein minimization using steepest descent (gradient <0.1) and conjugate gradient algorithms (gradient <0.01) of CHARMm force field. The active site determination includes within 10Å radius from the centre of the bound ligand. Docking was performed with monte-carlo simulations using the CFF95 force field. The grid resolution was set to 0.5Å (default), and the ligand accessible grid was defined such that the minimum distance between a grid point and the protein is 2.0Å for hydrogen and 2.5Å for heavy atoms. The grid extends from the defined active site to a distance of 5Å in all directions. The top 10 conformations were saved after rigid body minimizations of 1,000 steps. The scoring was performed using set of scoring functions (including Dock_score, -PMF, -PLP1 and –PLP2) implemented in LigandFit module. The combination of consensus scoring functions was employed to obtain the most preferable output conformation.
Molecular Dynamics (MD) Simulations
The molecular dynamics of the protein–inhibitor complex provides understanding to the flexibility associated with ligand conformational change, and thus provides an insight into molecular basis for inhibition. All of the simulations were carried out using the GROMACS 4.5.5 package with an identical protocol. The best orientation obtained out of docking of protein-ligand complex was used for simulation. We performed simulation for the novel inhibitors identified and their respective binding with human DNMT1. Here we separated the ligands from protein in order to prepare protein and ligand topology file separately. We used GROMOS96 43a1 force filed to generate topology file for protein. The ligand topology file was generated using PRODRG server employing GROMOS96.1 force field . The protein was solvated in a dodecahedron box with edges 1 nm in length using the explicit solvent–simple point charge model (SPC216 water molecules), which generated the water box. The next step followed the 5000 steps of steepest descent minimization and position restrained dynamics to distribute water molecule throughout in 20 ns. The simulation was carried out at 300 k of constant temperature and pressure of 5000 steps for 20 ns using Nose–Hoover method (nvt) and the Parrinello– Rahman method, respectively. Once the system was equilibrated with desired temperature and pressure the final step was to release the position restrains and run production of 500000 steps for 20 ns for data collection.
The binding free energy was calculated using the molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method implemented in Amber12 . For each complex, a total number of 2000 snapshots were taken from an interval of 10 ps in order to obtain low statically error. The complete simulation was carried out for 20 ns. The MM-PBSA method and nmode module of Amber 12 were implemented to calculate the binding free energy of the inhibitor and the detail is summarized as follows:
ΔGb = ΔEMM + ΔGsol − TΔS (1)
Where, ΔGb is the binding free energy in solution consisting of the molecular mechanics free energy (ΔEMM) and the conformational entropy effect to binding (−TΔS) in the gas phase, and the solvation free energy (ΔGsol). ΔEMM was evaluated as
ΔEMM = ΔEvdw+ ΔEele (2)
Where, ΔEvdw and ΔEele stands for van der Waals and electrostatic interactions in gas phase, respectively. The solvation free energy (ΔGsol) was calculated in two steps:
ΔGsol = ΔGpol + ΔGnonpol (3)
Where, ΔGpol and ΔGnonpol are the polar and nonpolar components of the salvation free energy, respectively. The ΔGsol was calculated with the PBSA module of Amber 12 suite. The nonpolar contribution of the solvation free energy was calculated as a function of the solvent-accessible surface area (SAS), as follows:
ΔGnonpol = γ (SAS) + β (4)
Where, the values of empirical constants γ and β were set 0.00542 kcal/ (molÅ 2) and 0.92 kcal/mol, respectively. The contributions of entropy (TΔS) to binding free energy can be evaluated as the sum of change in the translational, rotational, and vibrational degrees of freedom as follows:
ΔS = ΔStranslational + ΔSrotational + ΔSvibrational (5)
TΔS was calculated using classical statistical thermodynamics and normal-mode analysis.
Residue-Inhibitor Interaction Decomposition
The interaction between inhibitors and each residue of DNMT1 was calculated using molecular mechanics of Generalized Born Surface Area (MM-GBSA) decomposition process module of Amber 12 [42,43]. The binding interaction of each inhibitor-residue pair was evaluated in terms of electrostatic (ΔEele) contribution, van der Waals (ΔEvdw) contribution in the gas phase, polar solvation (ΔGpol) and nonpolar solvation (ΔGnopol) contributions.
ΔG inhibitor-residue=ΔEele + ΔEvdw + ΔGpol + ΔGnonpol (6)
The polar contribution (ΔGpol) to solvation energy was calculated by using the GB (Generalized Born) module. The hydrogen bonds (H-bonds) between protein-inhibitor complexes were analyzed using the ptraj module of the Amber 12 program. The linear and the angular distance cutoff between donor acceptor groups were ≤3.5 Å and ≥120° respectively.
Molecular docking of protein-inhibitor complex with DNA
The structure of DNMT1 bound to hemi-methylated DNA was retrieved from protein data bank (PDB). The PDB id of the complex is 4DA4. The complex constituting protein bound to ligand (SAH and zinc ion) and water molecules was clipped off using chimera software. The interaction between DNMT1-inbhitor to hemi-methylated DNA was studied using Hex 6.1 program . Hex is rigid macromolecular program and is based on basic principle of spherical polar Fourier (SPF) program in terms of shape and electrostatics and Van der Waals correlation. Thus in many approaches its working is based on conventional fast Fourier transformation (FFT) docking analysis. The docking analysis is based on Cartesian grid representation of the molecular shape, electrostatic properties and translational FFTs. It is used to perform rigid, flexible and semi-flexible docking. The set default parameter was used to study the interaction between wild and mutant protein-inhibitor complexes and DNA. The study was done to see the affinity of protein-inhibitor complex for substrate i.e. hemi-methylated DNA.
Results and Discussion
Identification of Novel Compounds on Screening of ZINC database
Virtual screening of compounds from chemical database waves path to decoy novel lead compounds. With the accomplishment of human genome project, pharmaceutical targets prediction has increased radically, where virtual screening method play significant role in pharmacogenomics study to identify novel lead compounds against these targets. The recently solved high-resolution X-ray crystallographic structure of DNMT1 (3PTA) presents epigenetic target for treatment of cancer. 5-azacytidine and 5-aza-2’-deoxycytidine approved by the Food and Drug Administration (FDA) of the United States have proved to be boon for treatment of many solid cancers. However, the cellular and clinical toxicity impose checkpoint to use those for further clinical applications. Fortunately, there are several phytochemicals such as tea polyphenol, (-)-epigallocathechin-3-gallate (EGCG) and curcumin which directly binds to the active site pocket of DNMT1, thus suppressing its activity. However, their specificity in action is still in controversy. In the present investigation, therefore, a novel set of compounds of natural origin were screened against DNMT1.
Of total compounds, 5120 FDA approved and commercially available compounds (mostly phytochemicals) were retrieved from ZINC database and were subjected to chemoinformatic analysis to identify molecular scaffolds having large chemical space coverage. These compounds were subjected to “Prepare Ligand” protocol of Discovery studio and 9746 molecules of isomeric and tautomeric forms were extracted for further screening analysis. The molecules were subjected to docking analysis into the SAH binding site of DNMT1. The docking algorithms constitute the sampling of molecules into the coordinate space of binding site and score each possible ligand orientation to filter best possible ligand. Here in our study we used different docking programs, such as CDOCKER, Glide_XP, autodock and LigandFit, each having different algorithm to handle ligand flexibility in the active site pocket of DNMT1. We initiated with the docking analysis of 9746 molecules using CDOCKER protocol of Discovery studio 4.0. The basic protocol was used to dock the molecules into the SAH binding pocket of DNMT1. On screening, 2753 molecules were identified having CDOCKER energy greater than SAH and 306 molecules possessed binding affinity greater than EGCG. To further prioritize the hits, top 306 molecules were subjected to docking analysis using Glide_XP in order to analyse the geometric and energetic complementarity with the catalytic residues of DNMT1. The compounds were ranked using scoring function of Glide-XP and the top 306 hits were visualised using Glide-XP visualizer. The scores obtained ranged from -14.33 to -0.25 kcal/mol for top 306 compounds. There was no correlation between the score and molecular weight rather the molecules were ranked on the basis of better interactions. These molecules scoring high were also subjected to binding energy analysis using Autodock program. The software offers partial flexibility to perform docking with high performance and accuracy. The binding for top 306 ligands were clustered with RMSD cut-off of 2 Å. In parallel, the Monte-Carlo based algorithm LigandFit was also used to study the interaction between DNMT1 and top 306 compounds. The software has docking accuracy and screening capacity through different scoring functions, such as PMF, PLP1, PL2, JAIN, LigScore1 and Lig- Score. The DockScore was evaluated and the binding affinity was calculated as sum total protein-ligand energy and the internal energy.
Finally, the screening based on different algorithms and scoring functions identify top 3 molecules from ZINC database: ZINC38517271, ZINC70455592 and ZINC31983781 to possess higher affinity for DNMT1. The detailed binding energy and their properties have been depicted in Table1.
Table 1. Top 3 Inhibitors Identified on Screening Compounds from ZINC Database and Binding Energy Calculated on Analysis with Glide_Xp, Autodock, CDOCKER and LigandFit.
Robustness of the Novel Compounds against Mutated DNMT1
One of the major challenges in the drug design is to identify inhibitor which can abide the mutations along the active site residues. The basic aim of the present investigation is to comprehend the mechanism behind protein-ligand interactions and to design broad spectrum inhibitor with minimal resistance. The investigation was initiated with the modulation and redesign of the active site pocket of DNMT1 and their subsequent influence on the catalytic efficiency and binding of novel inhibitors. The presence of Cysteine (Cys) residue along the motif IV of active site of DNMTs, which is responsible for nucleophilic attack on carbon-6 of target cytosine forming covalent intermediate, facilitate in addition of methyl group [45,46]. The mutational study of cysteine molecule in several Mtases supports the above mentioned fact. The cysteine residue is conserved along the active site pocket. The present study focuses on how the mutations along the cysteine residues, conserved along the active site pocket will affect the binding of the inhibitors to the SAH binding pocket. The Cys1191 residue was mutated to remaining polar, non-polar, acidic and basic amino acid residues. The efficiency of the novel inhibitors was investigated in terms of mode of interaction and evaluating the change in binding energy. The non-polar substitution of C1191I was detrimental as the binding energy nearly reduced to half while binding affinity nearly remained constant on basic amino acid substitution of C1191K. The detailed binding energy has been depicted in Table 2. The change in the binding energies was evaluated using Autodock 3.0.5. Here, we have evaluated the change in binding energies for top 3 inhibitors with respect to SAH.
This change in binding energy is a consequence of variations in non-covalent interactions, mainly the hydrogen bonding, van der Waals and electrostatic interactions along the active site pocket of DNMT1. The Discovery Studio visualizer was used to explore the detailed interactions. In the active site pocket of wild DNMT1 (3PTA), SAH forms 2 hydrogen bonds while the novel inhibitor ZINC38517271 establishes 10 hydrogen bonds. Similarly ZINC70455592 establishes 8 hydrogen bonds while ZINC31983781 forms 7 hydrogen bonds. The Pi-cation interaction is formed between Arg1312 and ZINC38517271. The detailed interaction of 3PTA with SAH and compounds along with hydrogen bond length is depicted in Figure 1.
Table 2. Binding Energy Calculated for Interaction of Top 3 Inhibitors with WT (3PTA) and Mutant C1191, C1191K DNMT1 using Autodock.
Figure 1. Schematic view of Hydrogen bond propagation on molecular docking analysis of a) WT (3PTA)-SAH, b) WT- ZINC38517271, c) WT- ZINC70455592 and d) WT- ZINC31983781 complex. The Cys1191 residue is labeled in yellow and hydrogen bond in Å is shown by green dotted line.
The mutation along the active site residue of DNMT1 at C1191I results in reduction of binding energy which is the consequence of decrease in non-covalent interactions. While SAH forms only 1 hydrogen bond, 5 hydrogen bonds are established by ZINC38517271 and ZINC70455592 each along the active site pocket. As compared to wild DNMT1, mutant form establishes 3 hydrogen bonds with ZINC38517271. The detailed interactions have been shown in Figure 2. However the binding energy nearly remains constant on mutation at C1191K, as the number of hydrogen bonds formed remains more or less the same. The detailed non covalent interaction of C1191K with the inhibitors has given as supplementary figure1. Thus from the above findings it is clear that small perturbation within the active site architecture influences the catalytic efficiency and dramatically change the affinity of inhibitors towards DNMT1.
Figure 2. Schematic view of Hydrogen bond propagation on molecular docking analysis of a) C1191I-SAH, b) C1191I- ZINC38517271, c) C1191I- ZINC70455592 and d) C1191I- ZINC31983781 complex. The Ile 1191 residue is labeled in green and hydrogen bond in Å is shown by green dotted line.
Finally, DNMT-inhibitor complexes were subjected to molecular dynamics simulations analysis in order to find the stability of the complexes.
Dynamic Characterization of Wild and Mutant DNMT1: Inhibitor Complexes
With an aim to authenticate the docking results and unravel inhibitor efficacy in inhibiting DNMT1, different molecular dynamics simulations were implemented using Gromacs 4.5.5. In order to maintain proper orientation of ligand distance, restraints were applied to inhibitors in the initial few picoseconds (ps) and then whole complexes were allowed to move freely. The docked conformations were analysed by examining their relative protein backbone root mean square deviation (RMSD), total hydrogen bonds, van-der Waals interaction, electrostatic interaction and root mean square fluctuation (RMSF) of active site residues.
In order to explore the effect of mutation on conformational stability of DNMT1 and novel inhibitors, RMSD values were calculated for DNMT1 Cα atoms during 20 ns production phase relative to initial (minimised and equilibrated) and plotted as shown in Figure (3 a-c). The RMSD plot for DNMT1-ZINC38517271 complex attains equilibrium much quickly at 10 ns with RMSD of 2.42 Å. Similarly in case of WT-ZINC70455592 and WT-ZINC31983781complexes; stabilityis achieved quicker than WT-SAH complex taken as reference. On analysis of interaction of these inhibitors with mutant DNMT1; i.e. C1191I, it is clear that that the complexes are unstable as the Cα deviation mounts higher. However, the trajectory of C199IK-DNMT1 inhibitor fluctuates comparatively quicker with the mean of 1.73 Å, 2.00 Å and 2.43 Å. Thus from the above findings it is evident that WTDNMT1 and C1191K-DNMT1 form stable complexes while the non-polar mutation involves decrease in stability of the complexes.
Figure 3. Monitoring the equilibration of the Molecular Dynamic trajectories for each system
a) Time evolution of the root mean squared deviation (RMSD) of backbone atoms for WT- DNMT1 catalytic core domain. b) Time evolution of the RMSD of backbone atoms for mutant C1191I DNMT1. c) Time evolution of the RMSD of heavy atoms for mutant C1191K DNMT1. The four simulations of ZINC38517271, ZINC70455592, ZINC31983781 and SAH bound to the WT, C1191I and C1191K mutant DNMT1 complexes are represented is orange, violet, green and black, respectively. The values reflect the equilibration of each of the systems relative to the initial structures.
The stability of an enzyme-ligand complex is characterised by ligand binding mode inside the active site pocket of the enzyme and the binding force includes strong hydrogen bonds, electrostatic and van der Waals interactions. The intermolecular hydrogen bond plots between WT and mutant DNMT1 and inhibitors are shown in Figure (4 a-c). From the graph, it is apparent that C1191K-inhibitors complexes form more number of hydrogen bond as compared to WT-inhibitor complex. On an average, C1191K-ZINC38517271 complex forms 7.6 hydrogen bonds while WT-ZINC38517271 complex forms 6.5 hydrogen bonds. However, on an average only 2.9 bonds are formed between the C1191I-DNMT1 protein and ZINC38517271 inhibitor. Similarly, the other two inhibitors (ZINC70455592, ZINC31983781) form more stable complexes on interaction with WT and C1191K-DNMT1 protein as compared to non-polar substituted C1191IDNMT1 protein. The protein-ligand complex with respect to SAH has been taken as reference which has forms least number of hydrogen bonds.
The root mean square fluctuation (RMSF) was plotted to compare the variation flexibility of amino acids on interaction with the inhibitors. The RMSF plot in Figure 5 (a-c) depicts that the amino acids along the N-terminal fluctuates higher however the active site residues (1000-1400) fluctuates with an average of 2 Å in case of interaction between WT-DNMT1-inhibitor complexes while the apoform C1191K-inhibitor complex fluctuates to 1.5 Å. However, the non-polar substituted amino acid C1191I fluctuates very high on interaction with an average of 4 Å, which indicated the formation of unstable complex.
Figure 4. Hydrogen bond analysis from Molecular Dynamic simulation study of interaction of inhibitors with (a) WT-DNMT1, (b) C1191I mutant DNMT1, (c) C1191K mutant DNMT1. The inhibitors ZINC38517271, ZINC70455592 and ZINC31983781 and SAH at SAH binding pocket are represented in orange, violet, green and black, respectively.
Figure 5. Comparison of the DNMT1 conformational flexibility.
The root mean squared fluctuations (RMSF) variations of Cα atoms were calculated for (a) WT-DNMT1 amino acid residues, (b) C1191I mutant-DNMT1amino acid residues, and (c) C119K mutant-DNMT1 amino acid residues. The fluctuation in residue is the consequence of interaction with ligands ZINC38517271, ZINC70455592, ZINC31983781 and SAH, represented in orange, violet, green and black, respectively. The lesser the fluctuation, the more stable is the complex.
From the above finding we can infer that WT-DNMT1 and the apoform C1191K-DNMT1 form stable complex on interaction with inhibitors as compared to the C1191I-DNMT1.
From the above molecular dynamics simulation studies constituting RMSD, hydrogen bond analysis and RMSF study, it is evident that the WT and C1191K-DNMT1- inhibitor complexes are stable while non-polar substitution of C1191I is detrimental to the interaction with the inhibitors. Further we performed the binding free energy analysis using MMPBSA method.
Thermodynamic Valuation of Binding Free Energy using MM-PBSA Method
Absolute binding free energies were evaluated using MMPBSA method in order to gain insight into continuous spectrum of binding energy of WT, C1191I and C1191K mutant DNMT1-inhibitor complexes. In this method, the interaction and solvation energy are computed for complex, receptor and ligand in order to investigate average binding free energy. The detailed binding free energies of WT (3PTA), C1191I and C1191K mutant DNMT1 and novel inhibitor complexes have been shown in Figure 6 (a-c). The binding free energy has been evaluated with respect to SAH interaction at the active site pocket of DNMT1. As shown in the Figure 6, WT (3PTA)-ZINC38517271 and C1191K-ZINC38517271 complexes have binding energy of -22.87 and -21.95 kcal/mol, respectively, while the C1191I-ZINC38517271 complex hasan average binding free energy of -8.34 kcal/mol. The binding free energy of ZINC70455592 on interaction with the wild and the apoforms of DNMT1 is -19.48, -18.63 and -7.06 kcal/mol respectively. Similarly the binding free energy for ZINC31983781 has been identified to be -15.05, -14.57 and -6.02 kcal/mol on interaction with respective wild and mutant form of proteins. The higher binding free energy of these inhibitors with respect to SAH in the active site pocket indicates the potential of the novel inhibitors to inhibit DNMT1 enzyme activity. Moreover, the inhibitors are resistant to C1191K mutation as the binding free energy nearly remain constant pertaining to WT protein. However, the lower binding free energy in case of the non-polar substitution of C1191I denotes that the mutation is unfavourable for inhibitors. This decrease in binding energy is a consequence of elevated entropy (−TΔS) parameter.
The differential binding free energy of WT and mutant inhibitor complexes has been carried out to elucidate their interaction mechanism. As shown in Figure 6, it is evident that for all DNMT1-inhbitors complexes, electrostatic and van der Waals interaction in the gaseous phase imparts favourable contribution to the binding of the inhibitors. The nonpolar
Figure 6. Energy components (kcal/mol) for the binding of the top 3 inhibitors at SAH binding pocket of WT (3PTA), C1191I, and C1191K. ΔEele, electrostatic energy in the gas phase; ΔEvdw, van der Waals energy; ΔGnp, nonpolar solvation energy; ΔGpb, polar solvation energy, TΔS, total entropy contribution; ΔGtotal = ΔEele + ΔEvdw + ΔEint + ΔGpb; ΔG = ΔGtotal−TΔS. Error bars in green solid line indicates the difference.
Figure 7. Decomposition of ΔG on a per-residue basis for the WT and mutant C1191I, C119K protein-inhibitor complex: (a) SAH (b) ZINC38517271 (c) ZINC70455592 and d) ZINC31983781. Wild, C1191I and C1191K mutant DNMT1 have been represented as cyan, pink and black, respectively.
solvation energies (ΔGnp) which is due to the embedded inhibitors in the solvent accessible surface area also favours binding affinity in active site pocket of enzyme. However, the polar solvation energies (ΔGpb) and the entropy (−TΔS) are detrimental to the binding free energy. The relatively lower value of non-polar solvation energy indicates the close packing of cavity regions.
Spectrum of Residue Interaction Contributions to binding free energy
In order to obtain the detailed thermodynamic description of residue contribution to the binding free energy, the interaction energies were further decomposed into individual residue contributions through MM-PBSA script in the AMBER 12 suite. This quantitative analysis is significant in understanding the effect of substitution rather mutation and binding mode of the inhibitors in active site pocket. The decomposition of free energy per residue constitutes molecular mechanics and solvation energy. The change in entropy parameter is not included. The differential binding energy spectrum for WT and mutant DNMT1 constituting residueinhibitor complex has been shown in Figure 7 (a-d). This method for residue decomposition aids in understanding the atomistic details of individual residue in protein-inhibitor interactions. From the graph in Figure 7, it is very clear that the inhibitor- residue pair interaction energy is very less for C1191I mutated DNMT1. However in case of WT and C1191K mutated DNMT1, binding energy across the residue is nearly same. On the basis of individual residue contribution to the interaction energy, the common residues which contribute to the binding are Phe1145, Glu1168, Met1169, Cys1191, Lys1191, Arg1574 and Val1576. As seen in Figure 7, the contribution to the binding energy of the above residues varies from -1.47 kcal/mol to -6.01 kcal/mol.
The difference in binding energy of ZINC38517271 bound WT and C1191K system and SAH bound complex is a consequence of individual residue contribution to binding free energy. From the Figure 7b, it is apparent that the high binding energy is due to the contribution of Glu698, Asp700, Cys1191I, Lys1191I, Glu1168, Phe1145 and Trp1170 residues. The binding energy of these residues ranges from -1.91 to 6.01 kcal/mol. The Cys1191 and Lys1191 residues have maximum contribution and the binding free energy for these residues is recorded as high as 6.01 and 6.46 kcal/mol, respectively. However, substitution to C1191I has feeble contribution to the interaction with ZINC38517271 inhibitor.
Similarly, the WT and mutant complex of DNMT1 with ZINC70455592 has elevation in binding affinity in comparison to SAH bound complex. The residues which contribute to the elevated binding energy are Met696, Glu698, Phe1145, Cys1191, Lys1191, Arg1574 and Asn1578. Figure 7c clearly depicts that Cys1191 and Lys1191 have maximum contribution to the binding energy with 5.78 and 6 kcal/mol, respectively. The residues, which contribute to higher binding energy for WT -ZINC31983781 and C1191K-ZINC31983781 complex, are Ala695, Phe1145, Glu1168, Met1169, Cys1191, Lys1191, Arg1570 and Arg1574. Here the Arg1570 has maximum contribution to binding free energy of 3.70 and 3.57 kcal/mol, for wild and mutant DNMT1, respectively.
From the figure 7 it is clearly reflected that amino acid substitution to basic residue C1191K outcomes increase in binding affinity while C1191I mutant substitution is detrimental to inhibitor interaction.
This decomposition of free binding energy (ΔG) per residue is the contribution of van der Waals (ΔEvdw) energy, the sum of electrostatic and the polar solvation energy and the non-polar solvation (ΔGnp) energy.
Contribution of residues to Hydrogen bonding
In order to further shed light on the contribution of residues to hydrogen bonding, cpptraj and ptraj module were performed using AMBER 12. The cpptraj module was used to identify the donor and acceptor groups while ptraj module was used to evaluate the distance, occupancy and angle formed by hydrogen bond. The bonds were counted at distance cutoff of 3.5 Å and angular cutoff of 120°. The details of donor, acceptor groups, distance and occupancy between residues of wild, mutant proteins and inhibitors have been shown in Table 3. The residues, which have maximum contribution to the hydrogen bonding, are Phe1145, Glu1168, Met1169 and Cys1191, Lys1191. These residues play significant role in anchoring the inhibitors with proper orientation into the active site binding pocket of DNMT1. In particular, highly stable hydrogen bond is formed on interaction between OE1 and OE2 atom of Glu and hydrogen atom of ligand. Both in case of WT and mutant C1191K DNMT1, Glu1168 has maximum occupancy of 99.9%. The N atom of inhibitors acts as donor group to Cys1191 and Lys1191 with an average occupancy of 95%. Beside the above mentioned residues, other hydrogen bonds are founds in between OD2 atom of Asp701 and OH group of ligand (ZINC38517271) with bond length of 2.64Å and occupancy of 84%, OD2 atom of Asp1143 and OH atom of ligand (ZINC70455592) having bond length and 82 % occupancy and NH of Met696 act as donor for the ligand (ZINC31983781) with bond length of 3.38Å and has occupancy 78.23%.
The wild and the mutant C1191K form of DNMT1 has an establishment of hydrogen bond with higher occupancy ranging from 99.9%-78.23%. The C1191I mutant DNMT1 offers least space for inhibitors to establish hydrogen bonding. Table 3 clearly reveals that the occupancy of hydrogen bond range from 43.14% to as lower range of 21.23%. It isapparent from the above findings that C1191K substitution of amino acid residue favours the binding of inhibitors in active site pocket of DNMT1 through active participation of hydrogen bonding while non-polar substitution of C1191I is lethal to the interaction of inhibitors to form stable complex.
Table 3. Hydrogen Bond Analysis from the result of Molecular Dynamic Simulation for interaction of top 3 inhibitors with WT and mutant C1191I, C1191K DNMT1.
Interaction of DNMT1-Inhibitor with Hemimethylated DNA.
The catalytic activity of DNMT1 is characterized by its affinity for its substrate, hemimethylated DNA. The N-terminus region of DNMT1 constitutes proliferating cell nuclear antigen (PCNA) domain and aids in recognition of hemimethylated DNA [47-49]. The CXXC domain of DNMT1 shows specificity to unmethylated CpG sites 50. The present study aims to determine the effect of the binding of DNMT1-inhibitor complex to the hemimethylated DNA. Figure 8 shows the structure of hemimethylated DNA retrieved from protein databank. The docking of DNA and DNMT1-inhibitor complex has been carried out using Hex program and binding energy of the complex was evaluated. Figure 9 depicts the interaction of WT and C1191K, C191I mutant DNMT1-inhibitor complex with DNA. Here we have analyzed the affinity WT and mutant DNMT1-SAH complex and DNMT1-ZINC38517271 complex for hemimethylated DNA. The detailed analyses for other complexes have been provided in supplementary figure 2. Figure 9 clearly reveals that the WT-SAH complex has efficient binding affinity for hemimethylated DNA of -345 kcal/mol.The affinity of C1191I, C1191K-SAH complexes for hemimethylated DNA decreases and the binding energies have been identified to -78.21 kcal/ mole and -56.21 kcal/ mol, respectively. However, the affinity of WT and mutant C1191I, C1191K-ZINC38517271 complexes for hemimethylated DNA decreases drastically and the respective binding energies have been identified to be as low as -21.15, -11.43 and -7.23 kcal/ mol. Moreover, the DNA binding site shuffles from the actual binding domain of DNMT1. It is very clear from figure 9 that the hemimethylated DNA interacts near to the catalytic domain of WT-DNMT-inhibitor complex. However, the hemimethylated DNA interacts away from the catalytic domain in C1191I, C1191K-inhbitor complex. From the above result, we can infer that both wild and mutant DNMT1- ZINC38517271 decreases the affinity for hemimethylated DNA.
Figure 8. Structure of hemimethylated DNA of pdb id 4DA4 retrieved from Protein data bank. The methylated cytosine residue in the hemimethylated DNA (labeled circle) is the identifying mark for the DNMT1 to come to the non-methylated cytosine which then protrudes out of the helix, and facilitates its firm interaction with DNMT1 enzyme and gets methylated.
Figure 9. Interaction of hemimethylated DNA with (a) WT-SAH, (b) WT-ZINC38517271, (c) C1191I-SAH, (d) C1191I-ZINC38517271, (e) C1191K-SAH, and (f) C1191K- ZINC38517271 complex. The WT, C1191I and C1191K residues are represented as brown, green and cyan, respectively.
Elaborate systemic virtual screening of 5120 phytochemicals from ZINC database lead to the identification of inhibitors ZINC38517271, ZINC70455592 and ZINC31983781 with highest binding affinity for DNMT1. Comparatively higher binding energy of these inhibitors with respect to SAH indicates their competitive nature at active site pocket of the enzyme. Besides, efficiency of these ligands was also affirmed by their ability to endure mutation, mainly polar substitution (C1191-C1191K) at DNMT binding pocket. Moreover, the molecular dynamic simulation in respect to RMSD, RMSF and average hydrogen bonding substantiated the above finding in correlation with the stability of the novel inhibitors towards the target enzyme. The energetic calculation of free binding energy by MMPBSA suggested that these inhibitors have higher binding energy for both wild and C1191K mutant variety of DNMT1. This high energetic interaction of inhibitors reduces the affinity of DNMT1 enzyme towards the hemimethylated DNA. Thus we summarize in this context that the novel inhibitors we identified, have highly efficient inhibitory effect, as its binding energy is high and appreciable for the target site. Moreover, these inhibitors are resistant to C1191K mutation. These inhibitors being endured to mutation and broad spectrum may accelerate the drug discovery program and open few waves of new path for treatment of cancer.
Granting those A. Shilpi, D. Sengupta, S. Kar, M. Deb and S. K. Rath are thankful to NIT-Rourkela for fellowships under the Institute Research Scheme. S. Parbin and N. Pradhan is thankful to DST, Govt. of India for INSPIRE fellowship. We would also like to acknowledge Dr. Vinod Devraji who helped with application manuals of Schrodinger software.