The crystallographic structure of human DNMT1 (PDB ID: 4WXX) was retrieved from the RCSB Protein Data Bank (PDB) available online: https://www.rcsb.org/ (accessed on 30 June 2023) (Berman et al., 2000). Among the different crystallographic structures of DNMT1 available on PDB we selected PDB ID: 4WXX because it contains a co-crystallized molecule of S-adenosyl-L-homocysteine (SAH), and was diffracted with a resolution of 2.62 Å. SAH is reported to be a potent inhibitor of both DNA and histone transmethylation (Halsted and Medici, 2016), therefore this 3D structure could be of interest as a model for inhibitory interactions (Alkaff et al., 2021). The protein preparation was made with default settings of the QuickPrep module of Molecular Operating Environment (MOE) v. 2022.02 (“Molecular Operating Environment (MOE). Chemical Computing Group Inc.: Montreal, QC, Canada, 2023”): addition of all the lacking hydrogen atoms, protonation state at pH 7, elimination of water molecules 4.5 Å farther from the protein and inside the SAH cavity, addition of missing amino acids residues (breaks of up to ten residues and terminal out gaps of up to five residues) and for larger gaps, neutralization of the endpoints adjoining empty residues and energy minimization. The parameters employed for the energy minimization stage were from the AMBER14:EHT forcefield [ff14SB (Maier et al., 2015) for the protein; MAB forcefield (Gerber and Müller, 1995), and AM1-BCC charges for SAH (Jakalian et al., 2002)]. The energy minimization of the protein in MOE is carried out with three successive nonlinear methods: steepest descent, conjugate gradient, and truncated Newton.

The 153 ligands with reported enzymatic activity against DNMT1 in a biochemical assay were obtained from ChEMBL API v. 32 (Davies et al., 2015; Mendez et al., 2019). Only molecules with binding assay type and unequivocally assigned IC₅₀ were selected. Compounds with nucleoside scaffolds (Supplementary Figure S1) were removed using RDKit library (Landrum et al., 2023) substructure search with SMARTS. Before docking (vide infra), the 153 ligands were built and their geometry was energy minimized using MFF94x forcefield implemented on MOE software. For every ligand, the dominant protonation state at physiological pH (7.4) was chosen (“Molecular Operating Environment (MOE). Chemical Computing Group Inc.: Montreal, QC, Canada, 2023”).

The file with the prepared ligands was split with the LeFrag module (Lephar Research, 2023), and Open Babel v.3.1.1 (O’Boyle et al., 2011) was used to convert to .pdb format. Protein and ligands were converted to.pdbqt with MGLTools v.1.5.6. The molecular docking was carried out with Vina v.1.2.3 (Trott and Olson, 2010; Eberhardt et al., 2021) with an exhaustiveness of 8 and 5 binding modes to output. The best score for each ligand was selected for further analysis, with the code freely available at https://github.com/DIFACQUIM/Docking. The grid box was centered in the coordinates: -47.673, 61.885, 6.256 (x, y, z) with a search space of 17 × 25 × 14 Å.

Docking with Ledock (Kirkpatrick et al., 1983) was carried out in the SAH cavity with the default settings of the software: the grid centered 4 Å around the co-crystallized SAH, twenty docking runs for every ligand and 1 Å for the root mean square deviation (RMSD) clustering. For further data analysis, the best score for every ligand was selected with the code available at https://github.com/DIFACQUIM/Docking.

Docking with MOE v. 2022.02 was centered on the SAH cavity and molecular docking was carried out with the default settings: placement (method: triangle matcher, score function: London dG) and refinement (method: rigid receptor, score function: GBVI/WSA dG) (Vilar et al., 2008). Using the “Triangle Matcher” method, the compounds were subjected to 30 search steps and the default values for the other parameters. The clusters with an RMSD <2 Å were visually explored. During the docking, the receptor was considered rigid and the ligands flexible. The conformations with the lowest binding energy were selected for additional analysis.

For this analysis, ligands with IC₅₀ equal to, or lower than 10 μM (pIC₅₀ ≥ 5) were labeled as ‘active’, otherwise they were considered ‘inactive’. Notably, a 10 μM value has been used as a general threshold to define active/inactive molecules in other large-scale studies (Sun et al., 2017; López-López et al., 2022). To develop the current consensus docking protocol, we made the approximation that the enzymatic inhibition assays and the activity values reported in ChEMBL are comparable. The RMSD between the docked and co-crystallized binding conformation of SAH was calculated with Open Babel v. 3.1.1 (O’Boyle et al., 2011). From molecular docking scores and the positive class probabilities, receiver operating characteristic (ROC) curves were generated using KNIME software version 4.6.0 (Gómez-García and Medina-Franco, 2022; Berthold et al., 2009). The results were recorded in a comma-delimited CSV file, which included the scores of each docking program and their ligand efficiency (LE) (vide infra).

Docked ligands were ranked according to their predicted scores in ascending order, compounds with higher rank have more negative values, thus better predicted affinity against DNMT1. pIC₅₀ values were also ranked in descending order since a higher value represents a more potent compound (pIC₅₀ ranking). Additionally, LE was calculated individually for each molecular docking score (obtained by Vina, LeDock, or MOE software) with the equation: L i g a n d   E f f i c i e n c y   L E = D o c k i n g   S c o r e   D S H e a v y   A t o m   C o u n t (1)

In Equation 1, the Heavy Atom Count for each ligand was calculated using RDKit library (Landrum et al., 2023). Correlations and graphs were obtained with SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), and seaborn (Waskom, 2021) libraries using Python programming language version 3.10.12.

Since there is not a single “best” scoring function and docking program, it has been established that combining results from different docking programs increases the likelihood of identifying correct docking poses and improve the performance of docking-based virtual screening (Charifson et al., 1999; Houston and Walkinshaw, 2013; Perez-Castillo et al., 2019; Blanes-Mira et al., 2022). In this study, docking scores and LE values from Vina, LeDock, and MOE were used to calculate seven data fusion metrics: maximum, minimum, arithmetic mean, geometric mean, harmonic mean, median, and Euclidean norm (Bajusz et al., 2019). The data fusion metrics were calculated employing SciPy (Virtanen et al., 2020).

Automated de novo design was carried out with alvaBuilder v.1.0.6 (Mauri and Bertola, 2023). Briefly, alvaBuilder combines structural fragments which are obtained from the training sets chosen by the user. The new sets of molecules constructed from the fragments are scored with a scoring function, also chosen by the user (vide infra). Two different training sets were selected as the source of fragments used as construction blocks. The first dataset was retrieved from ChEMBL 31 (Davies et al., 2015; Mendez et al., 2019) selecting compounds with IC₅₀ against DNMT1 equal to, or lower than 10 μM. The second is the diversity subset (PS6) of 5,000 compounds from Life Chemicals (“Diversity Screening Libraries, 2021”) (accessed in August 2021). Both datasets were curated with the same protocol. Briefly, compounds were standardized, the largest component was retained, and compounds were neutralized and reionized to generate canonical SMILES and remove duplicates, as previously published by our group (Sánchez-Cruz et al., 2019; DIFACQUIM, 2020). A random subset of 285 compounds from Life Chemicals was used to match the number of ‘active’ molecules from ChEMBL after curation. We set the scoring function with ranges of descriptors calculated from the molecules with reported biological activity, using alvaDesc 2.0.10 (Mauri, 2020) (the values used for the scoring function are in the Supplementary Table S1): molecular weight (MW), hydrogen bond donors and acceptors, consensus partition coefficient (logP), aqueous solubility (ESOL), synthetic accessibility (SAscore), topological polar surface area (TPSA). The aggregation method was an arithmetic mean with a population size of 70 and 100 iterations. For each training set, 700 molecules were computed. Finally, 1,398 compounds remained after curation.

De novo compounds were used for similarity searching with the commercial library from the Epigenetics Focused Set of ChemDiv (“ChemDiv, 2023”), with 25,883 compounds. Morgan fingerprints of radius 2 (Morgan2) and 3 (Morgan3) (Rogers and Hahn, 2010), along with MACCS keys (166-bit) fingerprint (Durant et al., 2002) were calculated for all compounds with RDKit (Landrum et al., 2023), and similarity was computed with the Tanimoto coefficient. Molecules from ChemDiv that exhibit one of the following similarity values to at least one compound de novo designed were selected for additional analysis: equal to, or higher than 0.30 for Morgan fingerprints radius 2 or 3; or equal to, or higher than 0.80 for MACCS keys. The selection of these thresholds was based on typical values of intuitive high structure similarity for each fingerprint (Medina-Franco, 2012). Similarity values, along with commercial availability criteria, were used to purchase compounds for further evaluation (vide infra).

Compounds obtained from ChemDiv were experimentally tested at the company Reaction Biology in an enzymatic inhibition methyltransferase assay (“Reaction Biology Corporation, 2023”) using the HotSpot^SM platform. Our research group has reported the methodology and results of this biochemical assay, including the identification of 7-amino alkoxy-quinazolines (Figure 1) (Medina-Franco et al., 2022). Briefly, HotSpot^SM is a low-volume radioisotope-based assay that employs tritium-labeled AdoMet (³H-SAM) as a methyl donor. The test compounds diluted in dimethyl sulfoxide were added using acoustic technology (Echo550, Labcyte, San Jose, CA, United States) into an enzyme/substrate mixture in the nano-liter range. The reactions were started by adding ³H-SAM and incubated at 30 °C. Total final methylations on the substrate (Poly dI-dC) were identified by a filter binding method implemented in Reaction Biology. Data analysis was conducted with GraphPad Prism software available at Reaction Biology (La Jolla, CA, United States) for curve fits. The enzymatic inhibition assays were carried out at 1 μM of SAM. The standard positive control was SAH. The compounds were tested in 10-concentration IC₅₀ (effective concentration to inhibit enzymatic activity by 50%) with a threefold serial dilution starting at 100 μM and 200 μM only for F447-0397. Activity percentage values and dose-response curve are reported as provided by the testing laboratory in Figure 8, Supplementary Table S4, respectively.

The validation of the molecular docking protocol was done with two approaches: RMSD values between the docked and co-crystallized binding conformation of SAH, and ROC curves (as detailed in the Methods Section).

The RMSD values for SAH were lower than 2 Å for all docking programs (Vina: 1.586 Å (second pose); LeDock: 1.291 Å; and MOE: 1.214 Å), Supplementary Figure S2 shows the 3D predicted pose of SAH with each software. The calculated values suggest that the docking protocols are able to identify the experimental 3D conformation of SAH found in the crystallographic structure.

Figure 3 shows the ROC curves for all three docking software using the docking scores and the LE. The ROC curves indicated that MOE’s binding scores led to better identification of true positives, in contrast with Ledock and Vina. However, the calculation of LE (Equation 1) is detrimental to the area under the curve (AUC) (Vina: 0.295; LeDock: 0.250; and MOE: 0.084), this suggests that LE does not contribute to discarding inactive molecules (Figure 3B). These results highlight the relevance of considering the ligand size, herein with the heavy atom count, to evaluate the performance of the docking programs.

To have an insight into the data distribution, correlation plots are shown in Figure 4: Vina (4A), LeDock (4B), and MOE (4C). In each plot, the horizontal axis represents the pIC₅₀ ranking. Docking scores, scores’ ranking, and LE are shown in the vertical axis for each docking software. Spearman correlation (ρ) was computed for each plot. ‘Active’ and ‘inactive’ compounds against DNMT1 are represented in different colors. Supplementary Figure S3 shows the correlation plots for all three docking programs plotting the pIC₅₀ values on the horizontal axis.

Although ρ values for docking scores and scores’ ranking are equal for each program, as expected due to the transformation to rank variables, the distribution of the data is more scattered when plotting the scores’ ranking. Despite the fact the highest correlation observed is low (0.55), the correlation plots in Figure 4 indicate that, overall, considering the ligand size improves the performance of the docking program. This is particularly noticeable in the results obtained with Vina and LeDock where the ρ values improved when considering the LE (Figure 4). The observations obtained with the correlation plots agreed with the conclusions obtained from the ROC curves (Figure 3).

As discussed in the Introduction, consensus SBVS could be more accurate at identifying active compounds as compared to individual methods (Wang and Wang, 2001; Houston and Walkinshaw, 2013). Data fusion also helps to rationalize the relationships between the chemical, physicochemical, and biological features that explain in more detail the possible binding mechanism of different kinds of inhibitors (López-López and Medina-Franco, 2023). Data generated with consensus docking has been useful in developing new drug candidates (Maia et al., 2020; Morris et al., 2022). The advantage of using different docking programs is that the variety of generated conformers is enriched because every program has its own conformer generation algorithm and scoring functions.

Analysis of consensus docking results with data fusion metrics has been shown to improve the results of individual docking (Bajusz et al., 2019; Triches et al., 2022; López-López and Medina-Franco, 2023). Table 1 summarizes the resulting correlations (ρ) of the pIC₅₀ ranking and different data fusion metrics implemented in this work (see Methods Section for details). The resulting correlations (ρ) with pIC₅₀ can be found in the Supplementary Table S2. The best performances were achieved with the median and the minimum rules for the docking scores and LE, respectively. There is a higher correlation with LE, in concordance with the results before the consensus.

Figure 5 shows the correlation between different consensus docking approaches obtained with different data fusion rules (described in the Methods Section) and the bioactivity of DNMT1 inhibitors reported in the literature. The two best correlations are shown as calculated with Spearman’s correlation: the pIC₅₀ values with the median docking score (ρ = 0.372) and with the minimum LE (ρ = −0.564). In agreement with the results discussed in Section 3.1, LE had the best correlations, which further emphasizes the convenience of accounting for the size of the ligand while doing docking analysis with DNMTis.

As observed from Figure 5, correlation values increased from individual docking scores in the case of Vina and LeDock. Minimum LE also shows a better correlation than LeDock and MOE alone. Nevertheless, the calculated correlation for Vina LE has a close value (ρ = 0.559).

The corresponding ROC curves are shown in Figure 6 emphasizing the improved performance of the consensus median and consensus minimum.

As a result of the similarity searching using the 1,398 molecules proposed with de novo design, six compounds from ChemDiv were purchased. The results of the similarity calculations are provided as supplementary .csv files. Figure 7 shows the chemical structures of the newly tested compounds with DNMT1, as well as the most similar de novo compound according to Morgan fingerprint of radius 2 (Morgan2) (Rogers and Hahn, 2010). Most similar de novo compounds according to Morgan fingerprint of radius 3 and MACCS keys are in Supplementary Table S5. The database of the de novo-designed molecules is available at https://github.com/DIFACQUIM/DNMT1-Protocol/tree/main/De-Novo_inspired. Compound F447-0397 had inhibition against DNMT1 in the enzymatic assays, with an IC₅₀ of 41.3 ± 2.1 μM. Dose-response curves used for the calculation of IC₅₀ are in Figure 8. The data used by the testing laboratory (Reaction Biology) to obtain the IC₅₀ is in Supplementary Table S4. It should be noted that, although one point was excluded from the curve fit, F447-0397 inhibited in more than 99% the enzymatic activity of DNMT1 at the highest concentration tested. Nevertheless, the Hill slope of the IC50 curve is not close to 1.0 as it occurs for the positive and internal control, SAH for which the assay conditions to measure DNMT1 were developed by the testing laboratory. Based on these results it is important to conduct additional biochemical and orthogonal assays (e.g., in a cellular context) to further confirm the activity of the compound F447-0397. This compound exhibits a novel scaffold, not previously published among DNMT1 inhibitors to our knowledge. This was shown as no matching molecule was found after the substructure search with the Murcko scaffold of F447-0397 as implemented in RDKit (Landrum et al., 2023), using the curated dataset of 743 molecules with biological activities against DNMT1 found in ChEMBL 33 (Davies et al., 2015; Mendez et al., 2019). Supplementary Table S3 summarizes the results of the enzymatic inhibition assays of the six compounds.

The compounds F447-0397, F447-0509, and F447-0644 are quite similar, in particular, F447-0397 and F447-0509 (an additional methyl group and substitution pattern in the sulfonamide phenyl ring). The docking scores and LE of the three molecules are also similar (Supplementary Table S3), as could be anticipated from their structural similarity. However, the percentage of enzymatic activity at 100 μM is quite different, with F447-0397 being the only inhibitor (12.76%). These are good examples of activity cliffs: compounds with similar chemical structures but very unexpected activity differences (Maggiora, 2006). Although the IC₅₀ of F447-0397 indicated that it is not a very potent compound (41.3 μM - and could be considered “inactive,” the scaffold is novel and could be an interesting starting point for optimization). The novel DNMT1 inhibitor has a “long scaffold” e.g., four-ring systems connected with one-to-three bond linkers. This is in line with other DNMT1 inhibitors with “long or extended scaffolds,” such as the 4-aminoquinoline SGI-1027 and its analogs (Datta et al., 2009; Gros et al., 2015) and glyburide (Juárez-Mercado et al., 2020) (Figure 1). However, unlike SGI-1027 and glyburide, F447-0397 was identified by a combination of de novo design and similarity searching.

Figure 9 shows the predicted binding mode of F447-0397 with DNMT1 generated with Vina, the docking program that had, overall, the best performance of all three docking programs (as shown in Figure 3). The predicted pose shows a hydrogen bond between Glu1168 and the piperazine ring of F447-0397. This could be a key interaction since the co-crystallized SAH also makes a hydrogen bond interaction with Glu1168. Re-docking of SAH, with the three software, predicted the same interaction (Supplementary Figure S4). The predicted binding mode with Vina also exhibits a hydrogen bond between Arg1310 and the oxygens of the sulfonamide from F447-0397. MOE predicted pose also showed the hydrogen bond with the oxygens of the inhibitor’s carboxylic acid (Supplementary Figures S5, S6). Interactions with Arg1310 and computational hits were previously observed in separate docking studies with DNMT1 (Bashir et al., 2023), and also between Arg1310 and EGCG (Figure 1) (Assumpção et al., 2020). Of note, LeDock and MOE predicted interactions between Asn1578 and F447-0397, the interaction with this particular residue could provide selectivity towards DNMT1 versus DNMT3A (Yu et al., 2019). This suggests that F447-0397 could be the starting point of an optimization project toward selective DNMT1 inhibitors.

Compound E760-5661 was an activator (142% enzymatic activity under the assay conditions), followed by the structurally related molecule L162-0591 (132% activity, Supplementary Table S3). Although this is an unexpected result (as we were looking for inhibitors), at least there is agreement that compounds structurally similar have similar (activation) profiles. The activation of DNMT1 also has clinical implications, as DNA hypomethylation has been related to various human diseases, like cancer, and cardiovascular diseases (Wilson et al., 2007; Pogribny and Beland, 2009). Other activators, also identified serendipitously, have been recently published (Rodríguez-Mejía et al., 2022). It would remain to confirm the capabilities of compounds E760-5661 and L162-0591 in a cellular context. To this end, a global human DNA methylation assay could be performed, as recently reported by Rodríguez-Mejía et al. that recently identified two DNMT1 activators (Rodríguez-Mejía et al., 2022). It also remains to explore, at the structural level, the activity cliffs identified in this work. Preliminary structural comparisons of the three compounds (F447-0397, F447-0509, and F447-0644) suggest a quite precise protein-ligand interaction of the active compound - F447-0397 - with DNMT1. The structural analogs could be binding in a different binding region that activates the enzymatic activity of DNMT1 at a certain level, possibly by relaxing allosteric autoinhibition of human DNMT1, as recently proposed for two activators of DNMT1. The mechanism of activation of DNMT1 is out of the scope of this study.