The dataset of imidazole compounds was collected from prior reports/literature. Their structures were drawn in Chem Draw professional (Perkin Elmer), and 2-Dimentional structures were converted into 3-Dimentional structures using Chem3D Ultra (Version 19.1.0.8, Perkin Elmer). The value of enzyme inhibition (experimental activity) was expressed as (IC₅₀) for the training dataset which was then altered to its positive logarithmic scale using the formula: pIC₅₀ = −log(IC₅₀) and defined as a dependent variable. The database of compounds was generated in Microsoft Excel as a CSV output file (comma delimited).

To demonstrate a hypothesis for 3D conformation, the Field Template module of Forge V6.0 software was used as no structural data was attainable for imidazole derivatives in their target-bound state. For this purpose, the information about field and shape was utilized by the template from the library of 84 compounds. The hypothesis was developed by generating the three-dimensional field point pattern and calculating the field points of bioactive conformation.

At the connexion point of a 3D grid, the 3D-QSAR method calculates various molecular properties as molecular descriptors. This methodology covers the complete data of aligned training set compounds. The pharmacophore template was transported into the Forge V6.0 software, followed by the alignment of compounds with the associated template. After the alignment of 84 compounds with known IC₅₀ values, the 3D-QSAR model was built using the Field point-based descriptors. While building the model, the maximum distance of sample points was set to 1.0Å, the maximal number of components was set to 20, Y scrambles were adjusted to 50, and volumetric as well as electrostatic fields were also used. For overall resemblance, 50% dice volume similarity and 50% field similarity were achieved using the Forge software. The experimental activity (IC₅₀) of compounds was changed to pIC₅₀ which is equal to the negative log of IC₅₀. The set of 84 compounds was divided into the training set and test set with a ratio of 80% and 20%, respectively, and one compound was selected as the reference drug to assess the QSAR modeling using the activity stratified method.

The model was verified by q² (cross regression coefficient), r ² (regression coefficient), and similarity score of conformers for every ligand. LOO technique (leave-one-out) was used to assess the derived 3D-QSAR model. The LOO cross-validation technique is thought to be one of the most efficient techniques for the validation of the regression model having a small training dataset. The data size of N-1 was used for training and the remaining one was tested; N identifies the complete dataset. In the LOOCV technique, the process of testing and training was repeated for the N number of times, and in this way, each data was passed through the testing method. Then, the test data which is not in the training set is used to derive the 3D-QSAR model.

The global aspect of training data was studied quantitatively by using the Bayesian approach. The hydrophobicity, electrostatics, and shape attributes, which lie beneath the SAR of a particular set of compounds, are better understood by this approach. These 3D models were viewed to achieve valuable information. The three types of interconnected biochemical evaluated data including regions explored analysis, activity cliff summary, and average of actives were revealed by the Activity-atlas study. The regions explored analysis exhibited the areas of aligned and fully explored compounds. The details about negative and positive electrostatic sites, appropriate and inappropriate hydrophobicity, and appropriate shape of actives were provided by the activity cliff summary. On the other hand, an average of actives helped in showing the common parts in active compounds which were selected.

The 3D structures of target proteins (PDB ID: 7CHM, 3PP0, 4UDD, 5NWH, 3HY3, and 4ZVM) ( Supplementary Table S2 ) were downloaded from the RCSB PDB database (https://www.rcsb.org/). The protein preparation was performed to accomplish various tasks such as identifying the active site, deleting alternate conformations, interpolating missing atoms in incomplete residues, protonating titratable residues, modeling the missing loop areas, and removing the water molecules and heteroatoms (Alam and Khan, 2019). The ligands of proteins were used to identify the active sites from the “Define and Edit Binding Site” option in Discovery Studio, and SBD_Site_Sphere was generated (Figure 3).

Active ligands along with reference Fulvestrant (standard drug) having known inhibitory potentials were collected from literature, and SDF files of some of the compounds were downloaded from PubChem while others were generated from the Chem3D software where they were optimized through MM2 and MMF9 force field. For ligand preparation, the “Open Babel” software was used. For this purpose, the input file was imported in the form of sdf–MDL MOL format and the output path was selected in pdb format. Then the ligand was converted into the desired form.

PDBQT files of proteins and ligands were prepared, and by using the Graphical User Interface program AutoDock Tools, grid box creation was accomplished. Fragmental volumes to the protein, polar hydrogens, united atom Kollman charges, and solvation parameters were assigned by AutoDock tools. After preparation, all files were saved as PDBQT. Preparation of the grid map was done by using a grid box via AutoGrid. The grid size was selected as 60 × 60 × 60 XYZ points with a grid spacing of 0.375Å. Default settings were used for all other parameters. Autodock was used for docking protocol and information about proteins and ligands was used along with grid box features in the configuration file. Both ligands and proteins are considered rigid when using Autodock. The results lower than 1.0Å in the root-mean-square deviation were assembled and were depicted by the result with the most suitable free-binding energy while the results with the lowest binding affinity or binding energy were extracted and subjected to further analysis (Azam and Abbasi, 2013).

A three-dimensional structure-activity relationship (3D-QSAR) was performed to throw more light on a series of imidazoles. For this purpose, a conformational hunt was carried out on these compounds (C1-C84). (Alam and Khan, 2017). The three-dimensional pattern of field points (Figure 4) was identified by illustrating the derived conception of bioactive conformation with its calculated field points. Four distinct molecular fields were calculated, namely, negative and positive electrostatic potential, hydrophobicity, and shape/van der Waal descriptors. To draw a pharmacophore template resembling the bioactive conformation (for further virtual screening), a molecular field-based similarity approach was employed.

Figure 4 shows that the compounds with similar filed points bind at the same target site. This characteristic provides a linear correlation between biological activity and structural similarity of ligands (Low et al., 2005). Based on the parameter, Figure 4 shows the active ligands, and their similarity metric was found within the range of 56% to 42%.

The ligand alignment in the protein context is required to use the 3D similarity metric for activity atlas model development. To ensure accurate model development, this alignment must be inspected. The compounds in the training set were aligned to ensure that the molecules being compared were in the same relative orientation (Figure 5). This alignment is necessary because molecules can adopt different conformations or spatial arrangements due to the freedom of rotation around single bonds. After it, the 3D-QSAR model was built by using the Field points-based descriptors. The activity interactive graph plot was used to represent the fitness of the derived 3D-QSAR model. This graph displays the comparison between predicted and actual activity with cross-validation data points. Fairly good activity-descriptors’ relationship accuracy of 81% was achieved by the derived 3D-QSAR model as the regression coefficient was r ² = 0.81. Similarly, as mentioned by the cross-validation regression coefficient (q² = 0.51), a high activity-prediction accuracy of 51% was attained. The derived 3D-QSAR model was proved to be very reliable to predict the anticancer and cytotoxic activity of imidazole derivatives as an MCF7 cell-line inhibitor (Table 1).

Figure 5 shows that there were displayed little spaces in ligands to accommodate small conformation changes and variations of moieties present in aligned molecules. High molecules have tight alignment to restrict the substitution or replacement of any group present in them. Low-active molecules, in comparison, lack steric tightness and have the capacity to substitute any moiety in the context of an activity enhancer (Low et al., 2005).

The QSAR model was also viewed in a 3D form to unveil the structure-activity relationship (SAR) mechanism of imidazole derivatives. The field points named coefficient and variance (associated with the bioactivity of training set compounds) were analyzed in a 3D structural form for the purpose. The derived model points for QSAR were contrasted with the reference compound for better comprehension of space field point localization (Alam and Khan, 2014). In a robust model, the high coefficient and variance field points were proved to be the highly essential correlating parameters. According to the results, electrostatic and steric coefficients both play a major role in modulating the anticancer activity as represented by the large size of red, cyan, green, and pink field points (Figure 6). Field points containing high steric and electrostatic variance indicated regions of high changes while the field points containing low variance represented the regions with less changes or no changes (Figure 6).

“View field contributions to predicted activity” study was done on imidazole derivatives. This evaluates the extent to which imidazole derivatives fit the derived field-based 3D-QSAR model and regions of structural field points governing the predicted activity. These field contributions (Figures 6A, B) were represented by purple, blue, and red color regions. According to the results, the orange- and purple-colored areas denote the region of electrostatics and steric field points, respectively, having the negative regulation capability on predicted activity (decrease anticancer activity). Whereas the green- and zinc-colored areas denote the regions of electrostatics and steric field points, respectively, with a positive regulation capability on the predicted activity (increased anticancer activity).

SAR study was practiced through the activity-atlas visualization technique and was used to unveil the key features of imidazole, regulating the anticancer activity and designing more novel drugs. For this purpose, an activity cliffs summary and an average of actives study were performed on imidazole derivatives.

On the basis of this model, imidazole compounds having a pIC₅₀ value higher than 6.4 were classified as active compounds while the rest of the compounds were considered inactive. This model (Figure 7; Fig. S1) represents the areas of high activity that reference drugs and active ligands have in common.

In Figure 7A, the positive and negative electrostatic regions represented by the red color sites and cyan color sites, respectively, correlate with the anticancer activity, i.e., more of these regions indicate more anticancer activity. Figure 7B shows the hydrophobic interaction regions of active compounds as indicated by the yellow color, and Supplementary Figure S1 shows the shape regions of active compounds as represented by the white color.

The activity cliff summary diagrams as indicated by Figure 8A correlate with the biologically active parts of imidazole drugs with the reference drug. The cliff summary of electrostatics in Figure 8A is visualized in two colors: red and cyan. The presence of the red color indicating the positive electrostatic field and the cyan color indicating the negative electrostatic field is favorable for high anticancer activity. Figure 8B shows the areas of favorable and unfavorable hydrophobics represented by the green and purple colors, respectively. Whereas in Figure 8C, the green color shows a favorable shape region, and the purple color shows an unfavorable shape region.

The descriptive features of compounds were explored in this model aside from their biological activities (Attiq et al., 2022). The more red and cyan colors indicating positive and negative fields, respectively, show the areas of strong SAR with the reference drug. The average regions explored also represent the areas of active compounds that would not take part in an anticancer activity (Figure 9).

Molecular characteristics regulating the active compounds as anticancer agents were retrieved for further prediction of their anticancer activity based on derived SAR models. Before that, prediction performance was analyzed on the test set and training set compounds by predicting their anticancer activity. This prediction was done by means of derived models and then the distance value (error) was compared. In order to perform the comparison, predicted activity and distance to models’ columns were examined for each derived model. The important ligand fields were illustrated for each derived model through this study and after that, these characteristics were utilized for virtual screening.

SBVS of selected compounds was performed to discover new valuable drugs in order to treat breast cancer (El Aissouq et al., 2021). AutoDock Tools provided notable results with overall binding energy of all the selected compounds ranging from −5.1 to −10.5 kcal/mol. Most of the compounds have binding free energy greater than the reference compound Fulvestrant when observed with all the selected proteins as shown in Table 2.

The docking results explain that compounds C6, C10, C21, and C76 are the most hit compounds giving excellent results, so only these compounds will be subjected to further study.

Molecular Interaction and Binding Mode. The top hit compounds represented by the shaded area in Table 2 were selected to evaluate the binding site interactions between the ligand and the target protein.

Compound C6 was fixed in the binding pocket of protein (PDB ID: 4ZVM) (Figure 10A) by undergoing electrostatic interactions (Pi-Anion) with GLU193 and hydrophobic interactions (Pi-Sigma, Pi-Pi Stacked, and Alkyl) with PHE17, TYR104, and PRO192 (Supplementary Table S3).

The docking of compound C6 with protein (PDB ID: 5NWH) is described in Supplementary Figure S2. C6 undergoes interactions with the protein by hydrophobic interactions (Pi-Sigma and Pi-Alkyl) with VAL 49, VAL29, and PRO86, Pi-Cation electrostatic interactions with ARG196, and other interactions (Pi-Sulfur) with CYS91 (Supplementary Table S3).

Supplementary Figure S3 describes the docking of compound C10 with the protein having PDB ID: 7CHM. The compound fits in the binding pocket of the protein through hydrogen bonding (Conventional Hydrogen Bond and Carbon Hydrogen Bond) with amino acids ASP608, SER611, GLN530, ASP674, and ILE607; hydrophobic interactions (Pi-Sigma, Pi-Pi Stacked, Alkyl, and Pi-Alkyl) with ILE531, ILE607, LEU654, ALA651, PRO673, VAL539, ILE663, and ALA551; and miscellaneous interactions (Pi-Sulfur) with MET602, CYS604, and MET671 (Supplementary Table S4).

However, the docking of compound C21 with PDB ID: 3HY3 provided different results (Supplementary Figure S4). C10 interacted by forming a conventional hydrogen bond with TRP109 and hydrophobic interactions (Pi-Pi Stacked, Pi-Pi T-shaped, Alkyl, and Pi-Alkyl) with TYR83, TRP109, PRO81, MET90, TYR152, TYR153, and LYS150 (Supplementary Table S4) (Elancheran et al., 2023a).

Compound C21 fits in the binding pocket of protein (PDB ID: 4UDD) (Supplementary Figure S5) through hydrogen bonding (conventional hydrogen bond and carbon-hydrogen bond) with GLN642, GLN738, and PRO637, hydrophobic interactions (Pi-Pi Stacked) with TRP557 and TYR735, and miscellaneous interactions (Pi-Sulfur) with MET745 (Supplementary Table S2). Docking results of Compound C76 with PDB ID: 3PP0 are described in Figure 10B. C76 interacted by forming conventional hydrogen bond with LYS753 and hydrophobic interactions (Pi-Sigma and Pi-Alkyl) with LEU726, VAL734, THR798, LEU852, ALA751, and LYS753 (Supplementary Table S2).

The RMSF graph provides insight into the flexibility of individual atoms or residues in the protein. It shows how much they deviate from their average positions during the simulation. The maximum value of RMSF indicates greater flexibility, while the smallest value denotes the system’s restricted motion across the simulation course. In all of our proteins, the RMSF graph showed a number of areas having high flexibility (Figure 12). The complex of C10 with H3Y3 showed a maximum number of peaks, indicating more flexible movements. All three docked complexes showed a maximum RMSF of 1.0. It indicated that the atoms or residues are, on average, deviating from their average positions by around 1 Å (Angstrom).

The B-factor, eigenvalue, deformability, covariance matrix, variance map, and elastic network model of the protein serves as a representation of its stability. The mobility characteristics of the docked proteins are determined by the deformability and B-factor. The peaks are associated with the protein regions with deformability whereas the areas with the highest peaks are those with the greatest deformability (Santra and Maiti, 2022). In B-factor graphs, the comparison between the PDB field and NMA of the docked complexes is provided. The B-factor graphs of 3H3Y-C10 and 7CHM-C10 complexes showed that the PDB data predicted higher B-factors compared to the NMA data. It suggested that the B-factor values predicted by the computational simulations using NMA showed lower mobility or flexibility than what was predicted by the experimentally determined B-factor values from the Protein Data Bank. Figure 13 illustrates the deformability and B-factor of 3HY3-C10, 4UUD-C10, and 7CHM-C10, respectively.

The eigenvalue indicates motion stiffness which is directly connected to the energy needed to deform the structure. If the eigenvalue is low, it means that the complex is more easily deformable. The eigenvalues of 3H3Y, 4UUD, and 7CHM complexes with C10 are 3.231587e-04, 7.624966e-05, and 2.479516e-04, repectively. It means that all of our docked complexes showed low eigenvalues, indicating a considerable amount of deformability and, hence, good flexibility and stability of the molecular motion. The individual variance is shown by purple-shaded bars in the variance graph of C10 with target proteins, while the bars with green shading show cumulative variance. The eigenvalue and variance graphs of protein-ligand complexes (H3Y3-C10, 4UUD-C10, and 7CHM-C10) are shown in Figure 14.

The covariance matrix indicates correlations among the pairs of residues in a protein-ligand complex (Figure 15). The red and white colors showed correlated and uncorrelated motion, respectively, while anticorrelations are represented by the blue color. Greater correlation means the formation of a better complex. The covariance matrices for H3Y3-C10, 4UUD-C10, and 7CHM-C10 complexes exhibited good correlations and minimal anticorrelations.

The elastic network model of docked proteins shows relationships between the atoms where the stiffer regions are indicated by the darker grey areas (Figure 16). All protein elastic maps yielded reliable results.

Toxicity results provided valuable information related to the toxicological profile of selected compounds after molecular docking studies and thus may be useful for drug designing (to select the dosage and preferred route of administration (Banerjee et al., 2018). However, all these results are preliminary and must be confirmed by experiment (Table 3).

As shown in Table 3, among all the hit compounds, C10 was found to be inactive in all toxicity parameters and, thus, it was considered the best-hit compound.

Very useful information about the compounds can be provided by Frontier molecular orbitals (FMO) such as electronegativity, stability, reactivity, and chemical hardness and softness (Elancheran et al., 2023b). The HOMO and LUMO parameters are used to compute the chemical reactivity descriptors and to assess the molecular reactivity (Al-Janabi et al., 2021). The energy values of HOMO and LUMO were determined by the DFT method as shown in Supplementary Table S6. Contour diagrams of FMOs are shown in Figures 17–20.

The geometry optimization was done to minimize the energy and to find the most stable atomic arrangement. Optimized geometries of all compounds are shown in Figures 17–20 along with the numbering system and the vector of dipole moment.

The DFT calculations for ionization energy, electron affinity, energy gap, electronegativity, chemical potential, electrophilicity index, chemical softness and hardness, additional electronic charges, maximum charger transfer index, and dipole moment were performed as explained in Supplementary Table S6).

A high E_HOMO indicates that the molecule is a strong electron donor and can easily donate electrons to the receptor, increasing biological activity. According to this concept, the activity ranking of the hit compounds is given below with an increasing E_HOMO value. C 10   >   C 21   >   C 76   >   C 6

A low E_LUMO value indicates that the compound can easily accept electrons from the donor molecule, increasing biological activity (Mu and Gao, 2022). According to this criterion, an increase in the biological activity of compounds is as follows: C 6   >   C 76   >   C 10   >   C 21

The third parameter is the energy gap (ΔE) between Homo and Lumo. If the energy gap is small, it indicates that the molecule is soft, biologically active, less stable, and has a high chemical reactivity. In other words, the biological activity increases with a decrease in the energy gap. The order of ranking should be C 21   >   C 10   >   C 76   >   C 6

A high chemical potential (CP) or Lower electronegativity X) value indicates electron delocalization. It means that the molecule can easily form bonds and coordinate easily with the biological system. So according to our calculated DFT data, an increase in the biological activity of compounds is given below: C 21   >   C 10   >   C 76   >   C 6

The dipole moment also affects the biological activity of the compound. The high value of the dipole moment indicates the strong ligand-protein interaction, thus increasing the biological activity (Sayin and Üngördü, 2018). Supplementary Table S6 shows that the dipole moment of our hit compounds are in the order of C 10   >   C 21   >   C 76   >   C 6

In order to determine a chemical mechanism, MEP maps and MEP contours play a very important role. Molecular electrostatic potential helps to determine the hydrogen bonding interactions and to interpret the nucleophilic as well as electrophilic reactions (Horchani et al., 2020). MEP can be used to indicate the shape of the molecule and the sizes of the negative, positive, and neutral electrostatic potential. The molecular structure of drugs along with the interaction among different physicochemical properties can be predicted by the MEP scale (Al-Janabi et al., 2021). The MEP of compounds C6, C10, C21, and C76 is determined under the basis set of B3LYP/6-311G. The negative charge is indicated by red and yellow areas that represent the electrophilic attack sites. The green color indicated a neutral charge while the blue region indicating the positively charged areas represented the nucleophilic reactivity (Bendjeddou et al., 2016). The MEP structures and MEP scales of compounds are shown in Figures 17–20.