A total of 7120 compounds were extracted from different databases during the research. OTVA Chemicals (https://www.otavachemicals.com/) provided 862 antiviral compounds, while the Enamine database provided 4201 compounds (https://enamine.net/compound-collections/real-compounds/real-database). Based on PubChem trials, 350 compounds were obtained, 21 compounds were obtained from PubChem records (https://pubchem.ncbi.nlm.nih.gov/docs/covid-19), and 1686 compounds were obtained from Life Chemicals (http://lifechemicals.com/). Using LigPrep (Schrödinger Release 2018), all ligand structures were optimized using the OPLS3e force field, maintaining the specified 3D structure and stereochemistry of the ionization states.

From the RCSB Protein Data Bank, a crystallographic construct of SARS-CoV-2 3CL was obtained with the inhibitor NK01-14, recognizable by the PDB ID 7TIA (Forouhar et al., 2022). Once acquired, this was incorporated into the Maestro suite’s Protein Preparation Wizard (Madhavi Sastry et al., 2013; Schrödinger Release, 2018). A hydrogen atom addition was performed, and adjustments were made to histidine residue protonation states. Additionally, hydrogen bonding patterns were analyzed to determine how water molecules should be removed. An OPLS3e force field was then used to conduct a restrained minimization procedure (Harder et al., 2016). Following the preparation of the protein, the co-crystallized ligand was positioned centrally in a grid box that was created using the grid generation wizard.

Once prepared, the assembled array of ligands was docked against the designed protein using Glide (Friesner et al., 2006). HTVS mode of Glide facilitates the virtual screening process. An initial procedure was conducted in which compounds were sorted according to the projected binding pocket and ranked based on Glide Scores (GScores). To enhance the precision of the screening, both Standard Precision (SP) and Extra Precision (XP) docking modes were implemented (Friesner et al., 2006).

The selected top 22 compounds from XP glide were evaluated and analyzed for their physicochemical and ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) characteristics. The SMILES notations of the significant hits were entered into the ProTox-II server (https://tox-new.charite.de/protox_II/) (Banerjee et al., 2018) to forecast toxicological properties. Furthermore, pharmacokinetic aspects such as absorption, distribution, metabolism, and excretion were evaluated through the SwissADME server (http://www.swissadme.ch/) (Daina et al., 2017). Additionally, BOILED Egg analysis was carried out using the SwissADME web-based tool.

Four leading molecules were selected for further analysis based on their Extra Precision (XP) docking scores, molecular interactions with the target protein, and ADMET parameters. When flexible ligand docking is implemented, protein structure can be anticipated using the induced fit docking (IFD) procedure. Compared to rigid ligand-receptor docking commonly encountered in structure-based virtual screening (Sherman et al., 2006a), this method provides a more detailed approach. In addition to the docking evaluation by Glide XP, free binding energy calculations, and ADMET analysis, IFD was performed on the chosen molecules. The process was executed in contrast to the standard protein SARS-CoV-2-3CL, adjusting the van der Waals energy to 0.50 and allowing Prime to generate and refine maximum poses within a radius of five degrees. Induced Fit Docking (IFD) relied on Glide for docking and Prime for refining binding poses following the standard parameters for these modules. A total of 20 poses were generated for each compound (Sherman et al., 2006b). Although this approach does not perfectly mimic biological conditions, it provides valuable insight into the molecule’s stability within a specific binding pocket when presented with various poses and frames.

The Molecular Mechanics Generalized Born Surface Area (MMGBSA) technique was used to compute the relative binding free energy (ΔG bind) for each ligand molecule (Jacobson et al., 2004). This evaluation aimed to determine the binding affinity of the ligand with the receptor, achieved through the Prime module. Using the VSGB 2.0 solvation model and the OPLS3e force field without modification (Brunt et al., 2022), the free binding energy of the top four molecules (hits) has been calculated. Binding affinity is a free energy that includes entropy and enthalpy components. Due to the generalized Born approximation, conformational entropy upon ligand binding was excluded from the MM/GBSA calculations. This provides a quicker but more approximate resolution to the Poisson–Boltzmann equation (Genheden and Ryde, 2015; Wang et al., 2018).

Following the MMGBSA analysis, molecular dynamics simulations were performed on the most promising molecule in order to gain further insight into their stability under biological conditions. Molecular Dynamics (MD) studies can provide insight into the stability of complexes within a solvent system (Bowers et al., 2006). The system was set up in an orthorhombic box with dimensions of 10 Å × 10 Å x 10 Å, derived using the buffer size approach, which also aids in minimizing the box’s volume. The Desmond system builder was used to establish a TIP3P solvent model. In these MD simulation studies, proteins and ions were modeled using the latest OPLS3e force field developed by Schrodinger Inc. (Mark and Nilsson, 2001; Shivakumar et al., 2010; Padhi et al., 2022). Sodium chloride, with Na + as the cation and Cl-as the anion, was incorporated at a concentration of 0.15 M.

A Desmond Molecular Dynamics module was used to operate the system for a runtime of 100 nanoseconds, producing approximately 1000 frames. A molecular dynamic simulation occurred at a stable temperature of 300 K and a pressure of 1 bar, using the NPT ensemble class. The simulation process began after the system model reached a state of relaxation. The molecules were analyzed using MMGBSA (Molecular Mechanics Generalized Born Surface Area) after the molecular dynamic studies had been concluded, and binding free energies were calculated every 10th frame throughout the entire 1000 frames of the simulation (Chaudhary and Dikshit, 2023).

As part of the High-Throughput Virtual Screening (HTVS) analysis, 5962 compounds from combined libraries were selected. The docking parameters of all molecules can be found in the Supplementary Table S1. As well as docking the selected compounds with the receptor, the native ligand of the receptor was docked separately, allowing for a comparison of the binding affinity between the compounds and the receptor. The initially anticipated binding site was utilized for high-throughput virtual screening (HTVS). In the end, 202 compounds were docked with Standard Precision (SP), and 22 compounds with Extra Precision (XP) docking were docked with SP (Figure 1). The dock cutoff score for SP docking by XP was set at −5.0 kcal/mol based on XP docking outcomes. The glide scores for HTVS and SP ranged from −7.492 to 1.148 Kcal/mol and −5.854 to −2.680 Kcal/mol, respectively, while those for XP ranged from −6.081 to −1.933 Kcal/mol. Most compounds from the Life Chemical library demonstrated better binding affinity than others. The compounds EN1036, F6548-4084, F6548-1613, and PUBT44123754 demonstrated high binding affinity, exceeding the standard ligand. Based on Glide XP scores, 16 molecules showed a better binding affinity than their native ligand (Table 1). The XP glide score for the native ligand was −3.275 kJ/mol, lower than the top four selected compounds. The compounds F6548-4084, F6548-1613, EN1036, and PUBT44123754 recorded docking scores of −5.772, −5.71, −5.78, and −5.003 Kcal/mol, respectively. Among these, compound F6548-1613 formed a hydrogen bond with Gln189 via the oxygen atom of a carboxamide attached to a furan ring. Distinct hydrophobic interactions are observed with amino acids such as Leu27, Met49, Leu50, Cys145, Met165, Leu167, Pro168, Val186 and Ala191. Polar interactions are also noted with Thr25, Thr26, Hie41, Ser46, Asn142, Gln189, Thr190, and Gln192 (refer to Figure 2). Furthermore, charge interactions, both positive and negative, are identifiable (Table 2). Alongside these, non-bonded interactions with amino acids Hie45, Asn142, Gly143, Cys145, Met165, Glu166, Val186, Arg188, and Gln189 are evident in both the native ligand and compound F6548-1613 (refer to Figure 2). The same applies to all the selected compounds, F6548-4084, F6548-1613, EN1036, and PUBT44123754, which showed common interactions with Leu27, Hie41, Met49, Gly143, Asn142, Cys145, Met165 Glu166, Arg188, and Gln189 residues (See Supplementary Figure S1). Induced fit docking (IFD) and MMGBSA calculations were performed alongside the native drug to examine how the selected compounds interact and conform within the protein target’s binding pose.

It was determined that the protein exhibits substantial conformational variation, which was analyzed through the Induced Fit Docking (IFD) strategy, utilizing GLIDE-XP results. To understand the potential binding conformational attributes of the top four inhibitors, EN1036, F6548-4084, F6548-1613, and PUBT44123754, IFD computations were performed (see Table 4; Figure 4 and See Supplementary Figure S4). Figure 4 illustrates how the inhibitor docks within the catalytic pockets of the 7TIA enzyme. Compound F6548-1613 displayed the most potent binding affinity from the shortlisted compounds, predicting 16 different poses with 3CLpro. The primary pose had the highest IFD score of −675.26 Kcal/mol, as indicated in Table 4. According to the investigation, pose 13, with an Induced Fit Docking (IFD) score of −668.47 Kcal/mol, exhibits the highest number of advantageous bondings and non-bondings with the primary protease. F6548-1613 established four hydrogen bonds with the protein. These included interactions with the 4th oxygen group of the benzene ring by Gly143, the NH group of the carboxamide substituted benzene ring by Gln189 with two groups, the NH group of the carboxamide substituted benzene ring, and the oxygen atom of the carboxamide linked with the furan ring, and the 3rd oxygen atom of the furan ring with the Thr190 residue. These interactions are illustrated in both 2D and 3D diagrams (refer to Figure 4). The IFD scores and 3D interactions can be found in Table 4 and Figure 4. Compound F6548-1613 was observed to interact with the same amino acid residue loops as the native ligand (co-crystal ligand), comprising residues like Met165, Glu166, Leu167, Pro168, Val186, Asp187, Arg188, Gln189, Thr190, Ala191, and Gln192. Furthermore, specific amino acid residues like Hie41, Met49, and Tyr210 showed common interaction in Compound F6548-1613 and standard drugs.

Compound F6548-1613, one of the top hits, displayed superior binding free energies compared to standard compounds, exhibiting the highest binding energy of −65.72 kcal/mol when interacting with 7TIA. A native ligand, on the other hand, showed ΔG binding energies of −58.59, −63.29, and −38.26 kcal/mol. The binding energies of GbindHbond, GbindvdW, GbindLipo, and GbindCoulomb were significant (see Supplementary Table S4). Due to its optimal balance of several types of interactions, F6548-1613 emerges as the most effective candidate for targeting SARS-CoV-2 3CLPRO protease. This compound’s MMGBSA dG Bind vdW score is the lowest (−54.55), indicating robust van der Waals interactions, contributing significantly to its high binding affinity. Despite Compound EN1036 having the lowest GbindCoulomb score, which suggests potent Coulombic or electrostatic interactions, Compound F6548-1613 still demonstrated the highest overall binding affinity based on a summation of its various interactions.

Furthermore, Compound F6548-1613 displayed the lowest GbindLipo score, indicating the presence of lipophilic or hydrophobic interactions that enhance protein-ligand binding and complex stability. Again, its GbindPacking score was among the lowest, indicating superior packing interactions and a well-fitted ligand within the protein’s binding pocket. However, despite F6548-1613 and PUBT44123754 showing relatively less negative GbindHbond scores, suggesting fewer or weaker hydrogen bonds, F6548-1613’s binding affinity remained superior. This dominance can be attributed to the combination of all its interaction types. GbindSolvGB (Solvation Energy) score, a measure of the free energy of solvation, indicates how favorably a ligand is solvated before entering the binding pocket. In this regard, F6548-1613 demonstrated a highly positive score (28.92), denoting a favorable solvation energy. All compounds, including F6548-1613, displayed positive GbindCovalent scores, a trait often associated with the absence of covalent bonds or the presence of destabilizing covalent interactions.

Nevertheless, F6548-1613 had the lowest positive value, suggesting the least destabilizing influence from potential covalent interactions. Molecular Dynamics (MD) was performed to analyze the stability of F6548-1613 within the protein target’s binding pose, providing insight into the molecule’s stability in the binding pocket against 3CL protease. This additional layer of analysis further substantiates F6548-1613’s potential as an effective inhibitor of the SARS-CoV-2 3CLPRO protease.

The presented graph (Figure 5A) depicts Root Mean Square Deviation (RMSD), a commonly employed statistic for evaluating disparities in atomic positions across various time points. The x-axis corresponds to the time in picoseconds (ps), while the y-axis displays RMSD values in nanometers (nm). The graph exhibits three distinctive simulation runs—Run-1 in black, Run-2 in red, and Run-3 in green. Observations from the graph indicate that RMSD readings peak during the initial 20,000 picoseconds (or 20 nanoseconds) for all three runs, suggesting potential initial fitting adjustments. Following this phase, each run’s RMSD tends to stabilize, fluctuating around specific values. This stabilization implies that the ligand’s positioning relative to the protein reaches a relatively constant state, with occasional deviations possibly arising from interactions within the simulation environment or conformational changes. Notably, Run-3 (green) displays the lowest average RMSD, indicating that, in contrast to the preceding runs, the ligand in this simulation maintains a more consistent position or conformation throughout. In comparison, runs 1 (black) and 2 (red) exhibit similar patterns with slightly higher RMSD values, implying increased conformational flexibility. This flexibility may be attributed to the ligand’s free end, lacking hydrogen bond formation with amino acid residues. Despite the initial instability observed in all three simulation runs, they eventually establish a dynamic equilibrium, where the ligand’s position fluctuates around a relatively constant value, as illustrated in the RMSD plot.

Graph (Figure 5B) visually represents the Root Mean Square Deviation (RMSD) of C-alpha atoms over time, obtained through least squares fitting to a protein in a molecular dynamics simulation. The y-axis reflects RMSD in nanometers (nm), while the x-axis corresponds to time in picoseconds (ps). The graph encompasses four distinct traces: the Apo form in black, the control run without any ligand; Run-1 in red; Run-2 in green; and Run-3 in blue, representing three simulation runs for the protein-ligand complex. Examining the graph, it is evident that the RMSD values for each run commence at relatively low levels, indicating minimal initial deviation of the protein’s C-alpha atoms from the reference structure, which is the starting point of the simulation. As time progresses, the RMSD values undergo shifts yet consistently maintain within a narrow range (approximately 0.1–0.4 nm). This suggests a lack of significant conformational changes in the protein structures throughout these simulation runs. Notably, the black trace representing the Apo form exhibits the least fluctuation, implying that the protein remains relatively stable in its original state or when devoid of any ligand.

In contrast, Runs 1, 2, and 3 display more significant fluctuations attributed to the presence of the ligand. Despite this increased movement, the variations observed in these runs remain confined within a limited range, suggesting the absence of substantial structural alterations. In summary, the RMSD graph provides insights into the stability of the protein’s secondary structure during the simulation. The Apo run demonstrates minimal deviation, while the other runs, influenced by interactions with the ligand, exhibit increased movement within a constrained range. This overall pattern indicates that the protein maintains structural stability throughout the simulation, even in the presence of a ligand.

The provided graph (Figure 5C) illustrates a protein’s Root Mean Square Fluctuation (RMSF), portraying the flexibility of specific residues during a simulation. Higher RMSF values indicate increased flexibility or mobility of these residues within the protein structure. The y-axis represents RMSF values in nanometers (nm), showing the extent of fluctuation for each residue. At the same time, the x-axis lists the residues of the protein sequence from the N-terminus to the C-terminus. The data from four distinct runs are depicted in the plot: Run-1 (red), Run-2 (green), Run-3 (blue), and Apo (black), with the Apo run representing the unbound or native state of the protein without a ligand. Notably, the Apo run serves as a baseline, displaying consistent fluctuation across residues. Across all runs, a similar fluctuation pattern is observed, suggesting that the overall flexibility of the protein remains relatively constant throughout the simulation runs. Peaks in the graph represent residues that exhibit higher flexibility, and these peaks are consistent across all runs, indicating inherent flexibility in specific regions of the protein. Most residues demonstrate moderate fluctuations ranging between 0.1 and 0.3 nm, with fluctuations typically staying below 0.4 nm. Towards the end of the residue sequence, prominent peaks emerge where RMSF values are notably higher. This is attributed to relatively lesser interaction in the terminal residues. Notably, residues such as Tyr154 and Arg222 consistently show maximum fluctuation across all three runs. Additionally, a peak at Leu50 in the Apo form, which is not as pronounced in the other runs, suggests that this residue is stabilized by adding the ligand in the binding pocket. While small differences exist between runs, the overall patterns are similar, indicating that the ligand or the conditions tested do not significantly alter the protein’s flexibility at the residue level. This suggests that the protein maintains a consistent level of flexibility across various simulation conditions.

Graph (Figure 5D) illustrates the radius of gyration for a protein, indicating how its mass is distributed from its center of mass and reflecting the molecule’s compactness. The x-axis represents time in picoseconds (ps), while the y-axis displays the radius of gyration in nanometers (nm). The four simulation runs represented by lines are Apo (black), Run-1 (red), Run-2 (green), and Run-3 (blue). Throughout all simulation runs, the radius of gyration exhibits variations within a relatively confined range, approximately 2.15 nm–2.35 nm. These variations suggest changes in the protein structure’s compactness over time, albeit within a relatively short range. Notably, the average radius of gyration for the Apo form appears to be the lowest, indicating that the protein is more compact in its unbound or reference state. Runs 1, 2, and 3 show comparable trends as the gyration radius gradually increases. This suggests that these proteins experience a loss of compactness, potentially due to unfolding or conformational changes induced by the interacting ligand. The overall plot indicates that the protein structures in simulation undergo some conformational change, although this change is not notably high. In summary, the radius of the gyration graph provides insights into the dynamic changes in the protein’s compactness over time in different simulation conditions. While there is observable variation, the alterations remain within a limited range, and the protein structures experience some degree of conformational change, particularly when interacting with the ligand.

Graph (Figure 5E) illustrates the Solvent Accessible Surface Area (SASA) of the protein, a metric measured in square angstroms (Ų) that reflects the surface area accessible to a solvent, typically water in biological systems. The y-axis represents SASA values, while the x-axis depicts time in picoseconds (ps). The graph encompasses four traces corresponding to different simulation runs: Apo (black), Run-1 (red), Run-2 (green), and Run-3 (blue). The dynamic nature of the protein’s conformational changes is evident as SASA varies over time for each run, reflecting alterations in the solvent-exposed area of the molecule. The Apo form (black) notably exhibits the highest average SASA, indicating a greater average exposed surface area when the protein is unbound.

In contrast, the Runs 1, 2, and 3 SASA values consistently appear somewhat lower than the Apo form, suggesting a more compact or less exposed structure in these scenarios. The similar fluctuation patterns across these runs imply a comparable behavior in terms of solvent accessibility. Despite dynamic fluctuations, the SASA plot reveals no sharp rises or falls for any of the runs, indicating that significant unfolding or refolding does not occur during the simulation period.

Furthermore, the SASA for each run consistently remains lower than that of the Apo form, implying an increase in the folding of the protein upon interaction with the ligand. In summary, the SASA graph provides a comprehensive view of the protein’s dynamic behavior, suggesting that while there are variations in solvent accessibility over time, the protein’s overall structure remains stable. The protein undergoes conformational fluctuations without experiencing substantial unfolding or refolding during the simulated conditions.

Graph (Figure 5F) illustrates the dynamic evolution of hydrogen bond formation over time in picoseconds (ps) across three distinct runs (Run-1, Run-2, and Run-3) involving compound F6548-1613 complexed with the target protein 3CLpro. The y-axis quantifies the number of hydrogen bonds formed between the protein and ligand, while the x-axis represents the temporal progression in picoseconds. Each run is differentiated by black, red, and green bars, corresponding to Run-1, Run-2, and Run-3. Examining the graph, it is evident that the total count of hydrogen bonds for each run exhibits notable variability over time, showcasing irregular fluctuations in bond numbers. The data portray a seemingly erratic pattern with no discernible trend or pattern across the simulation duration for any run. This suggests a certain degree of randomness in both the creation and disruption of hydrogen bonds, indicating continuous and variable occurrences of bond formation and breaking.

Moreover, the absence of a clear distinction between the runs implies that the conditions or factors influencing hydrogen bond formation were likely consistent across all three experimental runs. Exploring potential reasons for the absence of significant differences between each run could provide valuable insights into the intrinsic variability of the observed events. This detailed analysis could contribute to a deeper understanding of hydrogen bond dynamics’ continuous and dynamic nature under the specified experimental conditions.