The sequence of the GroEL protein (P61439), which comprises a total of 546 amino acids, was retrieved from the UniProt database in FASTA format (The UniProt Consortium, 2017). Subsequently, the 3D structure of GroEL was predicted using the AlphaFold program (Jumper et al., 2021) (Figure 1). To identify potential active site residues in GroEL, six complementary structure-based prediction tools were employed: COACH (Yang et al., 2013), TM-SITE (Yang et al., 2013), S-SITE, COFACTOR (Roy et al., 2012), FINDSITE (Brylinski and Skolnick, 2008), ConCavity (Capra et al., 2009 ), and CASTP (Tian et al., 2018). The outputs of these tools were used to generate a consensus list of residues for subsequent virtual screening experiments.

Virtual screening is a computational method that uses high-performance computing systems to identify, filter, and evaluate potential molecular conformations of chemical compounds from extensive databases. The structure of GroEL was prepared for virtual screening using the Protein Preparation Wizard in the Schrödinger Suite (Schrödinger, LLC, New York, NY, 2017-1) (Sahoo et al., 2024a; Sethi et al., 2024). The preparation steps included reconstructing the missing side chains and loops, assigning bond orders, adding hydrogen atoms at neutral pH, and removing non-essential water molecules. The structure was then subjected to energy minimization to resolve steric clashes and optimize the geometry. A compound library of 543,503 small molecules was obtained from the Life Chemicals Database. The ligands were processed using the LigPrep module of the Schrödinger Suite to generate low-energy conformations, assign ionization states, and optimize geometries. Docking grids were created around the predicted active-site residues of GroEL. SBVS was conducted using the GLIDE module. Three levels of precision were employed sequentially: High-Throughput Virtual Screening (HTVS) for rapid initial screening, Standard Precision (SP) for refinement, and Extra Precision (XP) for detailed ligand–receptor interaction assessment. The top-ranking compounds from XP docking were retained for further analysis. Molecular interactions and binding poses were visualized using UCSF Chimera (Pettersen et al., 2004) and the GLIDE 2D interaction diagram tool. Figure 2 presents the overall computational workflow adopted for the structure-based virtual screening of GroEL.

The absorption, distribution, metabolism, and excretion (ADME) properties of all selected chemical scaffolds were evaluated using the QikProp module from the Schrödinger Suite (Schrödinger, LLC, New York, NY, 2017-1). This extensive assessment provides valuable insights into their potential as drug candidates by examining essential pharmacokinetic parameters. Key molecular descriptors were considered, including compliance with the Rule of Five (Ro5) (Lipinski et al., 1997). Together, these evaluations provided an integrated understanding of both the pharmacokinetic behavior and potential safety liabilities, thereby supporting a more reliable assessment of the drug-like properties of the compounds for subsequent studies.

Quantum mechanical calculations were performed using DFT to investigate the electronic characteristics and geometric configurations of the target compounds using the Jaguar module of the Schrödinger Suite (version 8.7) (Bochevarov et al., 2013). This computational approach provides an efficient and accurate framework for determining molecular electronic structures, optimized geometries, and physicochemical properties at the ground state level. The calculations focused on identifying critical electronic descriptors, including the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), HOMO-LUMO energy gap (ΔE), and molecular electrostatic potential (MEP) surfaces. All of the computations employed the B3LYP hybrid density functional in conjunction with the 6-31++G(d,p) basis set (Lee et al., 1988; Bhatta et al., 2015). These comprehensive analyses provided valuable insights into the electronic reactivity, charge distribution, and potential binding characteristics of the designed molecules.

MDS were conducted to examine the structural stability and dynamic behavior of the top-ranked protein-ligand complexes. Simulations were performed using the GROMACS 2022.3 software package with the CHARMM36 force field, which is widely recognized for its reliability in all-atom protein-ligand modeling (Bekker et al., 1993; Van Der Spoel et al., 2005). Ligand topologies were generated using the CGenFF server (Vanommeslaeghe et al., 2010) and integrated with the protein topology files prepared using the CHARMM36 force field (Bjelkmar et al., 2010). Each protein-ligand complex was positioned within a cubic simulation box, maintaining a 1nm buffer between the solute and box boundaries. The system was solvated using the TIP3P water model to simulate a realistic aqueous environment (Jorgensen et al., 1983). A two-step energy minimization process was performed to remove steric clashes and achieve an energetically favorable conformation. The steepest descent algorithm was first applied, followed by the conjugate gradient method, with 50,000 minimization steps. Following energy minimization, the system was equilibrated under both the NVT and NPT ensembles for 100 picoseconds (ps) each. Temperature was controlled using the Berendsen thermostat (Berendsen et al., 1984), and pressure was regulated using the Parrinello–Rahman barostat (Parrinello and Rahman, 1981). After the equilibration phase, MDS were executed for 100 ns for each protein-ligand complex. The coordinates of each complex were recorded at regular intervals of 2 fs. The final resulting trajectory files were analyzed using built-in GROMACS tools to assess the stability of the GroEL-ligand complexes. Parameters such as the root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), hydrogen bonding interactions, and solvent-accessible surface area (SASA) were examined to determine the dynamic behavior and conformational stability of the complexes throughout the simulation, as described in our previous studies (Sahoo et al., 2024b, 2025).

PCA, a robust method for multivariate statistical analysis, was employed to identify and characterize the primary motions within the protein-ligand complexes. This approach involves calculating the dominant motions using eigenvalues and eigenvectors, as previously described (Sahoo et al., 2023; Samantaray et al., 2025). The analysis used the final 50 ns of the MDS trajectories to construct a covariance matrix, focusing on the backbone atoms of the protein in each complex. The cosine content for the leading eigenvectors was computed to evaluate simulation convergence, and eigenvalues were subsequently obtained (Amadei et al., 1993). FEL analysis was performed to visualize the energy minima and their corresponding conformational states, allowing the identification of stable and metastable regions within the protein–ligand complexes (Maisuradze and Leitner, 2007). This approach provides insight into the system’s response to ligand binding by elucidating the underlying energetics and quantifying the overall thermodynamic stability of the complexes. The FEL was constructed using the first two principal components (PC1 and PC2), obtained from PCA, using the g_sham module in GROMACS (Lindahl et al., 2001). The resulting FEL plots were generated and visualized using Mathematica software.

The Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) method was employed to calculate the protein-ligand complex binding free energies (BFE) (Genheden and Ryde, 2015). The MM-PBSA is a widely used computational approach in drug design that combines molecular mechanics force fields, continuum solvent models, and solvation energy terms to estimate the thermodynamic properties of protein-ligand interactions. The BFE is derived from three main energy components: non-polar solvation energy, polar solvation energy, and vacuum potential energy. To assess the stability of the complexes, MDS were run for last 50 ns of the trajectory information, and the g_mmpbsa tool was used with default settings and a solute dielectric constant of 2.0 (Kumari et al., 2014). This method allows the decomposition of the binding free energy into various contributions, such as van der Waals interactions, electrostatic forces, solvation energy, and entropy, aiding in understanding the individual interactions that influence binding affinity. MM-PBSA calculations are critical for evaluating ligand binding affinities, comparing ligand poses, and optimizing drug candidates, and can be validated by experimental binding data to ensure the reliability of the computational predictions.

Structure-based virtual screening (SBVS) remains a pivotal strategy in contemporary drug discovery, particularly for identifying novel inhibitors of key protein targets involved in pathogenic processes, such as biofilm formation (Oselusi et al., 2024). In the present study, SBVS was employed to explore potential inhibitors targeting the GroEL protein of Leptospira interrogans, a molecular chaperone known to be involved in biofilm development and host-pathogen interactions. Active site prediction tools collectively identified 26 candidate residues, including T29, L30, G31, P32, D86, G87, T88, T89, T90, S149, G413, and G414. Several residues were consistently predicted by all six methods, reinforcing their reliability. A comparison with previously reported conserved regions revealed a strong overlap. Specifically, residues T29–P32 were located near the apical domain (191–202), suggesting a role in substrate recognition; residues D86–S149 corresponded to the intermediate domain (residues 362–381), associated with conformational flexibility, and residues G413–G414 were mapped to the equatorial domain (401–415), a region that is essential for ATP hydrolysis and chaperone activity (Ho et al., 2021). The residual evolutionary conservation profile of GroEL, highlighting the active site residues and conserved variants considered for MD simulation, is illustrated in Figure 2A and Figure 2B. A compound library consisting of 543,503 small molecules from the Life Chemicals database was screened against the GroEL active site. The screening protocol utilized hierarchical docking within the GLIDE module, progressively filtering candidates using HTVS, SP, and XP protocols. This workflow enabled the systematic enrichment of compounds with high binding potential at the predicted active site of GroEL. The multi-step screening pipeline successfully reduced the compound library to the top 10% of molecules with the most favorable binding scores. Structural visualization highlighted key binding residues within the predicted GroEL active site. In particular, residues G413 and G414, which were identified by all prediction tools and located within a conserved functional domain, emerged as potential anchoring sites for ligand binding.

Following the XP docking stage, five lead compounds F3385-2019, F1243-0200, F3139-0927, F2801-0179, and F1864–0208 were shortlisted based on their favorable docking scores (all less than -8.0 kcal/mol), indicating strong binding potential and high affinity for the GroEL binding site (Table 1). Among these, F3385–2019 exhibited the most favorable binding characteristics, achieving a docking score of −10.343 kcal/mol. F3385–2019 formed three hydrogen bonds with the key residues Thr90, Asn152, and Asp397; and establishes hydrophobic contacts with Leu30, Gly31, Pro32, Lys50, Asp51, Thr89, Gly87, Ile149, and Ser150, indicating a robust and well-anchored interaction within the active site of the GroEL protein. Compound F1243–0200 closely followed a docking score of −9.668 kcal/mol, stabilized by four hydrogen bonds involving Asp86, Gly87, Thr90, and Asn152. Its hydrophobic interaction profile was extensive, involving multiple conserved residues, including Gly31, Lys50, Leu30, Pro32, Asp493, Val491, and others, suggesting deep and stable accommodation within the binding pocket. F3139–0927 achieved a docking score of −8.445 kcal/mol, forming hydrogen bonds with Gly31, Gly87, Thr88, and Ile149, and hydrophobic interactions with Arg394, Lys50, Asp86, and Leu30, indicating moderate but potentially meaningful binding. Similarly, F2801–0179 demonstrated a docking score of −8.364 kcal/mol and was stabilized through hydrogen bonding with Asn152 and Ala478, supported by hydrophobic contacts involving residues like Thr29, Gly87, Ile149, and Val491, contributing to favorable ligand positioning and binding strength. Finally, F1864–0208 displayed a docking score of −8.257 kcal/mol. Despite being the lowest among the top five, it formed two hydrogen bonds with Thr89 and Thr90 and engaged in extensive hydrophobic interactions with key residues, including Asp493, Ser150, Asp86, Gly414, and Asn477, thereby retaining its binding favorability. Overall, detailed molecular interaction analyses suggested that these five lead compounds effectively interacted with GroEL, with F3385–2019 being the top candidate. These findings reinforce the potential of these small molecules as inhibitors of GroEL-mediated biofilm formation by Leptospira interrogans. The two-dimensional chemical architectures of the identified scaffold compounds are shown in Figure 3. The binding interactions were systematically characterized, with particular emphasis on hydrogen bond networks, hydrophobic contacts, and π-π stacking arrangements that contribute to complex stabilization (Figure 4).

Evaluation of ADME properties is essential for identifying potential drug candidates with optimal pharmacokinetic and safety profiles. In this study, the QikProp module of the Schrödinger Suite was used to predict the key ADME descriptors and drug-likeness characteristics of the top-ranking hit compounds (Table 2). These computational predictions provide preliminary insight into the oral bioavailability, permeability, and systemic exposure potential of the selected molecules. All five compounds were assessed for compliance with Ro5 criteria, which is a benchmark for drug-likeness prediction based on molecular weight, lipophilicity, and hydrogen bonding potential. Notably, none of the compounds violated the Ro5 criteria, indicating favorable physicochemical properties and supporting their potential as orally available drug-like molecules. The predicted human oral absorption (HOA) values were above 80% for all compounds except F3139-0927, which showed a comparatively lower value of 52.24%. Particularly, F2801–0179 demonstrated the highest HOA (96.10%), suggesting excellent oral absorption potential. The remaining compounds, including F3385–2019 and F1243-0200, also exhibited high HOA values of 90.91% and 82.31%, respectively. Membrane permeability predictions using MDCK (QPPMDCK) and Caco-2 (QPPCaCo) models indicated strong cellular permeability for most compounds. F2801–0179 exhibited the highest permeability, with QPPMDCK and QPPCaCo values of 300.27 and 573.84, respectively, followed by F3385-2019 (QPPMDCK: 214.34, QPPCaCo: 461.25), highlighting their efficient potential for intestinal and blood–brain barrier passage. The logP (QPlogPo/w) values ranged between −0.559 and 3.481, indicating acceptable lipophilicity for all compounds. Interestingly, F3139–0927 had a negative logP, suggesting higher hydrophilicity, which may explain its lower HOA and permeability values. The hydrogen bond donor and acceptor counts were also within favorable limits, ranging from 1–4 and 5.25–7.65, respectively. These values suggest that balanced polar surface characteristics are essential for bioavailability and target engagement. Collectively, the ADME predictions demonstrate that F2801–0179 and F3385–2019 possessed pharmacokinetic and drug-likeness profiles, combining high oral absorption, strong membrane permeability, and no violations of drug-likeness rules. These findings reinforce the need for further preclinical studies.

DFT calculations were performed using the Jaguar module of the Schrödinger Suite to assess the electronic characteristics and possible reactivities of the identified hit compounds. The HOMO and LUMO energies were determined, and the HOMO–LUMO energy gap (ΔE) was analyzed to evaluate the electronic stability and chemical reactivity of each compound (Table 3). Among the top candidates, F2801–0179 exhibited the narrowest HOMO–LUMO gap of −0.186 eV, suggesting high chemical reactivity and enhanced potential for electron transfer interactions with the GroEL protein. This was closely followed by F1243-0200, which showed a ΔE of −0.173 eV, and F3139–0927 with −0.164 eV, both indicative of favorable reactivity profiles. F3385–2019 showed a moderately wider gap of −0.141 eV, implying slightly lower reactivity but still within a range conducive to bioactivity. Interestingly, F1864–0208 displayed an anomalously large gap of −1.9 eV (Figure 5), which may either reflect a computational artifact or suggest a very low electron transfer potential, possibly correlating with its relatively lower docking scores. FMO analysis highlighted the electronic diversity among the top-ranked compounds and supported the findings of molecular docking analysis. Molecules with smaller energy gaps, particularly F2801-0179, are more likely to undergo charge transfer interactions with amino acid residues in the protein active site, thus enhancing binding strength and biological efficacy (Figure 5).

To reinforce the reliability of virtual screening outputs, MDS, PCA, and MM-PBSA free energy calculations have proven essential for validating structural stability and refining predicted binding affinities. Multiple studies demonstrate the utility of this approach; for example, MDS and MM-PBSA have been applied to prioritize repositioned inhibitors targeting PfEMP1 in Plasmodium falciparum, while docking, MDS, and experimental assays have been used to validate repurposed FDA-approved drugs against FZD10 in nasopharyngeal carcinoma (Ngernsombat et al., 2024; Verma et al., 2025). These findings underscore the importance of MD-based refinement as a critical step following virtual screening. To further elucidate the binding mechanism and stability of the selected natural ligands within the GroEL binding pocket, MD simulations were performed for 100 ns. The trajectories provided insights into the conformational behavior and ligand–protein interactions over the simulation timescale. Various structural and energetic parameters, including the RMSD, RMSF, Rg, hydrogen bond formation, SASA, and ligand RMSD, were assessed. The backbone RMSD trajectories (Figure 6A) provide a general view of the overall structural deviation and equilibration of GroEL upon ligand binding. Most ligand–protein complexes exhibited RMSD values stabilizing within 0.3 to 0.6 nm, indicative of good structural convergence. Notably, the F2801–0179 complex displayed higher backbone fluctuation during the final 30 ns of the simulation, with the RMSD peaking at approximately 1.7 nm, suggesting possible local conformational adjustments. In contrast, the F1243–0200 and F1864–0208 complexes exhibited relatively low and stable RMSD values throughout the simulation, maintaining fluctuations within 0.3 to 0.4 nm, indicating high structural integrity. F3385–2019 and F3139–0927 exhibited intermediate profiles, with moderate deviations (0.4–0.6 nm), suggesting favorable stability under dynamic conditions.

The Rg values, which reflect the overall compactness of the protein structure, remained consistent across all complexes (Figure 6B). The average Rg values ranged from 2.6 to 3.0 nm, indicating no significant unfolding or structural destabilization. However, the F2801–0179 complex showed a mild upward trend in Rg values after ~80 ns, correlating with its higher RMSD values and suggesting a slight expansion in structural volume. The remaining ligands, particularly F1243–0200 and F1864-0208, maintained steady Rg profiles, reinforcing the conclusion of a stable tertiary structure. RMSF analysis (Figure 6C) provided insight into residue-level fluctuations and local flexibility within the GroEL structure. The F2801–0179 complex exhibited elevated fluctuations across flexible loop regions, especially between residues 200–350, where peaks exceeded 1.2 nm, suggesting enhanced motion or flexibility in peripheral regions. In contrast, the F1243-0200, F1864-0208, and F3139–0927 complexes showed lower fluctuations (≤0.5 nm), particularly in regions critical for ligand binding, suggesting that these natural ligands effectively stabilize key domains of the protein. Hydrogen bonds play a crucial role in the stability and specificity of protein–ligand interactions. Figure 6D shows the number of hydrogen bonds formed between GroEL and each ligand over 100 ns. F1243–0200 consistently formed the highest number of hydrogen bonds and maintained up to six concurrent interactions throughout the simulation. This sustained H-bonding profile supports its strong affinity and anchoring capability. F1864–0208 also showed a stable H-bond pattern (2–4 bonds), whereas F2801–0179 and F3139–0927 exhibited more transient bonding, with counts fluctuating between 1 and 3. The F3385–2019 complex formed fewer but consistent hydrogen bonds (1–2), which may indicate hydrophobic contributions that supplement the interactions.

The SASA (Figure 7A) showed that all complexes maintained average values between 250 and 270 nm², indicating stable protein surface exposure. F2801–0179 consistently exhibited slightly elevated SASA values, correlating with its increased backbone RMSD and Rg, potentially because of a looser fit or partial exposure of surface residues. In contrast, F1243-0200, F1864-0208, and F3139–0927 displayed more compact profiles with fewer fluctuations, suggesting tighter ligand-induced packing and reduced solvent exposure. Ligand RMSD analysis (Figure 7B) was used to assess how well the ligand remained within the protein binding pocket during the simulation. F1243–0200 showed the highest ligand RMSD, reaching ~0.28 nm, suggesting minor reorientation within the pocket, though without dissociation. In contrast, F2801-0179, F3385-2019, and F1864–0208 displayed ligand RMSD values under 0.15 nm, indicating a firm and consistent binding pose. F3139–0927 remained tightly bound as well, reflecting excellent spatial retention.

Collectively, the MD simulation results suggest that F1243–0200 and F1864–0208 form the most stable and well-integrated complexes with GroEL, as evidenced by their low protein RMSD, tight ligand RMSD, stable hydrogen bonding, and compact Rg/SASA profiles. Although F2801–0179 demonstrated increased flexibility and solvent exposure, its stable ligand orientation and moderate H-bonding suggest potential as a flexible binder. These insights complement the docking analysis, highlighting F1243–0200 and F1864–0208 as promising natural inhibitors of GroEL. Notably, the distinct dynamic behavior observed for F2801–0179 across multiple structural metrics suggests a fundamentally different interaction mode compared with the other compounds. The elevated backbone RMSD, increased residue-level fluctuations in flexible loop regions, and higher SASA values indicate that F2801–0179 induces greater conformational adaptability within the GroEL binding pocket. This behavior may arise from differences in its chemical scaffold or suboptimal complementarity with key binding-site residues, resulting in weaker structural anchoring despite maintaining a stable ligand orientation. These observations highlight that increased binding flexibility does not necessarily translate into enhanced complex stability, underscoring the importance of integrating dynamic analyses beyond docking scores alone.

PCA was performed to investigate the collective atomic motions and conformational transitions of the GroEL protein in complex with each of the five selected natural ligands over the 50 ns of the MDS. PCA reduces the dimensionality of atomic fluctuations and reveals dominant movement patterns by analyzing the eigenvectors derived from the covariance matrix of Cα atomic displacements. The eigenvalue distribution (Figure 8) shows that the first three eigenvectors accounted for the majority of the conformational variance in each complex. F2801–0179 exhibited the highest magnitude of motion among the first few principal components, indicating a greater degree of flexibility in the complex. In contrast, F1243–0200 and F1864–0208 displayed lower eigenvalue contributions, suggesting reduced internal motion and higher conformational stability during the simulation. The PC1 versus PC2 projection highlights the spatial distribution of the conformations sampled by each system. The F1243–0200 and F1864–0208 complexes formed tightly clustered groups in essential space, reflecting limited conformational drift and a more compact motion profile. In contrast, F2801–0179 exhibited an extended, scattered trajectory, indicating broader sampling of the conformational space. F3385–2019 and F3139–0927 showed moderate clustering, with some degree of structural flexibility. The individual PC1–PC2 trajectories (Figure 8) further visualized the dynamic motion of each complex. The GroEL complex with F1243–0200 followed a tightly grouped trajectory path with minimal variation, whereas F2801–0179 displayed a widespread, irregular motion path, suggesting higher conformational variability. F3139-0927, F3385-2019, and F1864–0208 occupied intermediate spaces, showing moderate fluctuation with overall stable dynamic behavior.

To further elucidate the energetic landscape and thermodynamic stability of the GroEL-ligand complexes, FEL analysis was conducted based on the first two principal components (Figure 9). The FEL maps provide quantitative insights into the relative stability of different conformational states by visualizing the energy barriers and minima across the sampled conformational space. The energy contour plots revealed distinct patterns for each complex, with color gradients representing free energy variations from stable (blue) to unstable (red) regions. The F1243–0200 complex (Figure 9B) displayed well-defined, deep energy minima with narrow basins, indicating highly stable conformational states with limited transitions between energy wells. Similarly, F1864-0208 (Figure 9E) exhibited concentrated low-energy regions with sharp energy gradients, thereby reinforcing its conformational stability. In contrast, F2801-0179 (Figure 9D) showed a broader and more diffuse energy landscape with multiple shallow minima, which was consistent with the higher flexibility and conformational diversity observed in the PCA analysis. The energy surfaces for F3385-2019 (Figure 9A) and F3139-0927 (Figure 9C) exhibited intermediate characteristics, featuring moderately defined energy basins with occasional higher-energy transitions. FEL analysis corroborated the PCA findings, demonstrating that F1243–0200 and F1864–0208 occupy the most thermodynamically favorable conformational states with minimal energy barriers for local fluctuations. The deeper energy wells observed for these complexes suggest stronger binding interactions and a reduced likelihood of dissociation. Conversely, the flatter energy landscape of F2801–0179 indicates greater conformational freedom but potentially weaker binding affinity.

Collectively, the PCA and FEL results reinforce the dynamic stability trends observed in the previous RMSD, Rg, and hydrogen bond analyses. The F1243–0200 and F1864–0208 complexes demonstrated the most stable dynamic profiles and favorable energetic landscapes, while F2801–0179 exhibited greater flexibility and energy dispersion, which may correspond to looser binding or adaptability within the GroEL binding site. These findings enhance our understanding of the structural behavior of GroEL-ligand complexes and support the potential of F1243–0200 and F1864–0208 as robust and stable GroEL inhibitors.

To quantitatively assess the binding affinity and thermodynamic favorability of the GroEL-ligand complexes, MMPBSA calculations were performed on the last 50 ns of the MD trajectories. The binding free energies and their individual energy components are summarized in Table 4, providing detailed insights into the energetic contributions governing complex stability. MMPBSA analysis revealed significant variations in binding affinities among the five compounds. F2801–0179 demonstrated the most favorable binding energy (-317.677 ± 43.984 kJ/mol), followed closely by F1864-0208 (-269.698 ± 19.107 kJ/mol) and F3385-2019 (-267.150 ± 41.947 kJ/mol). F1243–0200 exhibited a moderately favorable binding energy (-309.769 ± 24.983 kJ/mol), while F3139–0927 showed the least favorable interaction (-26.711 ± 62.260 kJ/mol) (Figure 10A). Analysis of the individual energy components provided mechanistic insights into the binding interactions. Van der Waals interactions contributed significantly to complex stabilization across all systems, with F1864–0208 showing the strongest vdW component (-169.804 ± 10.133 kJ/mol), followed by F1243-0200 (-152.576 ± 15.917 kJ/mol) and F3385-2019 (-111.635 ± 16.468 kJ/mol). Electrostatic interactions predominantly displayed unfavorable contributions, likely owing to desolvation penalties upon complex formation. However, F2801–0179 showed a relatively smaller electrostatic penalty (-94.313 ± 43.984 kJ/mol) compared to other compounds. The polar solvation energies were generally unfavorable for most complexes, indicating the disruption of favorable water-protein and water-ligand interactions upon binding. F3139–0927 showed the least unfavorable polar solvation energy (-18.628 ± 67.633 kJ/mol), while F1243–0200 exhibited the most significant penalty (379.061 ± 53.568 kJ/mol). The SASA component contributed favorably to binding in all cases, with F2801–0179 showing the most favorable contribution (-13.273 ± 3.611 kJ/mol). To identify key binding hotspots and understand the molecular basis of ligand recognition, a per-residue energy decomposition analysis was performed (Figure 10B). Energy decomposition profiles revealed specific amino acid residues that contribute significantly to ligand binding. Several residues showed consistently favorable interactions across multiple compounds, indicating critical binding determinants within the GroEL active site.

The per-residue analysis identified both stabilizing and destabilizing interactions, with energy contributions ranging from approximately -30 to +30 kJ/mol per residue. In particular, residues E17, D40, D51, E66, D82, D86, E101, D120, E128, E138, D154, D163- D195, E208, D325, E387, E390, D397, E408, E423, E44, E459, E465, E474, E484, and D493 primarily contribute to the binding of these compounds to the GroEL receptor. These residues exhibited strong, favorable interactions (negative energy contributions), which appear to be crucial for maintaining stable ligand binding. The differential per-residue interaction patterns among the five compounds provide insights into the selectivity and specificity of ligand recognition. The combined MMPBSA and per-residue decomposition analyses demonstrate that F2801–0179 and F1864–0208 form the most thermodynamically stable complexes with GroEL, primarily driven by favorable van der Waals interactions and optimized desolvation effects. These findings complement the structural stability observations from the MDS and support the identification of these compounds as promising GroEL modulators. Substantial binding free energies of −269.698 kcal/mol and −267.150 kcal/mol, respectively, further validate their suitability as potential inhibitors of GroEL. In contrast, F3139–0927 showed a significantly weaker binding energy of −26.711 kcal/mol, suggesting lower complex stability during the final simulation window. To contextualize the pharmacological relevance of the identified lead candidates, comparison with established antibacterial agents currently used in clinical practice is informative. Benchmarking docking scores, interaction profiles, and MD-derived stability parameters against reference drugs such as ciprofloxacin, amoxicillin, or gentamicin would allow a clearer assessment of the relative binding strengths and dynamic stability of the proposed compounds, thereby enhancing the translational significance of the present findings.