After a comprehensive literature review, 50 chemical compounds of S. album (Liu et al., 2008; Umdale et al., 2020) were retrieved in sdf format from the Pub-Chem database (Kim et al., 2021). The compounds were then cleaned, and energy was minimized using Avogadro software v1.2.0 with the help of mmf94 force field (Hanwell et al., 2012). The scaffold structure of the hit compounds shown in Figure 2.

The x-ray crystallography structure of SarA protein of S. aureus (PDB ID: 2FNP) was retrieved from the protein data bank (Bouley et al., 2015). With the aid of the software Discovery Studio v21.1.0.0 (https://discover.3ds.com/discovery-studio-visualizer-download), the protein structure was initially cleaned and heteroatoms were removed. Then, using the GROMOS96 43b1 force field and the SwissPDB Viewer software v4.1 (https://spdbv.unil.ch/disclaim.html) the energy of the cleaned protein was minimized and optimized (Guex and Peitsch, 1997). The active site and binding pocket of the target protein SarA (PDB: 2fnp) was identified by CASTp online server (Figure 3).

The molecular docking between phytochemicals from S. album and SarA protein of S. aureus was carried out using PyRx software v0.8 (https://sourceforge.net/projects/pyrx/) which works based on the Auto Dock vina configuration (Dallakyan and Olson, 2015) and docking was performed based on the previous methods with small modification (Muhammad and Fatima, 2015; Bhowmik et al., 2021; Dey et al., 2023). The protein structure was converted into a macromolecule and the ligands were converted to PDBQT format. The center and grid box size of the docked complexes were set as for center X: 3.8347 Å, Y: −0.1554 Å, Z: 8.7332 Å and for dimension X: 48.7233 Å, Y: 41.3090 Å, Z: 59.7350 Å respectively. The final docking calculation was conducted using PyRx and top molecules were selected based on lower binding energy. The binding interactions and poses were explored via Discovery Studio software.

YASARA (Yet Another Scientific Artificial Reality Application) Dynamics software v19.12.4 was used to conduct molecular dynamic simulation with the help of the Assisted Model Building with Energy Refinement (AMBER)14 force field (Supplementary Material S1) (Wang et al., 2004; Land and Humble, 2018). The hydrogen bond network was initially cleaned up and optimized together with the docked complexes. In order to reduce the protein complexes using a TIP3P water solvation model (0.997 g/L1, 25 c, 1 atm), the steepest gradient approaches were used (Harrach and Drossel, 2014). The physiological conditions were set at 0.9% NaCl, 310K, and pH 7.4 (Krieger et al., 2012) to neutralize the simulated system. The simulation time step was set as normal at 2 × 1.25 frames per second. The long-range electrostatic interaction was calculated by the particle mesh Ewald (PME) method with a cutoff radius of 8.0 Å (Essmann et al., 1995). The simulation trajectories were saved after every 100 ps and the final simulation run was conducted for 100 ns (Krieger and Vriend, 2015). The root-mean-square deviation, the solvent-accessible surface area, the radius of gyration, and hydrogen bonding were all assessed using the simulation trajectories (Baildya et al., 2021; Islam et al., 2021).

Calculating the binding free energy is a crucial approach for evaluating the strength of the interaction between a drug and a protein. This analysis provides insights into the energetic aspects of the drug-protein complex. To determine the binding free energy, various snapshots of the complex were analyzed using the MM-Poisson–Boltzmann surface area (MM-PBSA) method in YASARA software. The calculation involved the following formula:

MM/PBSA binding free energy = EpotReceptor + EsolvReceptor + EpotLigand + EsolvLigand−EpotComplex−EsolvComplex.

The calculations were performed using the AMBER 14 force field, and YASARA macros were utilized for efficient computation of the MM-PBSA binding energy.

Principal component analysis (PCA) was utilized to explore the overall variability among protein-ligand complexes, including comparisons with the apo form and a drug-protein complex, by considering all the structural features (Akash et al., 2023b; 2023a). This method allows for the identification and categorization of structural changes in protein-ligand complexes that occur during the simulation by comparing different variables. Through diagonalization of the covariance matrices and solving the eigenvalue and eigenvector problems, PCA was performed on the complexes. The eigenvalues provided information about the magnitude and direction of structural fluctuations, while the eigenvectors represented the directions of these fluctuations. The MD trajectories spanning 100 ns were pre-processed and standardized by scaling to unit variance and removing the mean (Ichiye and Karplus, 1991; Shukla and Tripathi, 2020). The PCA technique was implemented using Python v3.11 with the Scikit-learn v1.2 library, and the results were visualized using Matplotlib v3.7 (Supplementary Material S2).

The phytochemicals that passed the docking study were then subjected to absorption, distribution, metabolism, excretion, and toxicity (ADMET) analysis to see whether they have the necessary qualities to be considered as lead molecules. Therefore, the pkcsm (Pires et al., 2015) and SwissADME (Daina et al., 2017) web servers were used to analyze ADMET profiles and calculate molecules’ adherence to Lipinski’s rule of five respectively.

The pharmacophore mapping analysis of the top three ligands was performed using the online server PharmMapper (Wang et al., 2017). The ligands were obtained in sdf format from the PubChem server and then uploaded to the server. During the process, the “maximum number of conformations” parameter was set to 1,000. All available targets were selected in the “select target set” parameter, and the “number of reserved matched targets” parameter was set to 1,000. The fit score cut-off value in the advanced settings was set at 0. The default settings were used for all other parameters.

S. aureus is associated with skin cancer by secreting virulence factors which are regulated by SarA protein (Bromfield et al., 2023). Interestingly, there was no published study has been reported on the inhibitory activity of S. album phytochemicals on SarA protein. Therefore, in this current investigation, we performed in silico docking study targeting SarA (PDB ID: 2fnp) protein of S. aureus which is responsible for the expression of many virulence genes (Jiang et al., 2023), by using S. album derived phytochemicals. The active site residue of SarA protein are A:Phe110, A:Ser114, A:Thr117, A: Thr118, A:Lys121, A:Glu223, A:Leu224, B:Thr141, B:Thr142, B:Glu145, B:Asn146, B:His159, and B:Tyr162 (Table 1). In our study, we found all three compounds strongly bind with the most of the active site residues of SarA protein including A:Phe110, A:Thr117, A:Lys121, A:Glu223, A:Leu224, B:Asn146, and B:His159 (Table 2), which indicates that these three compounds can potentially significant for the inhibition of the target protein SarA.

In a molecular docking study, from 50 compounds of S. album top five compounds were chosen based on the lowest binding energy, interaction with target protein, pose and RMSD value. Compound CID: 5280441 found a maximum binding energy −9.4 kcal/mol followed by CID: 162350 and CID: 5281675 with binding energy of −9.0 kcal/mol and −8.6 kcal/mol respectively compared to the positive control ciprofloxacin CID: 2764 with −8.6 kcal/mol (Table 2).

Compared to the previous study our current compounds showed better binding affinity than compound hesperidin which interacts with active sites of SarA with the binding energy of −6.9 kcal/mol by revealing two hydrogen bonding interactions at Thr 117 and Lys 163 (Tong et al., 2015). Another study reported that benzimidazole type NHC precursors 1a-d molecules exhibited binding energy ranging from −5.04 to −5.46 kcal/mol with the SarA protein which is comparatively lower to our current findings and did not show any kinds of hydrogen bond (Üstün et al., 2021), the difference of the binding energy maybe due to presence of more hydrogen bond interaction between the protein and our current ligands (Figure 4). Some other studies also found that the SarA protein of S. aureus was inhibited by different compounds (de Oliveira et al., 2019; Cheruvanachari et al., 2023). However, in our study five hydrogen bonds were found from complex 2fnp + CID: 5280441 at B: Asn146 (2.29 Å), A:Leu224 (2.16Å), A:Glu223 (2.66Å), A:Lys121 (2.73Å) and B:His159 (2.66 Å) (Figure 4A). Similarly, Complex 2fnp + CID:162350 was stabilized by four hydrogen bonds at A:Glu221 (2.06 Å), A:Glu223 (2.24 Å), A:Lys121 (2.77Å) and B: His159 (2.40 Å) (Figure 4B). Moreover, the interaction between 2fnp and compounds CID: 5281675 revealed four hydrogen bonds in A chain of 2fnp protein at A:Lys121 (2.70 Å) and B chain of 2fnp protein at B: Asn146 (2.46 Å), B: Asn158 (2.18 Å) and B:His159 (2.07 Å) (Figure 4C). Additionally, the complexes 2fnp + CID: 114776 and 2fnp + CID: 5281654 showed almost similar binding energy with −7.7 and −7.4 kcal/mol respectively (Table 2), and also interact target protein with three hydrogen bonds where two residue such as B:Asn146 and B:His159 were same with different distance (Figures 4D, E).

In comparison with other three lead compounds, the positive control CID:2764 exhibited four hydrogen bond at A:Thr117 (2.38 Å), A:Asp120 (3.28 Å), B:Tyr162 (5.22 Å), B:Gln166 (3.79 Å). Overall, in silico docking study revealed that the three best compounds from S. album exhibited better inhibitory activity than positive control ciprofloxacin. This finding suggests that the binding interactions may serve as a potential mechanism accountable for the inhibition of S. aureus associated skin cancer.

Molecular dynamics (MD) simulations are performed on protein-ligand complexes to gain a detailed understanding of their dynamic behavior at an atomic level (Shukla and Tripathi, 2020). Which, enables a dynamic view of their behavior, providing valuable information for drug discovery, understanding biological processes, and elucidating structure-function relationships (Guterres and Im, 2020). In this molecular dynamics simulation, the three best docked complexes was performed to validate the docking conformational stability and its rigidity at 100 ns time-dependent manner which enables the uncovering of potent inhibitors. The root means square deviations (RMSD), solvent accessible surface area (SASA), the radius of gyration (Rg), and the hydrogen bonds of the SarA protein of Staphylococcus aureus and the best docked ligand complex are shown in Figure 5.

The Root Mean Square Deviation (RMSD) is a widely used metric for evaluating structural variations and determining the stability, precision, and conformational changes in protein-ligand complexes (Singh et al., 2023). In this study, we conducted an analysis of the RMSD values for the complex formed between drug candidate compounds and the SarA protein over a simulation period of 100 ns. The average RMSD value for complexes C1, C2, C3, Apo and PC was 2.57 Å, 2.07 Å, 2.03 Å, 2.10 Å, and 1.18 Å respectively (Table 3). According to the simulation result, the C1 complex showed the highest RMSD value between 20 and 40 ns which was 2.5 Å to 4.0 Å and exhibited increased RMSD value after 20 ns by following the upward movements. This might be due to the conformational variability of C1 complex after 20 ns. On the other hand, other two complexes C2 and C3 maintained a stable RMSD value from the very beginning to the rest of the simulation period which was almost similar to the RMSD profile of Apo protein and positive control (PC) ciprofloxacin (Figure 5A). The results of the analysis indicate that complexes C2 and C3 exhibited higher structural stability and maintained a consistent conformation throughout the simulation period in comparison to the C1 complex. The lower and more stable RMSD values observed for C2 and C3 suggest that they performed better in terms of overall structural integrity and robustness. Based on the RMSD findings, it can be concluded that C2 and C3 possess enhanced stability compared to the C1 complex.

The radius of gyration (Rg) is a crucial metric that provides insights into the overall size of a protein-ligand complex (Miraz et al., 2023). It is a fundamental measure used to assess the structural fluctuations occurring during molecular dynamics (MD) simulations. By analyzing Rg, we can quantitatively evaluate the degree to which the protein’s structure changes the simulation (Singh et al., 2021). The average Rg values for each complexes ranging from 22.05 to 22.43 (Table 3). The Rg profile of the three complexes showed a similar trend with the Apo protein and positive control throughout the whole simulation period ranging from 21 to 23 Å (Figure 5B). A lower Rg value is commonly associated with a compact and rigid structure (Islam et al., 2021). This suggests that all complexes are positioned relatively closer to their center of mass, indicating a folded or globular conformation, and did not change throughout the 100 ns simulation period.

The Solvent Accessible Surface Area (SASA) analysis is employed to determine the surface area of a molecule that is accessible to solvent molecules that is accessible to solvent molecules, can provide insights into the stability and folding of proteins (Baruah et al., 2022). The average SASA values for each complexes ranging from 13712.75 to 13832.02 (Table 3). In this study, the SASA profile of three complexes showed decreasing trend from the beginning of the simulation and after 20 ns had a stable profile that was almost similar to the Apo protein and ciprofloxacin (PC) (Figure 5C). It suggests that the complexes, Apo protein and ciprofloxacin share similar levels of accessibility to the surrounding solvent environment. This similarity in SASA values may imply similar levels of flexibility, exposure to functional sites, or potential for interactions. This implies that the ligand-binding regions in all complexes likely serve comparable roles or functions.

The analysis of hydrogen bonds in molecular dynamics (MD) simulations of protein-ligand complexes is crucial for understanding the nature and stability of their interactions (Takano et al., 2022). Hydrogen bonds play a vital role in determining the specificity and strength of binding between the protein and ligand (Wohlert et al., 2022). By monitoring hydrogen bond formation and breaking events during the simulation, one can gain insights into the dynamics of the complex, including transient interactions and the stability of the binding interface. The hydrogen bonds of C1, C2, and C3 complexes were similar to Apo protein and positive control (PC) ciprofloxacin, and did not change during the 100 ns simulation period (Figure 5D). These similar hydrogen bonds of all complexes suggest that the binding interactions between the protein and ligand are consistent and maintained throughout the simulation. It indicates that specific amino acid residues in the protein form stable hydrogen bonds with the ligand over time, contributing to the overall stability of the complex.

In our research, we also employed RMSF (Root Mean Square Fluctuation) analysis to investigate how individual atoms or residues in a biomolecular system behave in terms of flexibility and dynamics. This analysis also helped us pinpoint specific residues responsible for these fluctuations (Akash et al., 2023b). Accroding to the RMSF results, the average values ranging from 1.42 to 1.67 Å which is outstanding results (Table 3). Among the top three compounds we studied, their RMSF profiles closely resembled that of the positive control, ciprofloxacin, as well as the Apo protein. However, we observed that certain residues exhibited significant fluctuations, exceeding 4 Å, specifically residues numbered 85–89 and 204–210 (Figure 6). Despite the overall similarity in the RMSF profiles of these compounds, some residues within the complexes displayed a higher degree of flexibility and mobility during our simulations. This increased flexibility in these particular residues may have important functional implications, such as aiding in binding events, accommodating structural changes, or playing a role in molecular recognition processes.

The MMPBSA binding energy represents the net energy change associated with the formation of a protein-ligand complex compared to the unbound components. In our study, the binding free energy of protein-ligand interactions was determined using Molecular Mechanics/Poisson-Boltzmann Surface Area (MMPBSA) calculations. The MMPBSA binding free energies of all complexes were evaluated at a simulation period of 100 ns, as depicted in Figure 7. The average MMPBSA binding free energies for complexes C1, C2, C3 and PC were found to be −309.08 ± 1.92, −186.90 ± 2.40, −292.07 ± 2.02, and −215.193 ± 1.81 kJ/mol, respectively, herein it has been shown that the MMPBSA binding energy for C1 and C2 were better that the positive control (Table 4). This negative MMPBSA binding energy indicates a favorable binding interaction, indicating that the binding of the ligand to the protein is energetically favorable and stable.

The inclusion of structural information, both in the docked and simulated states is shown in Figure 8. The superimposition revelaed that all three lead compounds exhibited strong binding activity at the binding pocket whereas both docked and simulated complexes showed almost similar structure compared to the positive control.

The MD simulation snapshot of complexes C1, C2, C3 and PC at 0, 25, 50, 75, and 100 ns are shown in Figure 9. According to this results, both C1 nad C2 compounds remain sustain at the binding pocket of the target protein 2fnp throughout the simulation period where significant structural change were not observed for this both complexes (Figure 9). Similarly, C3 compounds also almost suatain with the target protein however there was slight structural changes of the protein was observed particularly at 75, and 100 ns. On the other hand, the positive control drug was splited out from the target protein at 50 ns however at rest of the time control drug was bounded with target protein. Overall, superimposion results suggest that all three compounds were strongly bounded with the target protein throughout the simulation period compared to the positive control.

PCA is commonly used in MD simulations to analyze and understand the conformational dynamics and structural variations of apo form and protein-ligand complexes (Paris et al., 2014; Zarezade et al., 2018; Kitao, 2022). It was conducted to gain insights into the alterations in conformational dynamics that occur upon the binding of the ligand. PCA helps uncover the most significant structural changes and collective motions occurring in the protein-ligand complex during a molecular dynamics (MD) simulation (Kaur Bijral et al., 2022). In our study, PCA was utilized to explore the overall variability among protein-ligand complexes, including comparisons with the apo form. The techniques employed included the diagonalization of covariance matrices and mathematical eigenvalues, which provided information about the magnitude and direction of structural fluctuations. Furthermore, the eigenvectors represented the directions of these fluctuations (Alom et al., 2023).

Figure 10 showcased PCA results derived from all the trajectories incorporating both coordinate information and pertinent structural features. The total explained variance ratio of Apo, C1, C2, and C3 was 35.89%, 39.51%, 38.74%, and 35.42% correspondingly. The conformational distributions of the apo form appear to be similar to those of the protein-ligand complexes. The results of the PCA analysis indicated that there is a widespread distribution of the apo form and protein-ligand complexes, suggesting that they exhibit conformational stability throughout the trajectory. Additionally, the distribution of the protein-ligand complexes closely resembled that of the Apo form, indicating a similarity in their overall structural characteristics. Moreover, the complexes did not display substantial structural variations throughout the MD simulation. The presence of low variance eigenvectors further confirmed the stability of the conformational states, suggesting that despite their dynamic nature, the complexes maintained consistent structural configurations.

ADMET (absorption, distribution, metabolism, excretion, and toxicity) analysis was performed for top three ligands to find out a lead compound. We found that CID: 5280441 (C1), CID: 162350 (C2) and CID: 5281675 (C3) followed Lipinski’s rule although having few violations. The bioavailability score of C1, C2, and C3 compounds was 0.55, 0.55, and 0.17 respectively, which indicates that these compounds were physiologically active as the bioavailability score of a compound determines its physiological activity (Hosen et al., 2023). According to the solubility scale, insoluble < −10 < poorly < −6 < moderately < −4 < soluble < −2 < very <0 < highly water soluble (Daina et al., 2014), the water solubility of compounds C1 (−2.845), C2 (−2.812) and C3 (−2.905) was found to be water-soluble. Compounds C2 showed the highest human intestinal absorption rate of 64.729% followed by C1 (46.695%) and C2 (43.733%). The Blood brain barrier and CNS permeability of all three compounds exhibited negative results (Table 5), indicaing that compounds are less likely to cross the BBB and less likely to permeate the central nervous system (CNS) as well (Shaker et al., 2023). Moreover, three screened compounds showed positive results in the human ether-a-go-go (hERG) l inhibitor test, and no toxicity was found in hepatotoxicity and AMES tests (Table 5). It suggests that these three compounds are suitable for further lab experiments. So, these compounds could be used as lead compounds to develop a drug against antibiotic resistant S. aureus.

In pharmacophore mapping, the fit score is a quantitative measure used to assess the degree of similarity or fit between a pharmacophore model and a given molecular structure or compound (Opo et al., 2022). Herein, we analyzed the pharmacophore mapping of all three compounds for SarA protein (Figure 11), and the fit scores of compounds C1 and C3 are 4.51 and 3.701, respectively, which are higher than that of C2 with 2.95 (Table 6). These higher fit scores of the compounds indicate a closer match between the pharmacophore model and the compounds, suggesting a higher likelihood of the compounds exhibiting similar biological activity to the target or reference molecule. On the other hand, the normalized fit score provides a standardized measure of the similarity between a compound and a pharmacophore model, enabling comparisons across different models and compounds in pharmacophore mapping studies (Thangavel and Albratty, 2022; Tyagi et al., 2022).

In our study, we revealed the highest normalized fit scores for compounds C1 and C3, with 0.7375 and 0.7403, respectively, followed by C2 with 0.6443 (Table 6). The normalized fit score of a compound is usually between 0 and 1, where 1 represents a perfect match between the compound and the pharmacophore model. This finding suggests that C2 and C3 were almost similar, and both compounds match the pharmacophore model well. Additionally, compounds C1 and C3 generated similar hydrophobic centers and one hydrogen bond acceptor each. On the other hand, C2 exhibited two hydrophobic centers with one hydrogen bond acceptor. Interestingly, except for C1, none of the other compounds exhibited a hydrogen bond donor. Moreover, all compounds C1, C2, and C3 did not generate any positively charged centers, negatively charged centers, or aromatic rings (Table 6).

To validate the findings of our in silico analysis, we conducted an in vitro assessment of the antibacterial activity of three key compounds, namely, Vitexin (CID:5280441), Isovitexin (CID:162350), and Orientin (CID:5281675) from S. album against S. aureus is shown in Figure 12. The result revealed that all compounds exhibited almost similar level antibacterial activity. The compounds Vitexin showed the strong activity with zone of inhibition 19.8 ± 0.11 mm followed by Isovitexin and Orientin with zone of inhibition 18.53 ± 0.27 mm and 18.16 ± 0.08 mm respectively compared to the positive control ciprofloxacin with 22.90 ± 0.05 mm. On the other hand, negative solvent control doesnot effect on the growth of S. aureus (Supplementary Figure S1). Hence, based on the comprehensive in vitro investigation, it can be inferred that these three compounds possess the potential to inhibit S. aureus associated with skin cancer. Nonetheless, further molecular studies are imperative to provide additional confirmation and insight into their mechanisms of action.