The SMILES approach was employed to evaluate the absorption, distribution, metabolism, and excretion (ADME) and Lipinski’s rule for five out of the seven compounds. To determine the properties, the pkCSM website (http://biosig.unimelb.edu.au/pkcsm/) was utilized.

The Protox II website (http://tox.charite.de/protox_II) was utilized to predict the toxicity of compounds. The toxicity chemical class prediction was categorized into the following six classes: Class 1: extremely lethal (LD50 ≤ 5), Class 2: fatal (5 < LD50 ≤ 50), Class 3: toxic (50 < LD50 ≤ 300), Class 4: harmful (300 < LD50 ≤ 2000), Class 5: possibly hazardous (2000 < LD50 ≤ 5000), and Class 6: non-toxic (LD50 > 5000) (OECD 2001).

QSPR strives to establish the quantitative correlation between molecular structure and physicochemical properties (Shoombuatong et al. 2017). Physicochemical parameters, including lipophilic, electronic, and steric parameters, were predicted using the Gaussian 09 program suite (Frisch et al. 2009). Lipophilic parameters, including LogP and LogS; electronic parameters, including E_TOT, E_HOMO, and E_LUMO; and steric parameters, including molecular weight (MW) and molecular volume (MV), were calculated. The total clearance (CL_TOT) value as a parameter of drug elimination was predicted using online pkCSM (http://biosig.unimelb.edu.au/pkcsm/prediction). The QSPR equation was determined using Microsoft Excel 2013 (https://office.microsoft.com/Excel).

The Gaussian 09 program suite (Frisch et al. 2009) was used to optimize the molecular structures and compute the total energies, dipole moments, and molecular orbital energies of the molecules (Fig. 1). All calculations were performed using the most popular hybrid density, B3YLP (Becke 1993; Lee et al. 1988; Sousa et al. 2007), along with the 6-31++G (d, p) basis set. The harmonic vibrational frequencies were determined using the same method used for optimization to ensure that all molecules were local minima structures in the absence of imaginary frequencies. Using the HOMO and LUMO energies, several molecular properties were calculated: The HOMO-LUMO energy gap, ionization potential (I), electron affinity (A), chemical hardness (η), molecular softness (S), electronic chemical potential (μ), electronegativity (χ), and electrophilicity index (ω) were calculated.

Target proteins and prospective ligands have been developed for accurate computational computation. Ligand preparation was performed using the LigPrep tool interfaced with the Maestro module of the Schrödinger suite. Using the optimal potential liquid simulation (OPLS4) force field, the 3D structures of all the ligands and the reference chemical oxypurinol, including all potential tautomers and ionization states at pH 7.0 ± 2.0, were created and reduced. Schrödinger’s multi-step Protein Preparation Wizard (PrepWizard) was used for protein preparation. First, the RCSB Protein Data Bank was used to obtain a high-resolution protein crystal structure of XO (PDB: 3NVY) at a resolution of 2 Å, which contains a quercetin molecule in its binding site. The allocation of bond orders and charges was accompanied by the addition of hydrogen atoms to the heavy atoms, while all water molecules and heteroatoms were eliminated, ensuring that quercetin remained in the active site. The OPLS4 force field was employed to optimize and minimize the energy of the final structures for both the ligand and protein. A formal representation of the binding pocket was generated using the centroid of the co-crystallized native ligand (quercetin), with default settings maintained for each case. The Glide XP module of the Schrödinger Suite was utilized to dock the identified compounds and reference standards into the active sites of the crystal structures (Schrodinger 2021).

Quercetin was subsequently docked into the pre-prepared XO binding site to validate the binding site. The molecule’s most optimal docked pose was subsequently compared to the conformation of the bound ligand as well as the X-ray crystal orientation, and the root-mean-square deviation (RMSD) was calculated.

Based on the obtained physicochemical and ADME properties, the analyzed compounds showed different water solubilities that varied from very soluble, highly soluble, moderately soluble, and soluble (Table 1), indicating that they are less likely to have solubility problems and thus have good bioavailability (Bergström et al. 2016). Five compounds were found to have high gastrointestinal absorption after oral administration, and five compounds did not permeate the BBB. Only two compounds are P-gp substrates that protect the central nervous system (CNS) from xenobiotics. Cytochrome P450 (CYP) enzymes are responsible for the metabolism of many drugs. Therefore, drug-drug interactions and significant side effects resulting from drug or metabolite buildup are caused by the inhibition of these enzymes (Testa and Kraemer 2007). Therefore, an assessment of the effects of potential drug candidates on these enzymes is essential. Four compounds were predicted to be non-inhibitors of any of the five CYP isoforms (Table 1). The physicochemical properties were calculated using Lipinski’s rule of five, which states that molecules with logP 5, molecular weight 500 Da, ten hydrogen bond acceptors, and five hydrogen bond donors are most likely drug-like. Molecular entities with ten rotatable bonds and a total polar surface area of 140 are predicted to exhibit high bioavailability, as indicated by our findings in Table 1. These substances are also anticipated to have molecular weights below 500 Da, which facilitates their movement, dispersion, and absorption compared to heavier molecules (Srimai et al. 2013; Pires et al. 2015).

According to the Protox II website and the obtained LCD50 values, only four compounds were predicted to be nontoxic (Table 2).

A theoretical density functional analysis was conducted following the prediction of the ADME properties. The highest occupied molecular orbital (HOMO) energy, E_HOMO, is associated with the electron-donating ability of the molecule of interest. Consequently, molecules with higher EHOMO values possess an increased capacity to donate electrons to the unoccupied molecular orbital. The value of the lowest unoccupied molecular orbital (LUMO), ELUMO, is directly proportional to the molecule’s capacity to accept electrons. In other words, molecules with lower ELUMO values exhibit a greater propensity for accepting electrons. The significance of the HOMO and LUMO, as well as their respective properties, cannot be overstated when it comes to predicting the most reactive positions in π-electron systems and demonstrating various types of reactions in conjugated systems (Kurt et al. 2011). Furthermore, the HOMO and LUMO energy levels indicate the relative chemical stability and biological activity of the molecules. The difference in energy between the two frontier orbitals, known as the HOMO-LUMO energy gap, can be used to predict the strength and stability of the molecule. Molecules with a smaller HOMO-LUMO energy gap are more polarizable and are typically associated with increased chemical reactivity (Kurt et al. 2011; Gunase Karan et al. 2008). Cpds 3 and 4 have good HOMO-LUMO energy gaps, indicating their stability. Inspection of Table 3 reveals that Cpd2 has the lowest HOMO-LUMO gap among the molecules (3.9484 eV). Qualitative molecular representations of the HOMO and LUMO frontier orbitals and their corresponding energy values, in addition to the HOMO-LUMO gaps of the compounds, are depicted in Fig. 2. The HOMO and LUMO are delocalized over the entire cpd2 molecule. In contrast, the HOMO is delocalized at the conjugated part of the carbon chain of the molecule, whereas the LUMO is localized over the carboxylic group in cpd5. The HOMO and LUMO energies were then used to calculate global chemical reactivity descriptors, such as ionization potential (I), electron affinity (A), chemical hardness (η), molecular softness (S), electronic chemical potential (μ), electronegativity (χ), and electrophilicity index (ω). The results are summarized in Tables 1, 2. As shown in the table, the ionization potential of the studied molecules was in the range of 5.9919–7.5402 eV, with cpd7 and cpd5 having the highest and lowest ionization potentials, respectively. In the case of electron affinity, cpd5 had the lowest value of 0.3891 eV, and cpd4 had the highest value of 2.0708 eV. These values reflect the electron-donating and electron-accepting abilities of these compounds. Cpds5 and 1, which have carboxylic acid functional groups, exhibited higher values of chemical hardness and lower values of chemical softness than the other studied molecules. These compounds contain ionizable functional groups (-COOH) and tend to be in ionic form and weakly polarizable. The remaining compounds (2, 3, 4, 6, and 7) contained hydroxyl group(s) and, as such, had a lower chemical hardness than those of the carboxylic acid-bearing compounds (1 and 5).

Figure 2. 

3D HOMO and LUMO frontier orbitals of the (1–7) cpd compounds calculated at the B3LYP-31++G (d, p) level of theory.

Molecular electrostatic potential (MEP) is a graphical representation employed to delineate the spatial distribution of electron density and the three-dimensional charge distribution in a molecule. The determination of charge distributions is crucial in assessing the strength of van der Waals interactions between molecules. MEP has been established as a valuable metric for elucidating hydrogen bonding and correlating reactivity with molecular structure (Murray and Politzer 2011, 2017). In order to align MEP with the molecular surface, we employed GaussView (Dennington et al. 2009) and used the optimized structures of the compounds depicted in Fig. 1. Fig. 3 presents a three-dimensional plot of the molecular electrostatic potential. The specific colors chosen in this instance are intended to depict regions of attractive potential in shades of red and yellow, which are indicative of electrophilic reactivity and correspond to the negative potential of MEP. Conversely, regions of repulsive potential are represented in blue, symbolizing nucleophilic reactivity and the positive potential of MEP, making them an ideal location for electrophilic attacks.

Figure 3. 

The 3D plot of the molecular electrostatic potentials (in a.u.) formed by mapping the total density and the corresponding contours of Cpd (1–7).

The MEP of cpd1 indicates a low electrostatic potential for the C=O bond and a high electrostatic potential for the hydrogen atoms in O-H bonds; thus, these hydrogen atoms behave like acceptors, whereas the hydrogen atom in the N-H bond is characterized by a relative abundance of electrons, so these bonds behave like donors.

The resulting overall MEP of cdp2 indicates that the oxygen atoms carrying protons are not attractive to the negative test charge. The MEP of cpd3 and cpd4 showed a relatively low potential and yellow color, characterized by a relatively abundant number of electrons at the oxygen atoms. The 3D plots of MEP for cpd5, cpd6, and cpd7 revealed that the oxygen atoms of the carbonyl group were subject to nucleophilic reactivity. In cpd5, the H atom in the O-H bond is characterized by the absence of electrons; therefore, this hydrogen atom behaves like an acceptor.

To observe all the MEP surfaces, we plotted each surface as a contour around the molecule. Fig. 3 shows surface contour plots for each molecule. Each contour curve around the molecule is the MEP surface; the outer contour has a lower isosurface value, whereas the inner contour has a higher isosurface value. This explanation is similar to that of the MEP for the selected contour.

Following the DFT analysis, QSPR was used to correlate the physical parameters that can influence the clearance of the studied compound. QSPR analysis was conducted for the four compounds based only on their physicochemical properties, toxicity profile, and DFT analysis. The calculated chemical descriptors and QSPR results are summarized in Tables 4, 5, respectively.

The regression analysis identified the optimal regression equation for this study, which demonstrated the strongest correlation coefficients and significance values among the five equations examined. Specifically, the ideal regression equation was as follows (Widiyana et al. 2018):

Log 1/CL_TOT = 0.93 E_HOMO + 6.90

Because it had the best correlation coefficient (R = 0.98), the smallest significance of 0.02 (0.05), the smallest standard error (SE = 0.12), and the highest F value (41.45), it was chosen (Table 4). Thus, the electronic parameter E_HOMO is the physical chemistry parameter that influences CL_TOT, according to the equation.

To clarify the interactions between proteins and ligands, a molecular docking study was performed to explore the binding modes of the isolated compounds and xanthine oxidase. Re-docking of quercetin in the active site of xanthine oxidase was used to test the accuracy of the docking process. In the absence of water molecules, Fig. 4 shows docked quercetin and co-crystallized quercetin in nearly identical positions among the receptors (RMSD = 0.39), confirming the validity of the docking process using the EP scoring function. Table 6 summarizes the findings of this study, including the predicted glide scores and the binding interactions of the docked compounds. The obtained results were consistent with those of the DFT analysis.

Figure 4. 

Docked (yellow colored) and co-crystallized quercetin (green colored) in the active side of Xanthine oxidase.

Molecular docking of all compounds revealed that Kaempferol-3-O-glucoside and Esculin gave the best glide scores (-8.391 and -10.184 kcal/mol, respectively), even more than the reference compound, with a binding score of -1.654 kcal/mol.

We noticed that the two compounds formed a hydrogen bond with the Glu802 residue as the binding of the reference ligand and the native ligand, which also formed other important hydrogen bonds with the residues Ser876 and Thr1010.

Esculin, which had the best docking score, formed two hydrogen bonds with the Glu802 residue and one hydrogen bond with each of the Lys771 and Ser876 residues. Moreover, it developed hydrophobic interactions and one aromatic π–π stacking interaction with Phe1013, which explains its high Glide score and binding energy (Fig. 5).

Figure 5. 

3D and 2D binding disposition of compounds 3, 5, reference standard, and native ligand after docking calculations in the active site of XO.

Kaempferol-3-O-glucoside showed significant stability at the active site, with a docking score that was higher than that of the reference ligand. It forms an important hydrogen bond with the Glu802 residue as well as two hydrogen bonds with the Ser876 residue. In addition, it developed different hydrophobic interactions and two aromatic π–π stacking interactions with Phe1013 (Fig. 5).