We conducted extensive computational analyses to gain insights into the characteristics of the synthesized naphthoquinone derivatives utilizing Gaussian16 software (Frisch et al., 2019). Molecular geometry optimization, executed at the B3LYP/6-31G (d,p) level of theory, yielded optimized structural parameters, as elaborated in Supplementary Tables S1–S16. In a prior investigation, Uddin et al. (2023) comprehensively explored the structural parameters, including bond lengths and angles, for seven compounds. The calculations were conducted at the B3LYP level of theory, utilizing various basis sets, namely, 6-31G(d,p), 6-311G(d,p), and 6-311++G(d,p). Their findings indicated that B3LYP/6-31G(d,p) exhibited good agreement with experimental results. Furthermore, Frontier Molecular Orbital (FMO) analysis was undertaken to elucidate the energy levels and spatial distributions of the HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital). These Frontier orbitals, positioned at the outermost boundary of the molecule, play a pivotal role in determining its electronic behavior. The HOMO energy reflects the molecule’s electron-donating ability, while the LUMO energy characterizes its electron-accepting potential. An ideal scenario is characterized by a high HOMO value and a low LUMO value, indicating favorable electron transfer properties. Gaussview 6 software (Dennington et al., 2016) was employed to visualize the electrostatic potential distribution across the molecules. Molecular electrostatic potential (MEP), represented through a spectrum of colors, reveals the net electrostatic potential arising from the total charge distribution at nucleophilic and electrophilic sites within a compound. Electrophilic regions, indicating the strongest attraction, are depicted in blue, while nucleophilic areas, signifying the highest repulsion, are colored red. Green zones represent a neutral electrostatic potential. Additionally, the dipole moments of the compounds were calculated using Gaussview 6. The FMO energy gap, representing molecular stability, is determined by the difference between HOMO and LUMO energies. The ligands’ overall chemical reactivity descriptors, including ionization potential (IP), chemical potential (µ), electron binding energy (H), maximal charge acceptance (∆Nmax), global chemical hardness (η), global chemical softness (σ), energy change (∆E), electrophilicity (ω), electronegativity (χ), and electron affinity (EA) (Alberty, 1960; Chamizo et al., 1993; Parr et al., 1999; Elkaeed et al., 2022), were computed using these HOMO and LUMO values. The formulas for these descriptors are as follows: E Gap   eV = E LUMO   −   E HOMO ;  IP  eV = − E HOMO ;  EA  eV = − E LUMO ;   μ   eV = IP + EA / 2 χ = −   μ ;   η = IP  −  EA / 2 ;  and  σ = 1 / η ;   ω = μ 2 / 2 η

The swissADME server tool (www.swissadme.ch) (Daina et al., 2017) was employed to assess the effectiveness of napthoquinone derivatives (1–13) and determine their physicochemical properties. Data were gathered from PubChem, the research database https://pubchem.ncbi.nlm.nih.gov/search/search.cgi (Kim et al., 2021), and a previous research paper. This tool provides accurate predictive results for identifying physicochemical properties, medicinal chemistry, and drug-likeness (Daina et al., 2017). Additionally, the AdmetSAR tool (http://lmmd.ecust.edu.cn/admetsar2/) (Yang et al., 2019) and ADMET predictor software were utilized to analyze ADMET properties. Using each compound’s canonical simplified molecular input line entry system (SMILES), various values, including multiple toxicity values, CYP inhibitors, and hERG pIC50 values, were determined. AdmetSAR offers both unrestricted and restricted values for each ADMET property. Furthermore, the Molinspiration online open-access server (https://www.molinspiration.com/cgi-bin/properties) (Molinspiration cheminformatics software, Choice reviews online. 2006) was employed to evaluate the correlation between the physicochemical properties and molecular activity of the compounds. The drug-likeness of the tested compounds was assessed as G-protein coupled receptor (GPCR) ligands, ion channel modulators (ICM), kinase inhibitors (KI), nuclear receptor ligands (NRL), protease inhibitors (PI), and enzyme inhibitors (EI).

Molecular dynamics simulations were carried out using Galaxy Europe to investigate the interactions between ligands 9, 11, 13, and ofloxacin as reference drugs and the chain of 3E0M. These simulations were conducted with GROMACS version 2021.6 (Van Der Spoel et al., 2005) and the AMBER99SB force field (Hansson et al., 2002; Showalter and Brüschweiler, 2007), offering a detailed account of atomic interactions. Molecular dynamics simulations are a robust computational technique, surpassing less computationally intensive docking methods. They provide high precision and offer insights into complex system behavior, often inaccessible through experimental means (Hansson et al., 2002; Showalter and Brüschweiler, 2007). In this study, we performed additional MD simulations for compounds 9, 11, 13, and ofloxacin with the A chain of 3E0M to validate previous findings. The protein topology parameters were generated using the Galaxy European Server (Bray et al., 2020), a widely used molecular modeling and simulation platform. The simulation parameters included the TIP3P water model, and hydrogen atoms were initially excluded from the GROMACS setup. Hydrogens were added to the ligand at pH 7.4, and during MD topology generation, the molecule’s charge was maintained at 0 with a multiplicity of 1. The gaff force field was applied for parameterization. The ligand and protein files were combined, and a structural configuration was established within a 1 nm triclinic box. The system was solved with SPC water molecules in a triclinic box, and sodium and chloride ions were introduced to attain standard salt concentrations and neutralize the system (Presti et al., 2016).

An equilibration process was conducted to ensure system stability, utilizing position-restrained dynamics (NVT) at 300 K for 3,000 ps, implemented with the leapfrog algorithm (Cuendet and van Gunsteren, 2007; Esha et al., 2023; Uddin et al., 2023). Following this equilibration phase, the system underwent a production run for an additional 3,000 ps under constant temperature and pressure conditions. Subsequently, the system was further simulated for 20 ns at the same temperature and pressure. Various GROMACS utility tools, including gmx rmsd, gmx gyrate, gmx rmsf, and ‘gmx hbond’, were employed to generate graphs for RMSD, radius of gyration, RMSF, and hydrogen bonds of the ligands at 300 K over the 20 ns period. Furthermore, principal component analysis (PCA) was utilized to assess the stability of protein complexes, proving highly advantageous (Biswal et al., 2023; Panda et al., 2023). PCA was carried out for the ligand-protein complexes (9, 11,13, and ofloxacin) at 300 K, as well as for compound 11 at temperature settings of 300 K, 305 K, 310 K, and 320 K, using the Bio3D package of GALAXY Europe server (Grant et al., 2006; Afgan et al., 2018; Kumar et al., 2022). The cosine content of the ligands was also evaluated. These MD simulations have provided valuable insights into the behavior of the protein-ligand complexes, aiding in a better comprehension of the underlying physical principles governing the structure and function of biological macromolecules.

Evaluating naphthoquinone and its derivatives (1–13) through fragment molecular analysis (FMO) has provided valuable insights into their chemical reactivity and stability. A broader energy gap (Egap) indicates heightened stability and reduced reactivity, whereas a narrower gap signifies softness, corresponding to increased reactivity and diminished stability. In addition to Egap, the energy levels of the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) play a significant role in understanding the electron-donating and -accepting characteristics of these compounds. To elaborate further, compounds with a narrow energy gap are classified as “soft” (σ), indicating high chemical reactivity and low stability, while those with a wider gap are labeled as “hard” (η). In this investigation, various molecular properties, such as ionization potential (IP), chemical potential (µ), electron binding energy (H), global chemical hardness (η), global chemical softness (σ), energy change (∆E), electrophilicity (ω), electronegativity (χ), electron affinity (EA), and dipole moment, were computed using the B3LYP/6-31G (d,p) method. The results of these calculations are presented in Table 1 and illustrated in Figure 3 and Supplementary Figures S1–S15.

Table 1 indicates that compounds 9, 10, 11, 12, 13, and D2 (ofloxacin) exhibited energy gaps (Egap) of 3.202 eV, 2.998 eV, 2.980 eV, 2.977 eV, 2.975 eV, and 4.369 eV, respectively. Among these, compounds 9 to 13 displayed low energy gaps, suggesting good reactivity and stability, requiring less energy for electron promotion from HOMO to LUMO orbital. Compared to D2, their energy gaps were lower, indicating greater stability and a willingness to participate in reactions. Compound 9 had the highest Egap of 3.202 eV, signifying lower reactivity but increased hardness and stability. Furthermore, compound 13 showed the lowest HOMO-LUMO gap (2.975 eV) with a dipole moment of 2.504 D, implying efficient biological activities. The energy gap values followed this order: compound 9 (3.202 eV) > compound 10 (2.988 eV) > compound 11 (2.98 eV) > compound 12 (2.977 eV) > compound 13 (2.975 eV). Compounds 9, 11, and 13 were chosen for further evaluation through additional computational analysis and molecular dynamics.

Figure 3 illustrates Molecular Electrostatic Potential (MEP) maps for compound 11. Electrophilic sites (depicted in blue) function as electron acceptors, while nucleophilic sites (shown in red) serve as electron donors. Partially nucleophilic regions are represented in yellow. The color-coded electrostatic potential highlights negative potential (in red, orange, and yellow shades) around electronegative atoms such as oxygen, indicating electrophilic reactivity. Conversely, positive potential (in blue) is observed over hydrogen atoms, signaling nucleophilic reactivity. Green areas denote the neutral potential. MEP maps provide insights into molecular stability by indicating consistent potential and are also helpful in studying molecular recognition and understanding interactions based on electrostatic complementarity.

Molecular docking is a vital computational approach for evaluating the binding affinity between a ligand and a protein or receptor. This method offers valuable insights into how compounds interact with a protein’s active site residues and influence cellular processes. In our study, we performed molecular docking experiments involving naphthoquinone derivatives (1–13) and three reference drugs—ciprofloxacin, ofloxacin, and vancomycin—against four specific proteins: MSR (3E0M), MERS-CoV (4RNA), EFGR (4UV7_L), and MDD (4DPT), as detailed in the table. To ensure the quality of the protein structure (PDB: 3E0M), we conducted various assessments, including PROCHECK, ERRAT, Verify3D, ProSA-web server, and ProSA-web energy server, as illustrated in Figure 4. These evaluations provided valuable insights into the protein’s structure. The structure met high-quality standards, as evidenced by the Ramachandran plot, overall quality factor, Verify3D score, and ProSA-web energy plot. In summary, the results confirm the validation and excellent quality of the protein structure, providing a solid foundation for our molecular docking Studies.

Table 2 illustrates the docking results, indicating that compounds 1–13 demonstrate notably higher binding affinities with the 3E0M protein compared to the other proteins in the study. A recent study reported that compound 11 exhibited a lower minimum inhibitory concentration in vitro. Our findings further highlight that compound 11 displays the highest binding affinity among all tested proteins, mainly showing strong affinities for Peptide methionine sulfoxide reductase msrA/msrB (PDB: 3E0M) −7.9 kcal mol⁻¹, aligning with our desired outcome. Compound 11 also exhibits binding affinity, Mevalonate diphosphate decarboxylase (PDB: 4DPT) −7 kcal mol⁻¹, The middle east respiratory syndrome coronavirus (MERS-CoV) papain like protease (PLpro) (PDB: 4RNA) −6.2 kcal mol⁻¹, and Epidermal Growth Factor Receptor (PDB: 4UV7) −6.5 kcal mol⁻¹. Table 2 in the SI shows that in the evaluation of MSR (PDB: 3E0M), all reference drugs demonstrated lower docking values ranging from −7.3 to −8 kcal mol⁻¹, in contrast to compounds 9 (−7.6 kcal mol⁻¹), 13 (−7.5 kcal mol⁻¹), and 11 (−7.9 kcal mol⁻¹). Additionally, the results indicate that most of the examined docking values for our compounds surpass those of the reference drugs used in antibacterial and cancer treatment. To summarize, compounds 11 and ofloxacin (−7.8 kcal mol⁻¹) exhibit strong binding affinities for all targeted proteins, while some compounds, such as 7 (−5.5 kcal mol⁻¹) and 8 (−5.8 kcal mol⁻¹), 12 (−6.2 kcal mol⁻¹), exhibit lower binding affinities, possibly attributed to differences in their chemical structures.

When treated with naphthoquinone derivatives, the functions of key enzymes like MERS-CoV papain-like protease (PLpro), Epidermal Growth Factor Receptor (EGFR), and Mevalonate Diphosphate Decarboxylase are potentially altered, impacting their roles in disease processes. Naphthoquinone derivatives, known for their broad biological activities, could inhibit MERS-CoV PLpro, disrupting the virus’s replication and ability to evade the immune system. This inhibition could be crucial in developing new antiviral strategies against MERS-CoV. In the case of EGFR, which plays a significant role in cancer cell proliferation, naphthoquinone derivatives might act as inhibitors, blocking EGFR’s signaling pathways that lead to tumor growth and metastasis. This suggests a promising avenue for cancer therapy, particularly in EGFR-driven cancers. Lastly, targeting Mevalonate Diphosphate Decarboxylase, a key enzyme in cholesterol biosynthesis, with these derivatives could reduce cholesterol synthesis. This action holds potential in treating disorders related to cholesterol metabolism, such as cardiovascular diseases. Thus, applying naphthoquinone derivatives to these enzymes offers a multifaceted therapeutic approach, spanning antiviral, anticancer, and cholesterol-lowering effects, demonstrating the versatility of these compounds in disease treatment and management.

In compound 11 and its interaction with the 3E0M protein, four distinct interactions have been identified, as outlined in Table 3 and Figures 5, 6 (Supplementary Figures S16–S18). These interactions include four hydrophobic interactions, specifically Pi-Alkyl interactions with MET A:116 and ARG A:73 on the residue and Pi-sigma interactions with LEU A:69. Additionally, one interaction is classified as a conventional hydrogen bond at the ALA A:165 residue. Compound 13’s protein-ligand interaction analysis with 3E0M revealed one conventional hydrogen interaction at the VAL A:465 residue, one halogen interaction with GLY A:367, and four Pi-Alkyl interactions involving residues ALA A:556, ARG A:415, ALA A:366, and VAL A:606. Compound 9 exhibited a binding affinity for 3E0M (−7.4 kcal/mol) with seven interactions in total, as shown in Figure 5 and Supplementary Figure S16. These interactions included two conventional hydrogen bonds (ARG A:261, ILE A:162), a Carbon Hydrogen Bond with ASP A:163, and Pi-Alkyl interactions involving MET A:116, ALA A:164, and ARG A:261.

The pivotal finding of this study highlights compound 11 as an exceptionally promising candidate for drug development due to its profound interaction with the 3E0M protein (Figure 7). Compound 11 demonstrates superior binding affinity to the 3E0M protein compared to other compounds and reference drugs like ofloxacin and ciprofloxacin, engaging four distinct amino acid residues. This comprehensive and specific binding profile underscores the remarkable affinity of compound 11 for the target protein. Based on the molecular docking and interaction analyses, compound 11 exhibits heightened binding affinity and forms a stronger interaction with the target protein compared to other compounds. This indicates its potential as a more effective inhibitor of proteins 3E0M and 4DPT, which are relevant in anti-cancer treatment and antibacterial activities. Consequently, this in silico study suggests that compound 11 holds substantial promise as an inhibitor of bacterial cell membrane proteins.

Two critical criteria, Lipinski and Veber, play a pivotal role in evaluating the physicochemical attributes of compounds. The investigation of these properties is essential to determine whether Naphthoquinone derivatives (1–13) conform to Lipinski and Veber guidelines. Lipinski’s criteria, crucial for preparing a compound for oral use, include five key conditions: a. Molecular weight (MW) of 500 g/mol or less b. Partition coefficient (log P) of 5 or less c. No more than 5 hydrogen bond donors (HBD) d. A maximum of 10 hydrogen bond acceptors (HBA) e. Topological polar surface area (TPSA) within 140 Å² (Veber et al., 2002; Lipinski, 2004). Veber’s guidelines add two more requirements for effective bioavailability: a. Fewer than 10 rotatable bonds (nrotb) b. TPSA not exceeding 140 Å², aligning with Lipinski’s rules. We utilized SwissADME to evaluate the compliance of compounds 1–13 with these criteria, particularly concerning their potential biological activity. Our analysis revealed that all compounds (1–13) adhered to Lipinski and Veber’s specific requirements. Furthermore, with a drug-like (bioavailability) score of 1, all compounds (1–13) were confirmed to meet the criteria for developing new medications, supported by strong theoretical evidence, with high agreement (less than 6) for drug compounds having a molecular weight below 500 g/mol, except vancomycin (MW = 1,449.2 g/mol), MLOGP of 4.15, and Log S (ESOL) values within the specified guidelines. Medicinal chemistry analysis, including the assessment of PAINS #alerts and Brenk#alerts, was also conducted using SwissADME for the studied structures. Among the 13 synthesized naphthoquinone derivatives (1–13), compound 11 exhibited the highest negative binding affinity, surpassing −7.9 kcal mol⁻¹, as shown in Table 2. Notably, compound 11 displayed favorable drug-likeness characteristics, with no PAINS #alerts and 0 Brenk #alert among the derivatives. Additionally, in silico ADMET profiling of compound 11, which demonstrated the highest negative binding affinity in molecular docking, met the requirements, indicating good drug-like properties (refer to Table 4). Consequently, based on their promising physicochemical properties, these compounds hold potential for further development as novel medications.

In our study, we utilized Molinspiration Chemoinformatics to assess the drug-like properties of naphthoquinone derivatives (1–13), targeting diverse biological receptors such as G protein-coupled receptors (GPCRs), ion channel modulators (ICMs), kinase inhibitors (KIs), nuclear receptor ligands (NRLs), protease inhibitors (PIs), and enzyme inhibitors (EIs), as detailed in Table 5. Notably, compounds 9 and 11 exhibited GPCR values of −0.19 and −0.16, respectively, while olioxacin demonstrated a GPCR value of 0.23. The lower GPCR value indicates better drug-like properties. Efficient development of novel drugs requires exploration of pharmacokinetic properties, encompassing absorption, distribution, metabolism, excretion, and toxicity (ADMET). In our investigation, we employed free tools such as SwissADME (http://www.swissadme.ch/index.php) (Daina et al., 2017) and AdmetSAR (http://lmmd.ecust.edu.cn/admetsar2/) (Yang et al., 2019) to evaluate the ADMET properties of all 13 active compounds (1–13). We assessed seven essential ADMET characteristics, including cytochrome P450 enzymes (CYP3A4 and CYP2C19) inhibition, hERG inhibition, blood-brain barrier (BBB) penetration, human intestinal absorption (HIA), and synthetic accessibility (SA) score, presented in Table 6. It is worth noting that BBB penetration is particularly crucial for drugs targeting the central nervous system (CNS) but may not be as critical for drugs with minimal CNS impact. BBB penetration is categorized into three ranges: high (>2), medium (2–0.1), and low (<0.1) adsorption (Veber et al., 2002). The results indicate that all compounds (1–13) meet ADMET standards for drug-likeness (bioavailability) and hold the potential for developing novel drugs.

Our findings reveal that, with the exceptions of compounds 8 and 10, most compounds exhibited moderate blood-brain barrier (BBB) penetration values, falling within the range of 0.625–1. Ciprofloxacin and ofloxacin displayed medium absorption values of 0.714 and 1, respectively, while vancomycin exhibited low absorption at −0.914. Moreover, the majority of naphthoquinone derivatives demonstrated favorable synthetic accessibility, ranging from 1 to 3.68, comparable to standard drugs. Compound 11, with a synthetic accessibility score of 3.97, and ofloxacin, with a score of 2.44, exemplified these characteristics. Our predictions regarding the potential risk of hERG activity inhibitors, with values ranging from −0.4 to −0.82 (Shadrack and Ndesendo, 2017), indicated that compound 11 achieved the desired hERG pIC50 value of −0.6441. In contrast, all other compounds remained within the reference range. These results affirm that all compounds (1–13) adhere to ADMET standards, signifying drug-likeness (bioavailability), and underscore their potential for novel drug development. Specifically, compound 11 exhibits favorable properties across all ADMET parameters, including moderate BBB penetration, high intestinal absorption, and a low risk of hERG inhibition. These characteristics align with the desired profile for a potential therapeutic agent.

We conducted molecular dynamics (MD) simulations spanning 20 ns for compounds 9, 11, 13, and ofloxacin in conjunction with the MSR protein (PDB: 3E0M). Previous molecular docking investigations revealed robust binding affinities of −7.6, −7.9, −7.5, and −7.8 kcal/mol for compounds 9, 11, 13, and ofloxacin, respectively, with the 3E0M protein, indicating their potential as promising drug candidates. The MD simulations aimed to assess the stability and interactions within the protein-ligand complex throughout the 20 ns trajectory. Various parameters, including RMSD (root-mean-square deviation), RMSF (root-mean-square fluctuation), Rg (radius of gyration) values, potential energies, temperature, hydrogen bonding, and principal component analysis (PCA), were employed for the simulation analysis. Additionally, the complex underwent simulations at four distinct temperatures (300, 305, 310, and 320 K) to explore its structural alterations under varying thermal conditions. This approach provided insights into the stability and adaptability of the protein-ligand complex under different physiological conditions.

The analysis of RMSD values obtained from the MD simulations provides insights into the stability and conformational dynamics of ligand-protein complexes (refer to Figure 8A and Supplementary Figures S19–S21). Notably, the complex involving compound 11 and protein 3E0M displayed the highest conformational stability, with RMSD values ranging from 0.2 to 0.55 nm (blue line), indicating minimal deviation from the initial structure. The complex of compound 13 and 3E0M exhibited the second-highest stability, showing RMSD values from 0.3 to 0.6 nm, although it demonstrated larger fluctuations, reaching up to 0.7 nm. In contrast, the complexes involving ofloxacin and 3E0M, as well as compound 9 with 3E0M, displayed lower stability, evident from increased scattering in their RMSD plots, with values spanning from 0.2 to 0.7 nm and 0.1 to 0.8 nm, respectively.

In addition to RMSD, root-mean-square fluctuation (RMSF) assessment was utilized to gauge the flexibility of individual amino acid residues within the protein-ligand complexes throughout the MD simulations. RMSF values offer insights into the dynamic behavior of specific amino acids, indicating their level of flexibility or rigidity. Analyzing the RMSF values for each compound interacting with protein 3E0M revealed variations in residual flexibility among the complexes. The observed fluctuations ranged approximately from 0.04 to 0.13 nm for ligand 11, 0.03–0.9 nm for the target protein 3E0M, and 0.05–0.9 nm for the protein-ligand complex (see Figure 8B and Supplementary Figures S22–S24). These values suggest that a majority of amino acid residues in the complex exhibited fluctuations. Notably, regions directly interacting with ligand 11 exhibited low fluctuations, while other regions containing amino acids without direct interactions displayed more significant fluctuations. The complex involving protein 3E0M and ligand 11 demonstrated the highest RMSF values, indicating greater flexibility compared to other complexes (refer to Figure 8B and Supplementary Figures S22–S24). This implies that ligand 11 may induce a more dynamic conformation in the protein compared to the other ligands. The complex of protein 3E0M with ligand 13 also displayed moderate flexibility, while the complexes of protein 3E0M with ligand 9 and ofloxacin-3E0M exhibited relatively higher rigidity.

The variations in RMSF values suggest diverse flexibility within the protein-ligand complexes, potentially influencing their biological functions. The analysis of RMSF provides valuable insights into the structural dynamics of protein-ligand interactions, complementing the information derived from RMSD analysis. These findings indicate the stability of both the ligand and the protein throughout the simulation, with minor deviations in atomic positions. Notably, the observed fluctuations in specific ligand atoms may carry significance, potentially related to particular functional or interactional aspects of the protein. However, further investigation is necessary to ascertain these fluctuations’ precise implications and potential impacts on the biological activities associated with 3E0M. In this computational study, it has been assessed that the radius of gyration (Rg) is a parameter used to characterize the size and shape of macromolecules like protein-ligand complexes. The Rg value determined for the rotein-ligand complex ranged from 2.12 to 2.45 nm, as depicted in Figure 9A (see Supplementary Figures S25–S28).

Notably, all findings indicate that including the ligand did not induce substantial alterations in the overall size or shape of the protein. This outcome aligns with the notion that the binding of the ligand to the protein occurred in a manner that preserved the protein’s tertiary structure without causing significant disruptions. Furthermore, the comparison in Rg values between the protein and the protein-ligand complex suggests that the binding event did not bring about notable structural modifications within the protein. Nevertheless, it is crucial to acknowledge that while Rg values provide insights into a protein’s general size and shape, they may not unveil the precise molecular-level structural alterations that have transpired. Further experiments, such as X-ray crystallography or NMR spectroscopy, may be indispensable for a more comprehensive understanding of the molecular interactions between the protein and the ligand. The results of the MD simulations emphasize the significance of the number of hydrogen bonds (HB) formed between the protein and the ligand—a critical factor in sustaining the stability of the protein-ligand system. These findings are depicted in the time-dependent assessment of intermolecular hydrogen bonds during the 20 ns simulation, as shown in Figure 9B and Supplementary Figure S29. The simulation conducted at the standard body temperature exhibited the most consistent HB profile, with the number of hydrogen bonds fluctuating between zero and four throughout the entire 20 ns simulation period. Notably, between the 5 ns and 15 ns timeframe (Figure 9B), ligand 11 demonstrated the highest occurrence of hydrogen bond formations with the residues of the target protein (PDB: 3E0M). Initially, the analysis of interactions within the docked complex did not clearly indicate the presence of hydrogen bonds, as seen in Figure 8B. However, the subsequent 20 ns MD simulation revealed the emergence of hydrogen bond interactions for ligands 9, 11, 13, and the reference ligand (ofloxacin).

Upon comparing the outcomes of simulations involving compounds 9, 11, 13, and ofloxacin with the MSR protein (PDB: 3E0M), it is evident that compound 11 establishes a more robust interaction with 3E0M. This suggests its potential as an inhibitor of the MSR protein, which holds significance due to its crucial role in cancer and bacterial inflammation. Previous in vitro and in vivo experiments have underscored the promise of compound 11 as a candidate for an anti-cancer and antibacterial drug. This is attributed to its potent anti-mycobacterial activity, low cytotoxicity, and unique lipophilic vehicle feature. The complex comprising compound 9 and MSR showed HOMO-LUMO combined results slightly lower than ofloxacin, which exhibited the most pronounced negative binding affinity and met all criteria for consideration as a drug candidate, as substantiated by the molecular simulation results at 310 K. To further corroborate our identification of compound 11 as a promising anticancer and antibacterial drug candidate, additional molecular dynamics simulations were performed at three different temperatures: 300, 305, 310, and 320 K for the complex of ligand 11 with 3E0M, as depicted in Figure 10A.

The RMSD analysis of the ligand 11–MSR protein complex revealed its optimum stability at 310 K, with minimal deviation and a peak value of 0.8 nm. At 305 K, the complex exhibited the second-highest stability, with RMSD values fluctuating between 0.1 and 0.4 nm and peaking at 0.5 nm. In contrast, the complex demonstrated reduced stability at 320 K and 300 K, manifesting more pronounced deviations. At these temperatures, the RMSD values ranged from 0.2 to 0.4 nm and 0.1 to 0.7 nm, respectively (see Figure 10B and Supplementary Figures S30–S44).

Furthermore, the system exhibited sustained stability at different temperatures (300, 305, 310, and 320 K) over the 20 ns simulation, as depicted in Figure 10A. Additionally, conformational principal component analysis (PCA) based on Cα atoms was applied to the molecular dynamics (MD) trajectories of the protein (PDB: 3E0M) and ligand 11 at various temperatures (refer to Figure 11 and Supplementary Figures S45–S48). This study utilized PCA to evaluate the protein’s variance, collective motions, and conformational changes during the MD simulations, as shown in Figure 11 (Table 7 and Supplementary Figure S17). The PCA analysis was performed on the MD trajectories of the target protein (ID: 3E0M) and the ligand 11 complex at 300, 305, 310, and 320 K, using the Bio3D package. The resulting plots illustrated eigenvalues versus eigenvectors, highlighting the primary motions extracted from the trajectories, with a focus on the first three eigenvectors (PC1, PC2, and PC3) and color dots representing the variance captured by the eigenvectors.

Table 7 provides a concise overview of the primary motions observed in the (ligand 11-3E0M) complex at various temperatures, specifically focusing on a subset and comparing the first three eigenvectors (PC1, PC2, and PC3). Notably, the complex simulated at 310 K (physiological body temperature) exhibited the most significant variation in PC1 (24.39) and the least in PC2 (19.19) and PC3 (14.14) concerning internal motion within the MD trajectory. This implies a robust interaction between compound 11 and 3E0M, a finding supported by both in vivo and in vitro investigations. However, further clinical applications are necessary to validate the efficacy and safety of compound 11 as an MSR-3E0M inhibitor.

In previous in vitro investigations, ligand 11 demonstrated significant potency and specificity as an inhibitor of biofilm formation and various bacterial diseases. Moreover, this ligand shows promise as a potential cancer treatment. Docking results revealed that ligand 11 exhibited the highest binding affinity (−7.9 kcal mol⁻¹) with MSR-3E0M, while its minimum inhibitory concentration (MIC) values against MSSA, MSRA-ATCCBAA1717, and MSRA-ATCCBAA-44 (Table 8) ranged from 1.25 to 2.5 (Song et al., 2020).

In comparison, ofloxacin, with a binding affinity of −7.3 kcal mol⁻¹, displayed MIC values ranging from 0.25 to 16. Ofloxacin is a promising candidate for treating cancer and bacterial diseases, given its efficacy and selectivity against methicillin-resistant staphylococcus aureus (MRSA). It targets the enzyme responsible for biofilm formation, causes cell membrane damage, chelates intracellular iron ions, and generates intracellular reactive oxygen species. Ligand 11 also demonstrated stable binding to methicillin-resistant staphylococcus aureus, with consistent RMSD and RMSF results. Additionally, it maintained a stable conformation at human body temperature. Ligands 9 and 13 also exhibited potential as inhibitors of bacterial diseases, boasting lower Egap values and higher docking affinity than ofloxacin. These ligands formed interactions with amino acid residues that had positive results in previous studies. Particularly, ligand 11 showed the highest binding affinity for the A chain of MSR-3E0M.

To validate our findings and determine whether 2-(4-butylbenzyl)-3-hydroxynaphthalene-1,4-dione (11) is optimal, we extended our investigation to include two additional naphthoquinone derivatives (9 and 13) and the standard drug, ofloxacin, which exhibited the most favorable docking results. Subsequently, we conducted 200 ns MD simulation revealed molecular dynamics simulations, and the outcomes are presented in Figure 12, comprising three RMSD panels: a) Merge protein-ligand complex of ligands 9, 11, 13 and ofloxacin, b) protein-ligand complex 11 and c) Protein-ligand complex of ofloxacin. Figure 12 shows the root mean square deviation (RMSD) analysis for the three naphthoquinones and the standard drug. In Figure 12A) the yellow line corresponds to ligand 11, while the grey, red, and purple lines represent ligands 13, 9, and the standard drug, ofloxacin, respectively. Naphthoquinone derivative 11 is the only one exhibiting an RMSD value close to 0.35 nm, whereas the others collectively yield a value over 0.60 nm except the reference drug ofloxacin (0.40 nm). This notably higher combined value indicates that these ligands possess dissimilar shape orientations and exhibit less structural similarity. Notably, the naphthoquinone derivative (11) displays a lower RMSD value than the other compounds, as evident in Figures 12A, B. This underscores its unique structural characteristics and suggests its potential as the most promising candidate. For more detailed information, including individual graphs for each naphthoquinone with respect to the reference molecule (PDB: 3E0M), please refer to the supplementary information (Supplementary Figure S45).

Figure 13 illustrates the merged radius of gyration for three naphthoquinones and the standard drug. In Figure 13A, we observe a radius of gyration range of 0.32–0.35 nm for naphthoquinone derivative 11, 0.37–0.41 nm for the standard drug (ofloxacin), and 0.38–0.42 nm for ligand 9 and 0.5–0.6 nm for ligand 13. It is important to note that higher values in the radius of gyration indicate lower structural similarity and fewer interactions among molecules. Remarkably, the ligand 11 exhibits the highest structural similarity and interactions among its molecules due to its lower radius of gyration. Detailed radius of gyration values and graphs for each naphthoquinone derivative can be found in the supplementary information (see Supplementary Figure S47). In conclusion, simulation results for ligand 11 provided valuable insights, indicating a robust interaction with 3E0M. These findings suggest favorable interactions with the protein, supported by both computational and laboratory experiments. Nevertheless, it is crucial to emphasize that further in vivo studies are necessary to validate the efficacy and safety of ligand 11 as an inhibitor of 3E0M. In conclusion, simulation results for ligand 11 provided valuable insights, indicating a robust interaction with 3E0M. These findings suggest favorable interactions with the protein, supported by both computational and laboratory experiments. Nevertheless, it is crucial to emphasize that further in vivo studies are necessary to validate the efficacy and safety of ligand 11 as an inhibitor of 3E0M.