Design, synthesis, molecular modeling and DNA-binding studies of new barbituric acid derivatives

Cancer disease is developing all over the world mainly in developing countries. We should learn more about DNA–ligand interactions to design new drugs that target biological activities like transcription, replication and translation of particular genes. To understand the mechanism of action and design-specific DNA binders, the evaluation of DNA–ligand interactions is critical. Novel barbituric acid derivatives based on (benzyloxy)benzaldehydes were synthesized and evaluated as DNA-binding agents. Among products, molecular docking studies revealed that 4j and 4m have the best interactions with the ctDNA via the minor groove binding. These results were approved by the quantum mechanics calculations. The interaction profiles of the selected compound (4j and 4m) with DNA were evaluated by UV–Visible titration. UV–Visible titration data confirm this interaction. According to the molecular modeling results, the Structure–Activity relationships for all synthesized barbituric acid derivatives were proposed. It was observed that N,N-dimethyl barbituric acid/4-hydroxybenzaldehyde derivatives have better DNA interactions than barbituric acid/vanillin and barbituric acid/3-hydroxybenzaldehyde derivatives. Supplementary Information The online version contains supplementary material available at 10.1007/s13738-022-02576-x.


Introduction
During the cell cycle, the DNA molecule loses its twisted state and its strands disintegrate. This event makes it possible to copy its genetic information [1].
Based on the difference in the rate of transcription of genetic information in normal and cancer cells, the effect of anticancer drugs on the DNA chain of cancer cells is greater than that of normal cells. This process is attractive to researchers and is the source of many anticancer drugs [2]. In recent years, many studies on the binding of organic compounds to DNA have been performed. Studying drug-DNA interaction mechanisms is crucial for designing the beneficial drugs that specifically target DNA [3]. Generally, drugs bind to the DNA through two types of covalent and noncovalent interactions [4].
Drug-DNA Covalent binding is irreversible, leading to complete inhibition of the cell division process and eventual cell death. Cyclophosphamide, Carmustine and Busulfan are examples of drugs that covalently bind to DNA. These drugs are often cytotoxic, non-selective, and have many side effects that can even cause new types of cancer in patients [5]. Consequently, drugs that act through non-covalent interactions with DNA are more attractive to researchers.
Non-covalent interactions cause some drugs to be located between adjacent nucleotide bases to form the DNA-ligand complex. Popular interactions that stabilize the DNA-ligand complex include van der Waals forces, hydrogen bonding, hydrophobic attractions, and charge-transfer forces. DNAbinding agents stop replication in the DNA chain of cancer cells and induce cell death [6].
Actinomycin, Doxorubicin, and Daunorubicin are examples of drugs that bind non-covalently to DNA [7].
One of the applications of spectrophotometric methods is to determine the constant binding of the drug to DNA. UV absorption spectroscopy is the simplest and most common method of studying DNA stability and its interactions with small molecules.
The easiest way to determine the interaction between DNA and drug is the comparison of maximum wavelength absorption (λ max ) of free DNA with the λ max of drug-bound DNA. The observed changes in λ max are considered a parameter to evaluate the strength of the interaction between DNA and the drug [8]. Molecules such as Berenil [9] and Netropsin [10], which have a long linear structure, can become curved and enter the minor groove of DNA.
This study aimed to: (a) synthesize new long linear structure derivatives of barbituric acid, (b) in-silico study their interaction with DNA through docking and simulation of molecular dynamics, and finally (c) evaluate this interaction by UV spectroscopy technique.

Materials and instrumental
Solvents and chemicals were purchased from Merck. ctDNA (Deoxyribonucleic acid sodium salt from calf thymus DNA) were purchased from Sigma.
Melting points are uncorrected and were measured on a Stuart SMP3 apparatus. IR spectra of products were recorded using the ALPHA II Compact FT-IR spectrometer on KBr disks. The NMR ( 1 H and 13 C) spectroscopy was performed using a Varian-INOVA 500 MHz and Bruker Ascend-400 MHz spectrometer. The mass spectra of products were recorded on an Agilent Technology (HP 5975C MSD) mass spectrometer operating at an ionization potential of 70 eV.
The UV-Vis absorption spectra were measured using a PerkinElmer apparatus with 1 cm quartz cuvette.
General procedure for the synthesis of (benzyloxy) benzaldehydes 3 In a 25-mL round-bottom flask, hydroxyl benzaldehyde 1a-c (1 mmol) and benzyl halides 2 (3 mmol) were dissolved in 3 mL DMF. Then 1 mmol potassium carbonate was added and the mixture was stirred at room temperature for 24 h. At the end of the reaction (followed by the TLC, n-hexane/ EtOAc, 4:1), the reaction mixture was poured into the mixture of water-ice. The precipitates were separated, dried and used for the subsequent reaction without further purification.

General procedure for the synthesis of compounds 4a-r
In a 25-mL round-bottom flask, 1 mmol (benzyloxy)benzaldehydes 3 and 1 mmol barbituric acid or N,N-dimethyl barbituric acid were added to 2 mL glacial acetic acid and 2 mL ethanol. Then, the reaction mixture was stirred at room temperature for 24 h. After filtration, the precipitate was washed with water and crystallized in methanol. All products are yellow solid.

Molecular modeling
To optimize the three-dimensional structure of the studied ligands, ab initio calculations were performed using the ORCA quantum chemistry package [12]. The BP86/ Def2-TZVP and MP2/Def2-TZVPP were utilized for optimization and energy calculation, respectively. Force field compatible partial atomic charges were calculated by the AM1-BCC method with the Antechamber package [13]. We used Molekel to visualize the molecular orbital structures [Ugo Varetto, Molekel 5.4]. The Matplotlib was used to draw the UV absorption spectra (https:// matpl otlib. org).

Molecular docking studies
Molecular docking simulations were carried out using AutoDock4 software with default settings to gain a molecular view of the interaction of ligands and DNA duplex. The docked poses for the two of the best-scored hybrids 4j and 4m, shown in Fig. 1, depict its binding along the minor groove of DNA. Crystallographic structures of two DNA strands 1CGC (CG-rich), and 1DNE (AT-rich), were obtained from the RCSB database (http:// www. rcsb. org/ pdb). The molecular docking simulations were performed based on the Lamarckian genetic algorithm of the AutoDock 4.2. The 60 × 60 × 60 Å cubic lattice box was centered on the DNA strand with the lattice distance of 0.375 Å by using AutoGrid [12]. Docking parameters were implemented as the maximum number of energy evaluations 2,500,000, population size 150, gene mutation rate 0.02, crossover rate 0.8, maximum generation 27,000, the genetic algorithm runs 200. The cluster analysis of similar conformations was accomplished with the RMSD (maximum root mean square deviation) of 2.0 Å.

MD simulation studies
Molecular dynamics (MD) simulation was performed by the AMBER99SB-ILDN force field accessible in the GROMACS 4.6.5 software package [14]. The predicted DNA-ligand complex was placed in the center of a dodecahedron box filled with TIP3P water molecules [15]. To neutralize the net charge of the system, enough Mg 2+ ions were added. After each update, the LINCS algorithm was used to reset all covalent bonds to their default lengths [16]. For both long-range electrostatic and van der Waals interactions, the Particle Mesh Ewald (PME) method with a cutoff length of 14 Å and cubic interpolation was used [17].
Energy optimization was performed by 5000 steps of steepest descent and 2500 steps of the conjugated gradient method. After the initial optimization, in the NVT condition (constant number, volume and temperature), while the heavy atoms of ligand and receptor were fixed, the temperature rose from absolute zero (0 K) to 300 K in 300 picoseconds (ps) and remained at this temperature for 200 ps to reach to equilibrium. As the separate temperature-coupling groups, the DNA-ligand and ions-solvent were treated using a modified Berendsen thermostat (V-rescale) [18].
After stabilization of temperature, the MD simulations were performed for 1000 ps in NPT condition (constant number, pressure and temperature), at 1 atm pressure using similar parameters together with isotropic pressure coupling [19]. When the system reached equilibrium, all limitations were removed and the 50 ns MD simulations were accomplished at the NPT ensemble.

Preparation of stock solutions
The stock solution of ctDNA was prepared by dissolving 1 mg of ctDNA in 5 mL of 0.1 M phosphate buffer (K 2 HPO 4 /KH 2 PO 4 ) at pH 7.4 with frequent stirring at room temperature (30 °C) for 24 h. The concentration of ctDNA was calculated by UV absorption spectrophotometry using the molar absorption coefficient (ε 260 = 6600 cm −1 mol 1− , as the molarity of phosphate groups) at 260 nm [20]. The adsorption ratio of ctDNA was measured at 260 and 280 nm to evaluate the protein content rate of the ctDNA solution. The A 260 /A 280 ratio was calculated to be 1.9 (A 260 /A 280 > 1. 8) showed that the ctDNA sample was appropriately free of protein content [21]. The stock solutions of the barbituric acid derivatives were prepared in DMSO with a concentration of 10 mM. Serial dilution of stock solution was carried out using 0.1 M phosphate buffer (pH = 7.4).

Spectroscopy measurements
In docking calculations, the best results were related to the 4j and 4m. Therefore, these products were considered for UV spectrophotometric studies. UV-Visible spectroscopic titrations were done by increasing the concentration of ctDNA (0-200 μl) versus the constant concentration of the ligands. Through titration, an equal amount of ctDNA was added to the barbituric acid derivatives and blank to remove its absorption spectrum. The binding constant (K a ) can be estimated based on the Benesi and Hildebrand equation [22]: In this equation, A 0 indicates the adsorption of barbituric acid derivatives in the absence of ctDNA, A indicates the adsorption of barbituric acid derivatives in the presence of ctDNA, Ka is the binding constant, ε G and ε G-H are the adsorption coefficients of free ligand and its complex with DNA, respectively [23].

Chemistry
The synthetic routes for the preparation of barbituric acid derivatives (4a-r) were depicted in Scheme 1. For this purpose, various (benzyloxy)benzaldehydes (3) were obtained by the reaction of hydroxyl benzaldehydes (1a-c) and benzyl halides (2) in the presence of K 2 CO 3 at room temperature in DMF. The Knoevenagel reaction between (benzyloxy)benzaldehydes (3) with barbituric acid/N,N-dimethyl barbituric acid was carried out using NH 4 HPO 4 as the base in 2 mL glacial acetic acid and 2 mL ethanol at room temperature to obtain the desired products 4a-r in moderate to good yields. The structure of all new barbituric acid derivatives were confirmed by IR, 1 H, 13 C NMR, and mass spectroscopy.
Barbituric acid products derived from 4-hydroxy benzaldehyde or vanillin are shown in Table 1 and the products derived from 3-hydroxy benzaldehyde are shown in Table 2.

In-Silico ADME prediction
To predict the drug-likeness properties of the compounds, the molecular weight factors, the number of hydrogen bonding donors, the number of hydrogen bonding acceptors, the lipophilicity index, and the number of rotating bonds of a compound must be calculated according to Lipinski's rule of five [24].
Physiochemical properties have significant effects on the behavior of compounds within a living system. Hence, it is vital to predict ADME properties through hit identification. We used 2D structures to predict the ADME data of products. Results of estimated parameters indicated that all barbituric acid derivatives (4a-r) met the criteria for druglikeness and followed Lipinski's law (Table 2S).

Quantum mechanics studies
Frontier molecular orbitals are one of the essential parameters for determining the biological activity of the molecules. The HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) can use to determine the stability, reactivity, and biological activity of the molecules. Table 3 shows the results of quantum mechanics calculations for 4j and 4m using the MP2/Def2-TZVPP||PB86/Def2-TZVPMP2 method.
The electrons fill the molecular orbitals from the lowest energy level orbitals up to the HOMO. HOMO energy is equal to ionization potential (I = -E HOMO ). The higher HOMO energy equals the lower ionization potential of a molecule, leading to the increase of the tendency of that molecule to lose an electron and participate in charge-transfer interactions. The LUMO electron density contains the parts of a molecule with a strong tendency to accept the electron. The electron affinity (A = −E LUMO ) is described as the amount of energy released by electron absorption. The greater the electron affinity of a molecule means the greater tendency to capture the electron [21].
The amounts of ionization potentials (I = −E HOMO ) and electron affinity (A = −E LUMO ) in the compounds 4j and 4m were almost equal, so they have a similar tendency to give electrons and participate in charge-transfer interactions.
I and A values of the Pyriproxyfen [25] are calculated and written in Table 2 compared with compounds 4j and 4m. 4j and 4m had lower I and A values than Pyriproxyfen. So, the charge-transfer interaction, as well as the electron affinity of Pyriproxyfen is more than 4j and 4m.

Molecular docking studies
To evaluate the molecular mechanism of DNA interaction in atomic details, molecular docking was performed to predict the best possible binding mode between two DNA strands for the new synthesized barbituric acid derivatives (4a-r). The binding affinities of the 17 molecules of barbituric acid derivatives (4a-r) with two DNA strands were investigated.
We considered the docking of recognized small molecules with two DNA chains, 1CGC and 1DNE. The affinity of small molecules to the CG sequence can evaluate by using 1CGC. 1CGC is a decamer d(CCG GCG CCGG) that contains the repeated CG sequence. Also, to assess the affinity of small molecules to AT sequence, we used 1DNE. 1DNE is a dodecamer d(CGC GAT ATC GCG ) with two repeats of AT at the central section of the strand. Table 3S shows the results of molecular docking of DNA interaction of the barbituric acid derivatives.
The free energy (ΔG) of the ligand-receptor binding indicates the amount of ligand tending to bind to the receptor. The docking process can calculate this parameter [26].
According to the free energy values of DNA-binding (Table 3S), the strongest DNA-bindings are related to the 4j and 4m derivatives. For these compounds, the free energy value for the CG-rich chain is −8.80 kcal/mol and for the AT-rich chain is -8.1 kcal/mol. These values for Pyriproxyfen are −7.40 kcal/mol for the CG-rich chain and -6.7 kcal/ mol for the AT-rich chain. Figure 1 shows that compounds 4j and 4m bind along the minor groove of DNA. Table 4 shows the graphs of DNA-interaction of the 4j and 4m with the CG-rich strand (1CGC) and AT-rich strand (1DNE).
Ligand efficiency (LE) is the free energy binding calculated by docking per heavy atom of the ligand structure. Scheme 1 Synthesis of barbituric acid derivatives (4a-r) Table 1 Barbituric acid products derived from 4-hydroxy benzaldehyde or vanillin It indicates the ability of the functional group (based on molecular weight) to interact with the receptor [27]. This parameter was calculated by the following equation: Table 5 compared the LE values of 4j and 4m with the Pyriproxyfen as a DNA-binding model.
The LE calculated for the 1CGC, for derivatives 4j and 4m, was 0.35 and 0.32, respectively, and 0.3 for the 1DNE chain. Therefore, these two compounds are optimal in terms of this parameter. The LE values for Pyriproxyfen were 0.3 for 1CGC and 0.27 for 1DNE, which were lower than 4j and 4m. These results indicate that in interaction with DNA, the role of atoms in compounds 4j and 4m is greater than that of Pyriproxyfen, which suggests the effectiveness of the design of this group of compounds. Considering the docking error in estimating the ΔG of binding (± 2.0 kcal mol −1 ), it is not reliable to predict the binding affinity for DNA strands only by docking results. Hence, in the next phase, we used docking results as the starting data to perform molecular dynamics simulation (MD) to get more plausible results.

Molecular dynamics simulation
To achieve a more realistic model of ligand-DNA interaction, 50 ns MD simulations were accomplished on the predicted docking model (for 4j and 4m derivatives). To analyze the system's stability during 50 ns of MD simulation, the average total energy and temperature of systems were calculated. The standard deviation (SD) and relative standard deviation (RSD) were calculated and used to evaluate the stability of the systems. The coefficient of variation was calculated by the RSD equation [22]:  In all systems, the rate of fluctuation (RSD) of energy was less than 0.40% and the temperature was less than 1.00%.
The low value of RSD shows that the systems were stable and the principle of total energy conservation was achieved in these systems.
The location of the atoms in each conformer was studied to evaluate the stability of the complexes during the 50 ns MD simulation. In this regard, the RMSD matrixes of 4j or 4m and DNA strands were calculated to display the changes in the conformation of each structure (Tables 11S and 12S).
The RMSD matrix shows conformational changes at any moment compared to all conformers that have been seen through the MD simulation. The observed subsquares in the diameter of the RMSD matrix square show the conformational change. Hence, inside each square, the conformation of the complex was stable. The blue to red color spectrum indicates changes in the RMSD matrix from 0 to more than 0.5 nm, respectively. Tables 11S and 12S show the RMSD matrix of the DNA strands in 1CGC and 1DNE in complex with 4j and 4 m. In the 4j complexes with 1CGC and 1DNE, the changes in the RMSD matrix of the 1DNE chain were more significant than those of the 1CGC strand. 4j showed a stable binding to the DNA strand during the MD simulation. Although the energy released to bind the 4m compound to the 1CGC strand is a smaller amount than 1DNE, given its stability in binding to both DNA strands, its selectivity cannot be stated [28]. The results of the molecular dynamics simulation data are shown in Fig. 2. A comparison of the docking data (Fig. 1) and the MDD data (Fig. 2) shows that although the system has undergone many changes during the simulation, the designed ligands are well placed inside the minor groove. Examination of the behavior of the 4j and 4m showed that they tend to be located between nucleotide bases, although this process was not observed in the simulation time.
To understand the dynamic stability of DNA-ligand complexes, the radius of gyration (R g ) values are calculated [29]. As seen in Table 6, the changes of R g were not statistically significant, and the DNA-ligand complexes have good stability MD.

UV-Visible spectral studies
To confirm the results achieved from MD simulation, the interactions of the 4j and 4m with DNA strands were investigated using the spectrophotometric method.  The UV spectra of titrated compounds with different concentrations of ctDNA are indicated in Fig. 3.
The binding constants of 4j and 4m (K a ) to the DNA strand were measured, and its values are listed in Table 7.
To in vitro investigate the interaction of the 4j and 4m with ctDNA, we measured the protein content of ctDNA to increase the accuracy of the titration analysis. For this purpose, the ctDNA adsorption ratio was measured at 260 and 280 nm. The A 260 /A 280 ratio was calculated to be 1.9. It indicates the absence of protein in the ctDNA sample. Therefore, the obtained K a s indicate the interaction between 4j and 4m and ctDNA.
The absorption spectra of 4j and 4m had a peak at 263 nm and 259 nm, respectively (Fig. 2). The constant values of binding to the DNA strand (K a ) were determined by drawing A 0 /A-A 0 versus 1/[ctDNA] and calculating the intercept ratio to the line slope.
Based on the results, the K a of the 4m was higher than 4j. Therefore, it performs stronger interactions with ctDNA than 4j.
Changes in the UV spectrum of compound 4m during titration had a redshift (bathochromic) displacement. In the spectrophotometric titration study of a ligand with ctDNA, redshift indicates intercalative interactions [8]. So, this compound tends to locate between the nucleotide bases of the ctDNA.

The structure activity relationships (SARs) survey
It can be concluded from the molecular docking data of 17 synthesized products that: Compounds derived from 4-(benzyloxy)benzaldehyde perform a better interaction with DNA than 4-(benzyloxy)-3-methoxybenzaldehyde and 3-(benzyloxy)benzaldehyde. Also, products containing N,N-dimethyl barbituric acid performed better than barbituric acid.

Conclusion
In conclusion, efficient access to the synthesis of novel barbituric acid derivatives with promising DNA interacting properties has been accomplished. Molecular docking studies revealed that among all products, 4j and 4m are the best interacting agents and are capable of interacting with the ctDNA through the minor groove binding. Based on the flexible shape, they match well with the topology of DNA strands. Quantum mechanics calculations predicted similar trends. The interactions of 4j and 4m with DNA were studied by UV-Vis spectroscopy. The data obtained from spectrophotometric measurements confirm this interaction. Based on the molecular modeling results, the Structure Activity Relationships (SARs) for these new barbituric acid derivatives were proposed and it was observed that the derivative of N,N-dimethyl barbituric acid/4-hydroxybenzaldehyde have better DNA-interactions than barbituric acid/vanillin and barbituric acid/3-hydroxybenzaldehyde derivatives.