Computational Methods and Molecular Modelling for Some Predicted Tautomers of Pyrimidine Compounds

Introduction: Several molecular modelling techniques and quantum chemical methods have been performed to correlate the chemical structures of the compounds with their physical molecular properties. Methods: Theoretical parameters were used to characterize the molecular structure of the investigated ligands and to study their equilibria mechanism. Molecular modelling data such as the bond length, bond order, bond angles and dihedral angles values were estimated for some pyrimidine compounds, where the data suggested the presence of tautomerism and dynamic equilibria were deduced between all the detected tautomers in the solid state. Results: Four main tautomers were predicted for barbital and thiobarbituric acid, while three tautomers were for thiouracil. Quantum chemical parameters such as the highest occupied molecular orbital energy (EHOMO), the lowest unoccupied molecular orbital energy (ELUMO), energy gap (ΔE), dipole moment (μ), sum of the total negative charge (STNC), electronegativity (χ), chemical potential (Pi), global hardness (η) and softness (σ) were calculated. Conclusion: These studies display a good correlation between the theoretical and experimental data. A R T I C L E H I S T O R Y Received: July 29, 2020 Revised: February 11, 2021 Accepted: March 01, 2021 DOI: 10.2174/2666001601666210402115815


INTRODUCTION
The molecular modelling calculations are widely increasing nowadays for the expectation of the mechanism of the reactions and identification of the products [1]. This saves time and money for understanding the biological phenomena and exploiting artificial models [2].
The common feature of molecular modelling [3] techniques is the atomistic level description of the molecular systems; the lowest level of information is individual atoms (or a small group of atoms), in contrast to quantum chemistry, where electrons are considered explicitly. The benefit of molecular modelling is that it reduces the complexity of the system, allowing many more particles (atoms) to be considered during simulations. Molecular modelling methods [4] are used to investigate the structure, dynamics and thermodynamics of inorganic, biological, *Address correspondence to this author at the Chemistry Department, Faculty of Science, Alexandria University, Alexandria, Egypt; E-mail: drmsmasoud@yahoo.com polymeric, drug discovery, medicinal chemistry and pharmaceutical research. The latter has been formalized as a field of study known as Computer Assisted Drug Design (CADD) or Computer Assisted Molecular Design (CAMD) [5].
The stability of DNA and RNA structure is due to the Hbond base pairing and the base stacking, which is actually an interaction between π-orbital of the aromatic rings of the bases and London dispersion forces [6]. The interesting information about their electronic structure and the weakly held molecular complexes can be obtained by quantum chemical treatment [7].
In general, ab initio methods are able to reproduce laboratory measurements for properties such as the heat of formation, ionization potential and molecular geometry [5]. Semiempirical methods utilize approaches that are similar to ab initio methods, but several approximations are introduced to simplify the calculations. These methods include the Huckel approach for aromatic compounds (in which the outer electrons in conjugated systems are treated, but the inner (or core) electrons are ignored) and the Neglect of Differential Overlap formalisms found in the CNDO (Complete Neglect of Differential Overlap) and INDO (Intermediate Neglect of Differential Overlap). Lesions in DNA are caused by electrons with high and low energy resulting in cancer cell formation [8]. Electron correlation is necessary to obtain accurate charge distribution and dipole moments. The valence adiabatic ionization potential (AIP) calculated by the difference between the energies of the appropriate cation radical and neutral species at their respective optimized geometries was reported for pyrimidines [9]. The technique known as energy minimization is used to find positions of zero gradient for all atoms, i.e., a local energy minimum. Lower energy states are more stable and are commonly investigated because of their role in chemical and biological processes. A molecular dynamics simulation computes the behaviour of a system as a function of time.
The phenomena of tautomerization of pyrimidines compounds were examined by Masoud et al. [10][11][12][13][14][15][16] to study the structural chemistry and the chemical equilibria of such compounds in the solid state. Applying software opened the door for the powerful use of theoretical chemistry to study the phenomena of tautomerization research. The possible tautomeric forms: imin-enamine, thiol-thion and keto-enol tautomerism were calculated as gaseous for the compounds under study. The corresponding geometries were optimized without any geometry constraints for full geometry optimizations. In the present manuscript, it is aimed to give spotlight for structural configurations from the view of molecular modelling data, including geometries (bond lengths, bond orders, bond angles, torsion angles); energies (total energy, heat of formation, binding energy, highest occupied molecular orbital energy [E HOMO ], the lowest unoccupied molecular orbital energy [E LUMO ], energy gap [ΔE]); electronic properties (dipole moments, charges, electronegativity, chemical potential, global hardness and softness) for different biologically active pyrimidine compounds.

Compounds
The color, m.p ( 0 C) and sources of investigated ligands are given in Table 1. Their chemical structures and predicted tautomers are presented (Figs. 1,3,5).

Simulations and Computer Programs
The molecules were built with the ChemBio3D Ultra 11.0., where geometrical structures and electronic properties of the free ligands (Barbital, SBA and TU) and their possible tautomers [10][11][12][13][14] were assigned by optimizing their bond lengths, bond orders, bond angles, dihedral angles and close contact using the MM2 calculation to minimizing energy, (Figs. 1-6) and (Tables 2-4). The structural chemistry was investigated by applying the HyperChem.v8.0.x Professional program using PM3 semi-empirical method to speed up the calculations of quantum chemical parameters such as the highest occupied molecular orbital energy (E HOMO ), the lowest unoccupied molecular orbital energy (E LUMO ), energy gap (ΔE), dipole moment (µ). Also, electronegativity (χ), chemical potential (Pi), global hardness (η), softness (σ) and the sum of the total negative charge (STNC) were calculated to study some nucleic acid bases and their predicted tautomers [15] (Fig. 1-6) and (Tables 5,6).

Molecular Modelling Geometries
The molecular modelling calculations [5] are based on neglecting the possibility of hydrogen bonding using different parameters like (Heat of formation, Gradient norm, Dipole, Charges, Solvation in water, Electrostatic potential, Molecular surface, Spin density, Hyperfine coupling constants and Polarizabilites) which are of major importance to control the modelling of the compounds. The optimized geometries of selected pyrimidines, their predicted tautomers and molecular modelling data are given (Figs. 1-6) and (Tables 2-4).
The magnitude of the bond order is associated with the bond length. The bond order concept is used in molecular dynamics and bond order potentials. [17] The bond orders of the same functional groups atoms show more similarities, where the ethyl groups atoms in different tautomers (1, 2, 3, and 4) ( Fig. 1) are approximately similar in their bond orders (Fig. 2). This lies between ~ 0.95-1 (Å) to prove that all ethyl groups atoms are held together via a single bond. The newly born bond orders of enol atoms that appeared in tautomers (2, 3 and 4), (Fig. 2) were approximately similar ~ 0.9 (Å) for O (7) Meanwhile, the bond angles between atoms in barbital for all its tautomers ( Fig. 1) and (Table 2b) pointed that most of them lie in the range: 105.3281 º -126.5950 º referring to distorted sp 2 hybridization due to distortion in electronic configuration. The perfect sp 2 hybridization at ~ 120 º was detected in tautomer (1)  From the dihedral (torsion) angle (τ) [18], a positive value of (τ) is assigned to the clockwise rotation of up to 180º necessary to bring the front atoms into an eclipsed position with the rear atom, while the negative value of (τ) is assigned to anti-clockwise rotation less than 180º to bring the front atoms into an eclipsed position with the rear atom. It is a very useful descriptive tool as it can describe all combinations of tetrahedral, trigonal and diagonal atoms. The most dihedral angles of barbital in all its tautomers ( Fig.  1) lie in the range: 153.0810º--179.8511º, (Table 2c), referring to the distortion in the linearity of sp 3 hybridization.
The approximately zero dihedral angles lie in the range -1.2883 º-29.2101 º in barbital tautomers pointing that some groups of angles are with +ve sign and others are with -ve sign of the same value, i.e., cancel other. The dihedral angles lie in the range -105.3030 º -149.9379 º in barbital tautomers, referring to deviation from sp 2 hybridization due to distorted electronic effects. However, the dihedral angles lie between 30.3667 º --97.9403 º referring to the strong deviation from the perpendicular angle attributed to the distortion effect.  (Table 2d), the tautomer (1) is more steric hindered between its atoms and less stable than other tautomers (2, 3 and 4), (Fig. 1).
The bond length and bond order values for three forms are given in (Table 4a) and (Fig. 6), respectively, where the data suggested the presence of tautomerism. The bond lengths of thione atoms, C(2)-S (8), are approximately similar for tautomers (1 and 2) (Table 4a) Table 4a). This is valid for tautomers (2 and 3) with the same bond length in both giving thione-enol and thiol-enol forms, respectively ( Table 4a). The formation of enol forms in tautomers (2 and 3) explained the disappearance of N(1)-H(9) bond length in that tautomers. However, this bond length is assigned in tautomer (1) only. All the above observations indicated the phenomenon of tautomerization between thione-keto in tautomer (1), thioneenol in tautomer (2) and thiol-enol in tautomer (3) (Fig. 5).
However, from the view of the bond angles between atoms in 2-TU for all its tautomers, (Fig. 5), it is pointed that most of them lie in the range: 109.2594 º -124.5419 º, (Table 4b), referring to distorted sp 2 hybridization due to distortion in electronic configuration. The perfect sp 2 hybridization at ~ 120 º was detected at C (5) (1) only. This is valid for tautomers (2 and 3) with approximately the same bond angle giving thione-enol and thiol-enol forms in tautomers (2 and 3), respectively. The formation of enol forms in tautomers (2 and 3) explained the disappearance of H(9)-N(1)-C(6) and H(9)-N(1)-C(2) bond angles in that tautomers, while these bond angles take place in tautomer (1) only. All the above observations indicated the phenomenon of tautomerization between thione-keto in tautomer (1), thione-enol in tautomer (2) and thiol-enol in tautomer (3). Similarly, the data are in favour of suggesting dynamic equilibria between all the tautomers of 2-TU in the solid state. [13] Meanwhile, the most dihedral angles of 2-TU in all its predicted tautomers (Fig. 5) lie in the range -170.7076 º --179.9917 º, (Table 4c) (Table 4d), referring to steric hindrance effects between those atoms in all these tautomers (Fig. 5).

Quantum Chemical Parameters of Pyrimidine Tautomers
Quantum chemical methods have been performed to correlate the chemical structures of the compounds with their physical molecular properties [19], using theoretical parameters to characterize the molecular structure of the investigated ligands and to study their tautomerization phenomena, (Figs. 1-6). The calculated quantum chemical parameters such as the highest occupied molecular orbital energy (E HOMO ), the lowest unoccupied molecular orbital energy (E LUMO ), energy gap (ΔE), dipole moment (µ) and those parameters that give valuable information about the reactive behavior such as sum of the total negative charge (STNC), electronegativity (χ), chemical potential (Pi), global hardness (η) and softness (σ) were calculated (Tables 5,6). The concepts of these parameters are related to each other [15] where: The inverse of the global hardness is designated as the softness σ as follows:

Barbital
Four main tautomers for barbital are expected depending on molecular orbital calculations [10] (Fig. 1). Their charge distributions are represented (Fig. 2). Most of the parameters given ( Table 5) are the same for the four tautomers of barbital. The major difference lies in the dipole moment values. The dipole moment, µ, was used to discuss and rationalize the structure [20], where the order here is in the sequence for the tautomers (4) > (3) > (2) > (1). The calculated quantum chemical parameters for the barbital tautomers calculated by the PM3 method are given ( Table 6).

Thiobarbituric Acid
Four tautomers are predicted [10,11,14] for SBA (Fig. 3) and their charge distribution given in Fig. (4). The four energy parameters (Table 5) are the same, where slight variations in the heat of formation values denote the presence of four tautomers. The electronegativity (χ), chemical potential (Pi), global hardness (η) and softness (σ) parameters were calculated by the PM3 method for the SBA tautomers ( Table 6).

2-Thiouracil
In a similar way, three main tautomers are expected for the entitled compound [13,14]: Keto-imino-thione Tautomer (1), Enol-thione-imino Tautomer (2) and Enol-thiol Tautomer (3), (Fig. 5). Their charge distributions are given in Fig. (6). Most of the energy parameters ( Table 5) are more or less the same, with slight variation in the three tautomers of thiouracil. The major difference lies through the heat of formation and the dipole moment values to assign the presence of different tautomers. The calculated quantum chemical parameters for the three tautomers of TU calculated by the PM3 method are given ( Table 6).
According to the Frontier Molecular Orbital theory, FMO, the chemical reactivity is a function of the interaction between HOMO and LUMO levels of the reacting species. [21] The E HOMO indicates the ability of the molecule to donate electrons to an appropriated acceptor empty molecular orbitals, and E LUMO indicates its ability to accept electrons. The lower the value of E LUMO , the more ability of the molecule is to accept electrons. [22] While the higher the E HOMO value of the ligand, the easier is its offering electrons to the unoccupied d-orbital of the transition metals center and the greater is its ability to complexation. Furthermore, the HOMO level is mostly localized on the pyrimidine moiety, especially with imino groups indicating that the preferred sites for an electrophilic attack at the metal center are through the nitrogen atoms. So, the pyrimidine compounds with high coefficients of HOMO density were oriented toward the metal center leading to easier complexation through the π-electrons of the pyrimidine ligands and the lone pair of nitrogen. The HOMO-LUMO energy gap, ΔE approach, which is an important stability index, was applied to develop theoretical models for explaining the structure and conformation barriers in many molecular systems. The smaller value of ΔE means that the compound has more probable complexation ability. [22] The recorded values (Table 6) pointed that tautomer (2) in SBA, TU and tautomer (4) in barbital have the smallest HOMO-LUMO gap (8.317487, 7.350207, 7.541028 and 9.209023 eV obtained by PM3 semi-empirical method) compared with the other tautomers. Accordingly, it could be expected that the above-mentioned tautomers (2 and 4) of molecules have more ability to complexation with the metal ion than the other tautomers of molecules.
Absolute hardness, η and softness, σ are important properties to measure the molecular stability and reactivity. A hard molecule has a large energy gap and a soft molecule has a small energy gap. Soft molecules are more reactive than hard ones because they could easily offer electrons to an acceptor. In the above system, the tautomer (2) in SBA, TU and tautomer (4) in barbital have the highest σ value (0.240457, 0.272101, 0.265215 and 0.217178 eV obtained by PM3 semi-empirical method) be more soft molecule compared with the other tautomers, i.e., more reactive than hard ones because they could easily offer electrons to an acceptor metal and possess high ability for complexation, ( Table 6).

CONCLUSION
Mulliken population analysis is mostly used for the calculation of the charge distribution over the whole skeleton of the molecule [22] to explain that the more negatively charged heteroatom is more able to complexation with the metal center through a donor-acceptor type reaction. Variation in the complexation ability of the pyrimidine compounds depends on the presence of electronegative Oand N-atoms. The calculated Mulliken charges [22] of the atoms are presented (Figs. 2,4,6). The sites of ionic reactivity can be estimated from the net charges on a molecule. Thus, the calculations of the summation total negative charge, STNG (Table 5), were calculated, where the ligand of more active sites has higher STNC and is more reactive toward the metal center.

ETHICS
APPROVAL AND CONSENT TO PARTICIPATE Not applicable.

HUMAN AND ANIMAL RIGHTS
No animals/humans were used for studies that are the basis of this research.

CONSENT FOR PUBLICATION
Not applicable.

AVAILABILITY OF DATA AND MATERIALS
Not applicable.

FUNDING
Declared none.