A compound library of 271 derivatives resembling MTX and PTX was employed in this study (also look at Supplementary Table S1). These compounds were specifically chosen as structural analogues of MTX and PTX, aiming to explore their potential as anticancer agents. The selection was based on their structural similarity to established anticancer drugs, focusing on compounds with analogous features and characteristics.

The scheme of study with the implementation of the DL model is shown in Figure 1. The library combined contained 271 structural analogues (as shown in the Supplementary Table S2) of MTX and PTX. The compounds library was prepared using DFT-optimized geometries and converted to pdbqt format for virtual screening using Auto Dock Vina (Trott et al., 2010). The crystallographic structure of the protein of interest was obtained from the Protein Data Bank (PDB) (https://www.rcsb.org/), accessed on May 3^rd, 2023. The target PDB identification code for the structure is 1U72 (Human Dihydrofolate Reductase Receptor). Subsequently, the macromolecules were prepared using MGL tools, involving the elimination of heteroatoms and water molecules, as well as the introduction of polar hydrogen atoms. The protein was assessed for any missing amino acid residues. Kollman’s charges were computed to neutralize the protein, ensuring its overall neutrality, while Gasteiger’s charges were calculated for the ligand, aligning with common practice in molecular simulations (Weiner et al., 1984; Dermawan et al., 2021). To conduct a virtual screening of the compound library using Auto Dock Vina, grid box size (x = 80, y = 80, z = 80) and dimensions (x = 31.10, y = 14.21 and z = −8.07) were adjusted for XYZ coordinates to cover the active binding pocket during the virtual screening using Auto Dock Vina. The process of virtual screening was conducted after the preparation of the specific protein, employing the script-based approach of Auto Dock Vina. The value for exhaustiveness was configured as 8, while the number of nodes was specified as 30. The virtual screening procedure was conducted twice to validate the precision of the docking outcomes.

The docking protocol was rigorously validated through a redocking procedure using the same algorithm and parameters employed for the investigated compounds. This validation method ensures the ability of the molecular modeling simulation to accurately replicate the ligand binding mode and residue-wise interaction patterns. Specifically, achieving a root mean square deviation (RMSD) value below 2 angstroms between the native and redocked poses confirms the reliability and biological significance of the adopted docking protocol (Katari et al., 2016). Following the conclusion of the virtual screening process, the resultant findings from the virtual screening module were subjected to analysis and docking scores of the drug candidates were compared to those of the standard drugs MTX and PTX.

DL models were used to predict drug-target interactions (DTI) in the current work. The training data consists of a diverse dataset of drug compounds and their respective protein targets. This dataset was carefully curated from publicly available sources and annotated with high-quality labels to ensure data integrity. Including diverse compounds, targets, and proteins helps ensure the model’s predictions apply to various scenarios. The model architecture was based on encoder and decoder architectures, widely used in deep learning for sequence-to-sequence tasks. This choice of architecture was made due to its effectiveness in handling the complex relationships between drug compounds and protein sequences. The encoder processed the Simplified Molecular Input Line Entry System (SMILES) string of the drug compound, while the decoder handled the amino acid sequence of the target protein. Such architecture allows the model to capture intricate features and interactions between drugs and targets. For performance evaluation, a rigorous and well-established protocol was followed. The dataset was split into training, validation, and test sets to assess the model’s generalization ability. The model was trained using over 17 state-of-the-art deep learning techniques, validated in various studies for their efficacy in DTI prediction. Measure the model’s performance using metrics such as precision, recall, F1-score, and receiver operating characteristic (ROC) curves. This extensive performance evaluation approach ensures the model’s predictions are robust and reliable (Huang et al., 2020a; Aziz et al., 2022a; Imran et al., 2023). This study utilized the DeepPurpose architecture (http://deeppurpose.sunlab.org/), a machine-learning approach for predicting the interaction between ligands and target proteins. The pre-trained model, MPNN-CNN-Binding Data Base-IC₅₀, which combines the Message Passing Neural Network (MPNN) with Convolutional Neural Network (CNN) techniques, was used to predict drug-target interactions (DTI) in terms of binding affinity (IC₅₀) and predicted pIC₅₀ (Huang et al., 2020a; Aziz et al., 2022a; Imran et al., 2023). This model incorporates advanced neural network approaches to enhance the accuracy of DTI predictions.

The inherent ability of AI to revolutionize the drug discovery process has already been established by the uptake and frequent use of absorption, distribution, metabolism, excretion, and toxicity (ADMET) predictive tools, virtual screening, and quantitative structure-activity relationship (QSAR) modeling (Leelananda and Lindert, 2016). Development and use of in silico QSAR models to forecast drug activity, has become increasingly targeted in drug development and drug discovery over the near past (Ekins et al., 2007). Simple ML algorithms like multiple linear regression (MLR) exhibit efficient utilization for model development of relatively small data sets. QSAR modelling aims to construct a model with strong robustness and predictive capacity (Veerasamy, 2011; Hammoudan et al., 2023).

The FDA-approved drugs MTX and PTX have demonstrated efficacy in cancer treatment by effectively inhibiting the DHFR pathway (Yamashita et al., 2020). The binding affinity (IC₅₀) and ADMET profile predominantly influence the drug’s effectiveness (El-Adl et al., 2021). Consequently, DL models have been utilized to forecast the IC₅₀, pIC₅₀, and ADMET characteristics of the most promising candidates obtained via virtual screening (Al-Jumaili et al., 2023b; Choudhary et al., 2023). This approach aims to facilitate a direct evaluation of the binding affinity of these candidates regarding the established standards of MTX and PTX. Applying in silico methods for predicting binding affinity and ADMET characteristics presents a promising alternative to experimental approaches (Brogi, 2020). In the present study, DL models were employed to predict interactions between drugs and their respective targets, commonly referred to as drug-target interactions (DTI). These models were constructed based on the encoder and decoder architectures (Yu et al., 2022). The DL model scheme uses the SMILES string (Table 1) and amino acid sequence of the specific protein of interest as its input, employing more than 17 cutting-edge DL techniques to forecast indicators of drug efficacy (Figure 1). In this study, the researchers used the MPNN-CNN deep learning algorithms to predict affinity and specifically applied the MPNN model for ADMET predictions (Huang et al., 2020a).

DFT methods are a reliable and proficient approach for the correct estimation of the electronic features of the compound (Hohenberg and Kohn, 1964). The complete geometry optimization for all the Zinc15 library database (https://zinc15.org/) compounds of MTX and PTX as shown in Supplementary Table S1 was conducted by applying the DFT methodology at the ground state level, employing three functional parameters. The Becke, Lee-Yang-Parr (B3-LYP) functional (Prieto-Martínez et al., 2018; Jemai et al., 2023) was used, consistently incorporating the GD3 correction for dispersion (DiLabio et al., 2013). The absence of imaginary frequencies observed during a harmonic vibrational analysis confirmed the existence of local energy minima in their character (Yang et al., 2019). The obtained results were utilized to calculate various electronic properties, such as the energy gap (Eg), for the top hits under investigation. The energy gap can be juxtaposed with specific molecular attributes, such as reactivity and electrical conductivity. The energies of the lowest unoccupied molecular orbital are commonly denoted as E_LUMO, while those of the highest occupied molecular orbital are denoted as E_HOMO. In addition to the aforementioned electronic properties, the values of physiochemical descriptors such as hardness (η), softness (S), chemical potential (µ), and electrophilicity index (ω) were determined employing Koopman’s theorem (Tsuneda et al., 2010). The relationship between chemical stability and reactivity can be elucidated by considering the parameters of chemical hardness (η), electronegativity (ꭓ), and softness (S) (El-Shamy et al., 2022). The basis set for all examined structures was 6-311+g(d,p). The calculations were executed using the Gaussian09 software suite (Frisch, 2009), while the Chemcraft suite was used to visualize the HOMO-LUMO Frontier molecular orbitals (Kanagathara et al., 2022). The Gauss View utility was also employed for visualization (Hanwell et al., 2012).

MD simulations provide information about the stability of the best complex (Rasheed et al., 2021). The following steps were employed to execute the MD simulation. The 3-dimensional (3D) models of the enzyme Human DHFR (PDB ID: 1U72) (Rose et al., 2010) in complex with methotrexate molecule were exported to the.pdb format using Pymol (DeLano, 2002). The dynamic behavior of the complexes was assessed through MD simulation using the GROMACS software package (version 2022.2) (Bekker, 1993; Ganesan et al., 2017; Thalla et al., 2020). The protein topology was generated using the CHARMM27 force field (Schmid et al., 2011) through the pdb2gmx tool. The ligand topology was also created using the SwissParam server (Van Aalten et al., 1996). The complexes were subsequently introduced into the system upon implementing the force field. The solvation of the system was carried out using the TIP3P water model (Mark and Nilsson, 2001). A cubic box with dimensions greater than 1 nm from the protein’s edge was employed and periodic boundary conditions were applied. The system was rendered inert by introducing Na⁺ ions, followed by the execution of energy minimization for a total of 50,000 steps utilizing the steepest descent algorithm. Subsequently, a 100 picosecond (ps) NVT simulation at a temperature of 300 Kelvin (K) was conducted, followed by 100ps NPT simulation to achieve equilibrium for the entire system. The Leapfrog algorithm was utilized in the constant-temperature, constant-pressure (NPT) ensemble to independently couple each component, including the protein, ligand, water molecules, and ions (Van Gunsteren and Berendsen, 1988). The temperature and pressure coupling constants for the Berendsen method were assigned values of 0.1 and 2, respectively. These values were chosen to maintain a stable environment for the system, with a temperature of 300K and a pressure of 1 bar (Berendsen et al., 1995). The molecular dynamics (MD) simulation was conducted for 100 nanoseconds under isothermal and isobaric conditions in an ensemble at 300 Kelvin. The time constant for pressure coupling was configured to 1 picosecond to ensure a constant pressure of 1 bar.

Additionally, the LINCS algorithm (Hess et al., 1997) was employed to enforce constraints on bond lengths. The Van der Waals and Coulomb interactions were truncated at a distance of 1.2 nm. The PME algorithm (Di Pierro et al., 2015), integrated into GROMACS, was employed to minimize the error resulting from truncation. The trajectory files were visualized using Visual Molecular Dynamics (VMD) 1.9.2 (Humphrey et al., 1996) and analyzed using the in-house developed tool HeroMD Analysis (Devi et al., 2021; Rawat et al., 2021) and Xmgrace 5.1.25 (Vaught, 1996).

This study has yielded significant findings in the realm of virtual screening and QSAR profiling utilizing DL techniques, alongside an evaluation of ADMET properties. The detailed results, including Supplementary Tables S1, S2, showcase the comprehensive outcomes of our investigation.

A total of eight compounds represented in Table 1 exhibited docking scores surpassing those of conventional drugs. The highest-ranking results underwent additional examination utilizing DL algorithms and QSAR. Thus, the utilization of DL models facilitated the prediction of drug affinity and the assessment of protein-ligand complex stability (Rasheed et al., 2021; Aziz et al., 2022a).

The structures of SBVS are presented in Figure 2, which highlights the most significant hits. It is important to note that each occurrence exhibits a shared pharmacophore with established MTX and PTX, recognized as inhibitors of DHFR. A common pharmacophore indicates that these compounds can bind to the active site of DHFR in a way that is analogous to MTX and PTX.

The implementation of a shared pharmacophore with established DHFR inhibitors presents numerous benefits. Initially, it creates a foundation for these compound’s logical development and improvement to augment their inhibitory efficacy and specificity towards DHFR. Through comprehension of the fundamental structural components accountable for binding, alterations can be implemented to enhance the potency and effectiveness of the compounds. In addition, a common pharmacophore among these compounds implies that they may operate through comparable pathways such as MTX and PTX, thereby reinforcing their viability as DHFR inhibitors for addressing DHFR-linked malignancies.

Figure 3 represents the 3D protein structure of 1U72 and catalyic binding pocket. The molecular docking results of the top hits with the lowest binding energy are illustrated in Table 2. The 3D hydrogen bonding interactions of top hits with residual amino acids within the binding pocket of 1U72 are shown in Figure 4 whereas the 2D interactions of top hits and reference compounds are presented in Figures 5–7.

The Figure 5 and Table 2 illustrate that compound 9 exhibited binding interactions with amino acid residues VAL08, ALA09, VAL115, and TYR121, which are also known to interact with MTX and PTX. It is worth noting that compound 9 displayed a greater number of H-bonding interactions with their corresponding distance in angstroms, reflecting that ALA9 (2.256 Å), LYS54 (2.2.798 Å), THR56 (2.991, 2.959Å), VAL115 (3.255 Å), VAL115 (3.302 Å), TYR121 (2.767 Å), and VAL8 (3.620 Å), thus suggesting the possibility of better anticancer potential compared to other selected top hits. However, it is important to acknowledge that these findings are based on in silico analysis, and further experimental validation would be necessary to confirm the anticancer efficacy of compound 9. The sulfonyl oxygen in compound 27 formed H-bonds with ALA9 (2.841 Å), VAL8 (3.436 Å), and GLY17 (3.735 Å), and aromatic rings displayed π-π/π-σ interactions with ILE16, PHE34, ILE60 and LEU67 suggested a favorable binding mode for compound 27. In comparison, MTX and PTX exhibited similar H-bonding interactions with amino acid residues, emphasizing the potential of compound 27 as a promising anticancer candidate. The sulphone moiety of compound-41 could not interact with the receptor due to the steric effect of adjacent ring systems. However, the ring system is engaged in H-bonding interactions with GLU30 at a distance of 2.170 Å and hydrophobic interactions with ALA9, ILE60, and PRO61 demonstrating a favorable binding profile for compound 41. In comparison, MTX and PTX primarily formed H-bonds with specific residues, highlighting the distinct binding characteristics of compound 41. Compound 68 formed H-bonds with ASP145 at the relevant distance of 3.17 Å. PHE34, ILE16, PHE34, and ALA9 depicted the hydrophobic interactions thus displaying favorable binding interactions. These H-bonding interactions contribute to its potential as an effective anticancer agent. MTX and PTX also exhibited H-bonding and hydrophobic interactions, demonstrating similar binding patterns to compound 68.

The results in Figure 6 and Table 2, Compound 74 showed hydrophobic interactions with ALA9, VAL115, PHE34, LEU22, and ILE16, indicating its ability to interact favorably with the target residues. The hydrophobic interactions observed for compound 74 suggest a potential role in disrupting the hydrophobic core of the target protein. MTX and PTX similarly engaged in hydrophobic interactions, demonstrating comparable binding characteristics to compound 74. Compound 85 formed H-bonds with THR56 and LEU22 at the distance of 3.027 and 4.249 Å, highlighting its favorable binding interactions. Additionally, it exhibited hydrophobic interactions with PHE34, LEU22, ILE16 and LYS55. These binding features suggest a promising potential for compound 85 as an anticancer agent. MTX and PTX also demonstrated H-bonding and hydrophobic interactions, indicating similar binding profiles. Compound 99 was found to be engaged in H-bonding interactions with ASP145 (3.378 Å), THR56 (2.961 Å), and GLY17 (3.439Å) with different bond distances, suggesting a favorable binding affinity. It also displayed hydrophobic interactions with GLY17, PHE34, LEU22, ALA9, ILE16, and LYS55. These interactions contribute to the potential of compound 99 as an effective anticancer candidate. Comparatively, MTX and PTX exhibited similar H-bonding and hydrophobic interactions, aligning with the binding characteristics of compound 99. The amino acid residues TYR 121 LYS68 interacted with compound 180 by electrostatic interactions whereas GLN35 (3.269 Å), showed H-bonding interactions with aromatic amino nitrogen. The hydrophobic interactions via pi-sigma bonding and π-π stacking with ILE16, LEU22 and ILE60 are crucial for favorable electron correlation. Our analysis unveiled a robust connection between the binding modes of studied ligands, MTX and PTX, highlighted by hydrogen bonds that play a crucial role in stabilizing interactions. In addition, the distances and angles of the hydrogen bonds were thoroughly analyzed to confirm the optimal hydrogen bonding for both ligands. By aligning these interactions with the known binding mode of MTX and its conformation, the biological significance of the binding poses of studied ligands and PTX has been established. This explanation offers valuable insights into the structural foundation of ligand-protein interactions and highlights the potential therapeutic significance of studied ligands and PTX as a new version of MTX. Overall, these interactions indicate a unique binding profile for compound 180 and exhibited binding similarities with MTX and PTX. MTX and PTX binding interactions are shown in Figure 4 (3D interactions) and Figure 7 (2D interactions).

The findings from the molecular docking analysis demonstrated the varied binding interactions between the top-hit compounds and the specific amino acid residues of the target. The observed compounds exhibited a synergistic interplay between H-bonding and hydrophobic interactions, suggesting their promising utility as efficacious agents for combating cancer. Although MTX and PTX displayed different binding patterns, the top-hit compounds exhibited comparable or superior binding characteristics. This implies that the compounds with the highest activity show considerable promise as innovative agents for combating cancer, thus justifying the need for additional research and refinement.

The current study employed the multiple linear regression (MLR) method to conduct regression analysis (Liu and Long, 2009), as illustrated in Table 3. As shown in Supplementary Table S1, the dataset was partitioned into separate training and test sets. The analysis was conducted using the substitution groups and the inhibitory activity of compounds within the dataset. The training dataset was utilized for constructing the QSAR model, while the test dataset was employed to assess the predictive ability of the developed models. A dataset comprising 271 compounds, shown in Supplementary Table S1 from the Zinc15 library, was used in a DL model to forecast ADMET properties. Subsequently, the DL-based ADMET predictions were employed to construct a model.

The selection process was guided by a focus on structural similarity to established anticancer drugs, emphasizing compounds with analogous features and characteristics. Lipinski’s rules were applied to evaluate the resemblance of the investigated molecules to pharmaceutical compounds (Abdelrheem et al., 2020). A drug candidate that adheres to no more than one of the rules above will probably be pursued for development as a potential oral medication. Furthermore, we assessed the drug-like properties of the recently developed molecules using deep-learning ADMET techniques. This analysis aimed to identify the molecular structures exhibiting favorable oral drug administration characteristics. To validate this choice, an in silico assessment was conducted to evaluate the pharmacokinetic parameters of the selected compounds using ADMET prediction (Rifaioglu et al., 2019).

By comparing the binding affinities of the compounds to the DHFR cancer target protein, it is possible to observe variations in their predicted IC₅₀ values and pIC₅₀ values, as shown in Table 4. Compound 68 exhibits notable characteristics concerning binding affinity, as evidenced by its predicted IC₅₀ value of 78.66 nM and pIC₅₀ value of 7.10. The data indicates that compound 68 exhibits a robust and efficacious binding affinity towards the DHFR oncogenic protein. Compound 74 demonstrates a notable binding affinity, as evidenced by its anticipated IC₅₀ value of 62.46 nM and pIC₅₀ value of 7.20. This observation suggests a robust binding affinity between the compound and the target protein, similar to compound 68.

In contrast, it can be observed that compounds 27, 41, and 180 exhibit moderate binding affinities, as their anticipated IC₅₀ values fall within the range of 1,394.69 nM to 1,423.03 nM, thereby yielding pIC₅₀ values that vary from 5.85 to 5.93. Although the binding affinity of these compounds is comparatively weaker than that of compounds 68 and 74, they display significant interactions with the DHFR protein, which is a target for cancer treatment. Compounds 09, 85, 99, MTX, and PTX exhibit reduced binding affinities, as indicated by their predicted IC₅₀ values, which range from 129.91 nM to 2,187.48 nM and correspond to pIC₅₀ values ranging from 5.80 to 6.89. The studied compounds display lower binding affinities towards the target protein than those above. To summarize, the binding affinities of the compounds 68 and 74 are the highest, while compounds 27, 41, and 180 exhibit moderate affinities. Compounds 09, 85, 99, MTX, and PTX demonstrate reduced binding affinities with DHFR target protein when they were subjected to the DL prediction model.

Table 5 presents a comparative analysis of the ADMET (absorption, distribution, metabolism, excretion, and toxicity) outcomes generated by the DL, highlighting notable patterns and variations across the compounds currently under study using the MPNN model. The solubility values, expressed in logarithmic units of mol/L, range from −3.93 to −5.45. The solubility of the compounds can be inferred from their respective negative values, where higher solubility is associated with more negative values. Compound 180 has the highest solubility with a value of −5.45, while compound 09 has the lowest solubility of −3.93. As a result, it is possible that compound 180 exhibits superior dissolution properties in comparison to the remaining compounds. The range of lipophilicity, as determined by the logarithmic ratio, falls between 0.19 and 2.69. The lipophilicity of compound 74 is the highest among the tested compounds, as indicated by its value of 2.69. This suggests that it has a greater tendency to dissolve in lipids or fats, with higher values indicating increased fat solubility. On the other hand, it can be observed that compound 09 displays the least degree of lipophilicity (0.19), thereby suggesting a comparatively reduced capacity for dissolution in fats.

The absorption percentages, specifically related to human intestinal absorption (HIA), transportation by P-glycoprotein (Pgp), and Bioavailability-F20, have been included in the study for reference. These values indicate how well a compound is absorbed in the human intestines, the extent to which it is transported by P-glycoprotein, and its potential bioavailability, with a focus on Bioavailability-F20. These parameters are valuable for assessing a compound’s suitability for further development and its potential as a drug candidate. The values mentioned above denote the proportion of the substance assimilated following ingestion via the oral route. The oral absorption rates of the compounds were found to be high, as evidenced by the absorption (HIA) values ranging from 93.6% (compound 09) to 99.12% (compound 99). The absorption rate of compounds across the gastrointestinal tract through P-glycoprotein, an efflux transporter, exhibits a range of variability from 12.5% (compound 74) to 33.3% (compound 41). The bioavailability of the administered compounds, ranging from 76.3% (compound 09) to 78.83% (compound 99), pertains to the extent of absorption and represents the proportion of the dose that enters the systemic circulation. Distribution (BBB) and Distribution (PPBR) percentages are presented as distribution values. The distribution (BBB) parameter denotes the capacity of the compound to traverse the blood-brain barrier, with a range of values observed between 43.47% (compound 180) and 74.68% (compound 85). Elevated values indicate superior distribution across the blood-brain barrier. The distribution of a compound, as measured by its plasma protein binding rate (PPBR), indicates its affinity for binding to proteins in the bloodstream. The observed PPBR values for the investigated compounds range from 35.3% (compound 09) to 94.11% (compound 180).

The metabolic rates of specific cytochrome P450 enzymes (CYP2C19, CYP2D6, CYP3A4, CYP1A2, CYP2C9) that participate in drug metabolism are indicated as percentages in metabolism values. The metabolic profile of Compound 41 marks a significant metabolism rate by CYP3A4 (28.66%) and CYP1A2 (79.34%), implying that these enzymes are playing a crucial role in the compound’s metabolic pathway. The results indicate that Compound 180 undergoes considerable metabolism by CYP2C9, with a rate of 61.74%, suggesting a significant contribution of this enzyme in the metabolic process of this compound. The column “half-life” indicates the duration required for the concentration of a given compound within an organism to decrease by 50%. Compound 6 demonstrates a brief half-life of 6.53 h, implying a swift elimination rate. Conversely, compound 68 displays the lengthiest half-life of 8.33 h, indicating a more gradual elimination process.

The clinical toxicity values indicate each compound’s toxicity rating or score in the last column compound 180 exhibits the most favorable clinical toxicity ratings, with a value of 8.3%, respectively. Conversely, compounds 09, 27, 41, 68, and 99 demonstrate clinical toxicity ratings ranging from 31.64% to 35.2%. The clinical toxicity ratings of compounds MTX and PTX are 31.78% and 30.12%, respectively, suggesting the presence of potential toxicity. Methotrexate, in particular, is known to be a highly toxic drug, classified as a cytotoxic agent. The level of toxicity can vary significantly based on factors such as dosage, duration of administration, individual patient characteristics, and other relevant determinants. Additionally, toxicity can indeed vary in terms of its impact and severity on patients. In brief, these DL model findings presented in Table 5 indicate that compound 180 manifests the most outstanding solubility, compound 74 displays the highest lipophilicity, compound 99 exhibits the highest rate of oral absorption, compound 85 demonstrates the highest distribution across the blood-brain barrier, compound 41 undergoes substantial metabolism by CYP3A4 and CYP1A2, compound 68 has the lengthiest half-life, and compounds 180 exhibit the least clinical toxicity ratings.

The physiochemical descriptors of the top hit compounds obtained through virtual screening and the standard anticancer drugs MTX and PTX were determined using DFT analysis and are presented in Table 6. The aforementioned descriptors in the table offer valuable insights into the electronic properties and reactivity of the compounds, thereby playing a crucial role in comprehending their behavior and potential applications. The electronic energy values obtained for the compounds ranged from −2,810.59 atomic units to −1,079.78 atomic units. Electronic energy is a comprehensive measure of the energy the electrons possess within a given system, and the molecular structure and composition of the compound determines its magnitude (Closs et al., 1986). The dipole moment, denoted in Debye units (D), is a quantitative measure of the spatial displacement between the positive and negative charges within a given molecule (Minkin, 2012). The dipole moment values presented in Table 6 exhibited a range spanning from 0.92 D to 9.31 D, which can be attributed to the polarity and asymmetry characteristics of the compounds. The E_Homo values, which denote the highest occupied molecular orbital energy, exhibited a range from −6.24 eV to −5.50 eV. These values determine compounds’ electron donation ability, where more negative E_Homo values indicate a greater propensity to donate electrons in chemical reactions. The values of E_Lumo, which denote the energy of the lowest unoccupied molecular orbital, exhibited a range spanning from −2.73 eV to −1.62 eV.

The E_Lumo plays a crucial role in the acceptance of electrons. Compounds exhibiting lower E_Lumo values demonstrate a heightened capacity for electron acceptance. The energy gap, which is determined by the disparity between the E_Homo and E_Lumo, explains the electron transfer capabilities of the compound. The observed energy gap values range from 2.85 eV to 4.21 eV, signifying the compounds’ diverse capacities for electron donation or acceptance, as shown in Table 7. The ionization potential is a fundamental concept in chemistry, which can be mathematically expressed as the negative of the highest occupied molecular orbital energy (E_Homo). It signifies the minimum amount of energy needed to extract an electron from a molecule that is in a neutral state. This study’s range of ionization potential values spanned from 5.50 eV to 6.24 eV. This range suggests that the compounds examined in this research exhibit varying degrees of stability and resistance to remove electrons. The electron affinity of a molecule is determined by the energy change that occurs when an electron is acquired, and it is quantified as the negative value of the energy of the lowest unoccupied molecular orbital (E_Lumo). The observed electron affinity values ranged from 1.62 to 3.09 eV, signifying the compounds’ inherent capacity to accept electrons. Supplementary Table S3 showing the optimized XYZ coordinates for top-hit anticancer compounds and reference drugs.

The chemical descriptors of the compounds were virtually screened as top hits, and the standard anticancer drugs MTX and PTX were assessed using DFT calculations. The results are presented in Table 8, which showcases the outcomes obtained from the DFT calculations. These descriptors provide valuable understandings of the compounds’ electronic structure, reactivity, and stability, elucidating their chemical behavior. The chemical potential (μ) is determined by the arithmetic mean of the ionization potential (I) and the electron affinity (A), and it signifies the inclination of a compound to either donate or accept electrons. The chemical potential values in Table 8 ranged from −4.52 eV to −3.65 eV. These values indicate the compounds’ electronic stability and capacity to engage in electron transfer phenomena. Electronegativity (χ) measures an element’s ability to attract electrons, calculated as the average ionization potential and electron affinity. It is a quantitative indicator of a compound’s capacity to attract electrons. The compounds exhibited a range of electronegativity values, ranging from 4.52 eV to 3.65 eV, indicating their diverse propensities to attract electrons during chemical reactions. The chemical hardness (η) is determined by taking half the difference between the ionization potential and the electron affinity. It characterizes the ability of a compound to resist alterations in its electronic configuration. The chemical hardness values in Table 8, range from 1.43 eV to 2.11 eV, suggesting the compounds’ inherent stability and capacity to endure electronic perturbations. Global softness (σ), which is the inverse of chemical hardness, offers valuable insights into the reactivity and vulnerability of a compound to electron transfer. The global softness values exhibited a range spanning from 0.48 eV⁻¹ to 0.70 eV⁻¹, whereby lower values were indicative of heightened electron transfer resistance and chemical reactivity. A compound’s electrophilicity (ω) can be determined by evaluating the square of its chemical potential divided by twice its chemical hardness. This parameter measures the compound’s capacity to accept electrons and engage in electrophilic reactions. The electrophilicity values exhibited a range spanning from 6.64 eV to 10.22 eV, thereby indicating the diverse propensities of the compounds to function as electron acceptors. Most of the selected top hits exhibited ionization potential, electrophilicity, and electron affinity values similar to those of the reference drugs, suggesting they behave electrochemically similarly to MTX and PTX under physiologic conditions (Avendaño, 2008). However, compounds 27 and 180 showed greater ionization potential and electrophilicity, which may contribute to improved interactions and potentially enhanced therapeutic efficacy. The literature also showed similar results of related studies with the DFT method use (Graffner-Nordberg et al., 2000; Zia et al., 2019; Aziz et al., 2022a).

The physicochemical and chemical descriptors of the top hit compounds identified through virtual screening and the standard anticancer drugs MTX and PTX were thoroughly examined using DFT analysis. The discoveries above contributed to an enhanced comprehension of the compounds’ electronic characteristics, reactivity, and stability, which are pivotal in assessing their potential efficacy as anticancer agents. The results obtained from this study thus established a valuable resource for future drug design efforts, as they contribute to the development of new compounds that exhibit improved therapeutic effectiveness and targeted anticancer properties and also contributed in alignment with the existing literature (Hobani et al., 2017; Mammen et al., 2020; Rana et al., 2020).