Ligand-based virtual screening to discover potential inhibitors of SARS-CoV-2 main protease

Abstract

The main protease (M^pro, also known as 3CL^pro), a pivotal enzyme of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has been considered a prime target for drug development due to its crucial role in viral replication and transcription. Importantly, a high degree of conservation in more than 13 million SARS-CoV-2 sequences affords M^pro as a promising target for antiviral therapy to impede the genetic evolution of SARS-CoV-2. In this work, ∼16 million compounds from various small molecule databases were screened using ligand-based virtual screening (LBVS) with boceprevir as the reference compound to identify new small molecule inhibitors of M^pro. Boceprevir [hepatitis C virus (HCV) drug] has been repurposed as a drug candidate against M^pro activity (IC₅₀ = 4.13 ± 0.61 μM). The lead compounds exhibiting higher binding affinities (−9.9 to −8.0 kcal mol⁻¹) than boceprevir (−7.5 kcal mol⁻¹) were identified from a library of 850 compounds using molecular docking. Furthermore, molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) analysis depicted ChEMBL144205 (C3), ZINC000091755358 (C5), and ZINC000092066113 (C9) with binding affinities of −65.2 ± 6.5, −66.1 ± 7.1, and −67.3 ± 5.8 kcal mol⁻¹, respectively, as high-affinity binders to M^pro. The identified compounds displayed a favourable drug-likeness profile without violating Lipinski's rule of five. Molecular dynamics (MD) simulations revealed the higher structural stability and reduced residue-level fluctuations in M^pro upon binding of C3, C5, and C9 as compared to apo–M^pro and M^pro–boceprevir. Notably, conformational clustering and FEL analyses depicted hydrogen bond interactions of C3 with Thr26, oxyanion hole residues (Asn142 and Gly143), the catalytic residue (Cys145), and Glu166 of M^pro, suggesting its strong binding affinity and potential inhibitory effect. The integrated computational methodology employed in this work identified promising lead compounds against M^pro activity, which warrants further experimental validation to develop them as antiviral agents against SARS-CoV-2.

---

1. Introduction

The fatal coronavirus disease (or long COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has adversely infected ∼778 million people worldwide and resulted in nearly 7.1 million mortalities as of July 13, 2025.¹ COVID-19 patients exhibit symptoms of a cold, cough, chest uneasiness, fever, etc.² Several repurposed drugs such as baricitinib,³ favipiravir,⁴ hydroxychloroquine,⁵ molnupiravir,⁶ paxlovid,⁷ and remdesivir⁸ have been used to treat the symptoms of SARS-CoV-2. However, these antiviral medications have low efficiency and cause unwanted side effects on patients.⁹ Thus, researchers are focusing on new therapeutic candidates against various druggable SARS-CoV-2 proteins.

SARS-CoV-2, belonging to the genus Betacoronavirus and family Coronaviridae, is an enveloped positive-sense and single-stranded RNA virus.¹⁰ Likewise, bat coronaviruses, SARS-CoV, and MERS-CoV, SARS-CoV-2 have a similar genetic code comprising accessory, structural, and non-structural proteins.¹¹ Notably, SARS-CoV-2 main protease [M^pro, also known as 3CL^pro (3-chymotrypsin-like protease)] is the most prominent drug target among several other promising targets of SARS-CoV-2 due to its highly conserved structure with no human homologs and crucial role in virulence.^12–15 Ambrosio et al. noted a high degree of conservation in more than 13 million sequences of SARS-CoV-2 and reported that M^pro is a promising target for antiviral therapy to impede the genetic evolution in SARS-CoV-2.¹⁶ M^pro is also known as cysteine protease and plays a key role in processing the proteolytic cleavage of viral polyprotein 1ab at 11 different sites.¹⁷

M^pro is a homodimer and consists of two protomers arranged nearly perpendicular to each other.¹⁸ Each protomer is comprised of three domains of 306 residues and contains a catalytic dyad (His41 and Cys145) situated in a cleft between domain I (10–99) and domain II (100–182). In the proteolytic cleavage of polyprotein 1ab, His41 acts as a base and accepts a proton from the sulphur of Cys145, generating a thiolate intermediate that behaves as a nucleophile when it reacts with polyprotein.¹⁹ The catalytic domain is connected to domain III (198–303) with a longer loop (183–197). M^pro consists of conserved binding subpockets, S1 (Phe140, Leu141, Ser144, Cys145, His163, His164, Met165, Glu166, His172), S2 known as the catalytic center (Thr25, Thr26, Leu27, His41, Gly143, Cys145), S3 (His41, Met49, Tyr54, His164, Met165, Asp187, Arg188, Gln189), S4 (Met165, Glu166, Leu167, Pro168, Asp187, Arg188, Thr190, Gln192), and S5 (Thr24, Thr25, Thr26, His41, Cys44, Thr45, Ser46, Met49).^20,21 M^pro functions as an asymmetric dimer following a flip-flop mechanism with only one monomer active at a time.²² It likely undergoes a half-site catalytic cycle, alternating between acylated and deacylated states in its two subunits.^22,23 The dimerization of M^pro is of paramount importance as it significantly influences both enzymatic activity and viral replication.^18,24–26 Consequently, inhibition of M^pro can potentially stop viral replication.

Several small-molecule inhibitors against M^pro activity have been developed.^27–31 Ma et al. reported boceprevir [hepatitis C virus (HCV) drug] as a potent inhibitor of M^pro by employing the drug repurposing approach.³² The in vitro studies displayed significant inhibition of M^pro activity (IC₅₀ = 4.13 ± 0.61 μM). The inhibitory mechanism and non-covalent interactions of boceprevir with M^pro scrutinized using molecular docking and molecular dynamics (MD) simulations revealed that boceprevir interacted with the catalytic dyad and binding site residues of M^pro through hydrogen bonds and hydrophobic contacts.²¹ MD simulations highlighted the structural stability of M^pro and depicted the inhibition of M^pro activity in the presence of boceprevir.²¹ Numerous other studies employed MD simulations³³ to illuminate the inhibitory mechanism of boceprevir against M^pro and evaluated small-molecule inhibitors based on boceprevir using in vitro studies.³⁴

Recently, Singh et al. identified six remarkably effective M^pro inhibitors (bexarotene, diacerein, KT185, ledipasvir, simeprevir, and WIN-62577) from the virtual screening of a library of 8000 compounds.³⁵ The screened molecules efficiently inhibited the M^pro activity with IC₅₀ values ranging from 0.64 to 11.98 μM and EC₅₀ values from 1.51 to 18.92 μM. Notably, the binding of minocycline (one of the compounds screened against the dimeric interface of M^pro) to the allosteric site in close vicinity of the dimeric interface of M^pro interferes with both allosteric and active sites of M^pro, disrupting its dimeric stability and impairing its catalytic activity. Khamto et al. elucidated the inhibitory potential of 22 flavonoids against M^pro, including cytotoxicity toward Vero cells.³⁶ The in vitro enzyme assay highlights that tectochrysin (IC₅₀ = 24 μM), compound 9 (IC₅₀ = 19.87 μM), panduratin A (IC₅₀ = 13.28 μM), and genistein (IC₅₀ = 17.98 μM) showed strong anti-proteolytic activity towards M^pro in comparison to the control baicalein (IC₅₀ = 86.57 μM). Among these flavonoids, genistein extracted from Millettia brandisiana had low cytotoxicity on Vero cells with cytotoxic concentration 50% (CC₅₀) > 50 μM. The molecular mechanics Poisson–Boltzmann (generalized-Born) surface area [MMPB(GB)SA] analysis depicted the favourable binding of flavonoids to M^pro with binding free energy from −35.54 to −22.41 kcal mol⁻¹. Jin et al. performed a structure-based virtual screening of 18 263 traditional Chinese medicines (TCM) against M^pro.³⁷ The in vitro experiments confirmed the three TCM compounds (CAS No. 18085-97-7, 521-61-9, and 490-31-3) as M^pro inhibitors with IC₅₀ and EC₅₀ values of (4.64 ± 0.11 and 12.25 ± 1.68 μM), (7.56 ± 0.78 and 15.58 ± 0.77 μM), and (11.16 ± 0.26 and 29.32 ± 1.25 μM), respectively. The MD simulations depicted the stability of the three complexes during the last 100 ns and non-covalent interactions with key residues His41, Asn142, Gly143, Met165, Glu166, and Gln189 of M^pro.

Ambrosio et al. employed virtual screening and in vitro studies to discover potent inhibitors of M^pro activity as repurposed drugs from 8000 FDA-approved and investigational drugs.¹⁶ The in silico and in vitro studies highlighted eight promising inhibitors against M^pro, and betrixaban (an oral anticoagulant, IC₅₀ = 0.9 ± 0.0053 μM), potentially blocked the proteolytic activity of M^pro without any cytotoxic effect. Samanta et al. virtually screened remdesivir analogues as M^pro inhibitors from the PubChem database and deciphered their inhibitory mechanism against M^pro activity using molecular docking and MD simulations.³⁸ The binding free energy analysis using MMPB(GB)SA methods depicted three compounds with PubChem IDs [134133102, 76314404, and 58059494] as potent inhibitors of M^pro activity. Li et al. illuminated the binding interactions of lopinavir, saquinavir, ritonavir, and PF-07321332 (nirmatrelvir) with M^pro using conventional and replica exchange MD (REMD) simulations.³⁹ Among all inhibitors, PF-07321332 displayed the highest binding affinity (−26 kJ mol⁻¹) with M^pro. Importantly, PF-07321332 specifically binds to the catalytic site residues through multiple hydrogen bonds with His163 and Glu166 in comparison to other inhibitors, as they bind to multiple sites in the M^pro structure. Notably, PF-07321332, co-formulated with ritonavir (known as Paxlovid),⁴⁰ is in clinical use as an M^pro inhibitor and displayed ∼89% efficacy against severe COVID-19.⁴¹

Although several potent inhibitors of M^pro activity have been discovered,^27,42 however, their failure in clinical trials emphasizes the urgent need to find new potent therapeutic candidates against M^pro activity. In comparison to the significant time, costs, and efforts required for de novo development, the drug repurposing of existing inhibitors is a strategic choice for the rapid development of antiviral agents, leveraging their available pharmacological applications, efficacy, and safety profile. Drug repurposing (also referred to as therapeutic switching) has promise for accelerating the discovery of effective treatments by re-evaluating clinically approved molecules for their activity against viral targets such as M^pro. The global response to the COVID-19 pandemic has been marked by rapid and multifaceted efforts, including the repurposing of diverse classes of therapeutic agents and the accelerated development of various types of emergency-use vaccines.^43,44 The drug repurposing technique has been widely used to identify potent inhibitors against Alzheimer's disease, Parkinson's disease, and type 2 diabetes.⁴⁵ Thus, a ligand-based virtual screening (LBVS) method using boceprevir as a reference compound was employed in this work to identify the lead compounds from various small-molecule databases against M^pro activity, which can be repurposed as antiviral agents. Notably, multiple databases have been utilized in this work to expand the range of identification of potential non-covalent inhibitors of M^pro activity, in contrast to previous studies,^16,35,37 which restricted their virtual screening efforts to a single or double database. Importantly, a top-hit compound C3 identified as a potent inhibitor of M^pro in this work, previously displayed inhibition against cell proliferation, triggers cell cycle arrest, and induces apoptosis of human breast cancer cells.⁴⁶ Importantly, previous studies highlighted that C3 at <20 nM effectively suppresses the malignant phenotype in established neuroblastoma (NB) cell lines as well as primary NB cells,⁴⁷ which highlights its potential as a repurposed drug candidate against M^pro activity.

2. Computational details

2.1. Library screening using SwissSimilarity

A library of small molecules was virtually screened by SwissSimilarity.⁴⁸ SwissSimilarty is a web server used to screen small molecule databases like ChemBridge,⁴⁹ ChEMBL,⁵⁰ ChEBI,⁵¹ DrugBank,⁵² Ligand Expo,⁵³ GLASS-GPCR,⁵⁴ HMDB,⁵⁵ ZINC,⁵⁶etc. To get both 2D and 3D shape similarity, a single approach named combined (a combination of electroshape and FP2) in the SwissSimilarity tool has been applied to screen the small molecules from various databases.

2.2. Drug-likeness assessment

The evaluation of drug-likeness was performed using the SwissADME⁵⁷ web server employing Lipinski's rule of five by taking into account several parameters. These parameters include molecular weight (MW) below 500 Daltons (Da), a partition coefficient (ilogP) below five, less than ten hydrogen bond acceptors (HBA ≤ 10), less than five hydrogen bond donors (HBD ≤ 5), and a topological surface area (TPSA) below 140 Ų.

2.3. Structure preparation of M^pro and screened compounds

The X-ray crystal structure of M^pro homodimer (PDB ID: 6Y84) was chosen as the starting structure in this work and was retrieved from the Protein Data Bank (PDB).⁵⁸ The M^pro structure was reported at room temperature with a resolution of 1.39 Å. The cleaned M^pro structure without any heteroatoms and water molecules was used for docking. The structures of screened compounds were saved in SMILES format from the SwissSimilarity tool (Table S1), which was changed to SDF format using Open Babel.⁵⁹ Furthermore, the structures were converted into PDBQT format by performing energy minimization with the MMFF94 force field⁶⁰ using Open Babel.⁵⁹ The generated PDBQT files of screened compounds were submitted to molecular docking using AutoDock Vina.⁶¹

2.4. Molecular docking

To elucidate the binding affinity of 850 compounds with M^pro, molecular docking was performed employing the AutoDock Vina version 1.1 package.⁶¹ For this purpose, the size of the docking grid was set at 48 Å × 44 Å × 46 Å with the grid center defined at 11.615, −5.694, and 22.049 in x, y, and z dimensions, respectively. The Lamarckian genetic algorithm (LGA)⁶² was employed in molecular docking that performs a global search with the genetic algorithm and a local search using the Solis & Wets algorithm.^63,64 The exhaustiveness for the global search was kept at 100 with a default root-mean-square deviation (RMSD) of 1.0 Å. The cut-off for binding affinity was kept ≤−7.5 kcal mol⁻¹ to screen the top hits as the reference compound, boceprevir, displayed a binding affinity of −7.5 kcal mol⁻¹ with M^pro. The top-ranked binding poses from Vina [i.e., the pose with the lowest (most negative) binding affinity] were chosen for subsequent analysis, including MD simulations. The docking analysis was performed using AutoDock tools (ADT),⁶⁵ PyMOL,⁶⁶ and LigPlot+.⁶⁷ Based on highly negative binding affinity and a large number of binding interactions with the key residues of the binding site of M^pro, the top nine compounds were retained in comparison to boceprevir (−7.5 kcal mol⁻¹). The docked poses of M^pro complexes with a binding affinity ≤−7.5 kcal mol⁻¹ and interactions with key residues of the binding site of M^pro were chosen for further studies (Table S2). To validate the AutoDock Vina results, redocking of the screened compounds against M^pro was performed using Glide.⁶⁸

2.5. Binding free energy evaluation to screen top hit compounds with M^pro using molecular mechanics Poisson–Boltzmann surface area (MM-PBSA)

From molecular docking screening, the top nine compounds were chosen for binding free energy evaluation using the MM-PBSA method. The MM-PBSA was employed to screen the top hits out of nine compounds obtained from molecular docking. The top nine docked M^pro complexes and M^pro–boceprevir were subjected to energy minimization employing the steepest descent algorithm. NVT and NPT equilibration were performed sequentially during the equilibration phase, each for 1 ns. Then, all ten complex systems were run for 10 ns MD simulations. After MD simulations, the resulting trajectories of complex systems underwent binding free energy (ΔG_binding) calculations between M^pro and screened compounds using the g_mmpbsa tool.⁶⁹ Among the various sets of atomic radii in the g_mmpbsa tool, the bondi set was employed for the MM-PBSA calculations. The bondi set of atomic radii has been widely used in previous studies.⁷⁰ The contribution of conformational entropy was ignored during the binding free energy calculations, following the literature.⁷¹ The top three M^pro complexes were chosen, with highly negative binding free energies as compared to the M^pro–boceprevir complex (−51.8 ± 6.5 kcal mol⁻¹). Furthermore, the MD simulations of apo–M^pro and M^pro–boceprevir were performed along with the top three shortlisted complexes to assess their inhibitory potential against M^pro activity.

2.6. MD simulations and analysis protocol

MD simulations were performed using the GROMACS package⁷² to scrutinize the conformational stability and binding interactions of M^pro with the top three compounds screened from MM-PBSA. A total of five MD simulation systems were prepared, namely: apo–M^pro, M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 (Table 1 and Fig. 1). All the MD simulations were performed using the OPLS–AA/L force field,⁷³ and parameters of this force field for the reference and screened compounds were generated using the LigParGen Server.⁷⁴ All systems were solvated by TIP3P water⁷⁵ within a cubic box of dimensions 11.69 nm × 11.69 nm × 11.69 nm, with a minimum 1.0 nm distance of the receptor from the box edge. The number of water molecules added to each system is listed in Table 1. The protonation states of M^pro residues were predicted using PROPKA⁷⁶ at a physiological pH of 7.4 and assigned during the system preparation in GROMACS (Table S3), where His41 was set as neutral with only protonated delta position (HID type) according to Nutho et al.⁷⁷ The overall neutrality in each system was maintained with the addition of a stipulated number of Na⁺ and Cl⁻ counterions along with the 0.15 M NaCl (Table 1). The energy minimization of each system was completed by the steepest descent algorithm and then equilibrated under NVT and NPT conditions (temperature at 310 K and pressure at 1 bar).⁷⁸ Both NVT and NPT equilibration steps were performed sequentially during the equilibration phase, each for 1 ns, resulting in a total of 2 ns of equilibration before the production MD run. First, the system was equilibrated under NVT to stabilize the temperature, followed by NPT equilibration to stabilize the pressure and density of the system. The equilibrated systems were subjected to MD simulations (100 ns each), preserving a target pressure of 1 bar using a Parrinello–Rahman barostat⁷⁹ and a temperature of 310 K using a Nosé–Hoover thermostat.⁸⁰ The LINCS algorithm⁸¹ was used to constrain all the bonds involving hydrogen atoms, and the SETTLE algorithm⁸² to constrain the bond lengths in the solvent (water) molecules with a 2 fs integration step. The long-range electrostatic interactions were evaluated using the particle mesh Ewald (PME) method, and a 1.2 nm cut-off was used to estimate the short-range van der Waals interactions.⁸³ To ensure the reproducibility of MD simulations, repeat simulations of apo–M^pro, M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 have been performed using different initial velocities. However, data from one representative trajectory were analyzed to illuminate the binding interactions of boceprevir and top hit compounds with M^pro.

Table 1 MD simulation details of apo–M^pro and M^pro complexes

System	Simulation time^b (ns)	Box dimensions (nm)	Total number of water molecules in the simulation box	Total number of Na^+ ions added in the simulation box
a SARS-CoV-2 M^pro (PDB ID: 6Y84). b The simulations were performed with an all-atom OPLS force field and TIP3P water model with 0.15 M NaCl concentration.
apo–M^pro^a	100 × 3	11.69 × 11.69 × 11.69	50 269	8
M^pro–C3	100 × 2	11.69 × 11.69 × 11.69	50 261	8
M^pro–C5	100 × 2	11.69 × 11.69 × 11.69	50 265	8
M^pro–C9	100 × 2	11.69 × 11.69 × 11.69	50 267	8
M^pro–boceprevir	100 × 2	11.69 × 11.69 × 11.69	50 255	8

---

Fig. 1 Chemical structures of boceprevir and top hit compounds as M^pro inhibitors identified using LBVS, molecular docking, and MM-PBSA from various small molecule databases.

The MD trajectories were analyzed using various GROMACS utilities.⁷² The structural variations in apo–M^pro and M^pro complexes were analyzed using the RMSD, radius-of-gyration (R_g), and root-mean-square fluctuation (RMSF) by employing GROMACS tools, i.e., gmx rms, gmx gyrate, and gmx rmsf, respectively. The GROMACS tools “gmx hbond” and “gmx sasa” were employed to evaluate the number of hydrogen bonds and solvent accessible surface area (SASA), respectively. The conformational clustering was performed using the “gmx cluster” utility of GROMACS employing the Daura et al.⁸⁴ algorithm at 0.11 nm RMSD cut-off. The center of mass (COM) distances between M^pro residues and compounds were evaluated using “gmx distance” utility.

The principal component analysis (PCA) was performed for M^pro in the absence and presence of small molecules to examine the conformational motions of the protein.⁸⁵ The eigenvectors with corresponding eigenvalues of Cα atomic positions were extracted using “gmx covar”, and “gmx anaeig” was used to analyze the projections of the first two eigenvectors. The conformational dynamics of M^pro and M^pro complexes were scrutinized using the first two eigenvectors (PC1 and PC2) with the highest eigenvalues.⁸⁶ The free energy landscape (FEL)⁸⁷ was generated using the Boltzmann relationship with “gmx sham”:

G_(PC1, PC2) = −k_BT ln P_(PC1, PC2) (1)

where G refers to the Gibbs free energy, and T, k_B, and P represent the absolute temperature, Boltzmann constant, and probability distribution along PC1 and PC2, respectively. The Origin 9.1 package was used to plot the data obtained from all analyses.⁸⁸

3. Results and discussion

3.1. LBVS of small molecule databases to identify top-hit compounds against M^pro activity

Boceprevir was chosen as a reference compound for LBVS of small molecule databases to identify top-hit compounds against M^pro activity, as it exhibits potent inhibitory activity against M^pro.³² A total of 1535 small molecules with a score greater than 0.70 or a similarity index ≥70% with boceprevir have been shortlisted from the virtual screening of 15 small-molecule databases comprising ∼16 million compounds (Fig. 2). A total of 1137 compounds were selected for further screening studies after removing duplicates.

Fig. 2 Flowchart of the integrated computational methodology involving LBVS, SwissADME, molecular docking, MM-PBSA, and MD simulations employed in this work to identify potential M^pro inhibitors from various small-molecule databases.

3.2. Evaluation of drug-likeness properties and molecular docking of top-hit compounds with M^pro

The drug-likeness assessment of 1137 compounds obtained after removing duplicates was done using the online available SwissADME web server.⁵⁷ According to Lipinski's rule of five, compounds exhibiting two or more of the five parameters (hydrogen bond acceptors, hydrogen bond donors, logP, molecular weight, and TPSA) that lie outside the cutoff values have poor membrane permeability or oral bioavailability. A total of 850 compounds out of 1137 displayed no Lipinski violations. Molecular docking using AutoDock Vina⁶¹ was performed to evaluate the binding affinity and key interactions of 850 compounds with M^pro. The binding energy of screened compounds with M^pro evaluated using AutoDock Vina ranges from −9.9 to −5.1 kcal mol⁻¹ (Fig. 3(a)). Boceprevir binds to M^pro with a binding affinity of −7.5 kcal mol⁻¹, and it displays hydrogen bonds with Thr26 (3.0 Å), Asn142 (2.7 Å), Gly143 (2.4 Å), and Cys145 (2.9 Å) of M^pro (Fig. 3(b)). Furthermore, Thr25, Leu27, His41 (catalytic dyad residue), Cys44, Thr45, Met49, Phe140, Leu141, Ser144, His163, Met165, Glu166, and Gln189 of M^pro displayed hydrophobic contacts with boceprevir (Fig. 3(c)). The binding energy of boceprevir with M^pro is consistent with the binding energy for the M^pro–N3 complex (−7.5 kcal mol⁻¹) reported by Khamto et al.⁸⁹

Fig. 3 Distribution of 850 compounds as potential binders to M^pro over a range of docking energies (panel (a)). The docked pose of M^pro–boceprevir displaying the hydrogen bonds of boceprevir with the binding pocket residues of M^pro (panel (b)). The 2D interaction map of M^pro–boceprevir docked complex displaying hydrophobic contacts of boceprevir with M^pro residues (panel (c)). The docked poses display the binding of the top nine compounds and boceprevir in the S1, S2, S3, S4, and S5 subpockets of M^pro (panel (d)).

The top nine compounds were shortlisted from the library of 850 compounds based on their higher negative binding affinity and a large number of binding interactions with M^pro as compared to boceprevir (Fig. 1 and Fig. S1). The shortlisted compounds displayed the binding affinity with M^pro in the range of −9.9 to −8.0 kcal mol⁻¹ (Table S2). The docking analysis depicted the binding of all shortlisted compounds with the catalytic dyad (His41 and Cys145), including the other active site residues of M^pro. The docked poses of the top nine hits displayed hydrogen bond interactions (Fig. S2) and hydrophobic contacts (Fig. S3), preferably with the binding site of M^pro. Additionally, the top nine shortlisted compounds interacted inside the S1, S2, S3, S4, and S5 subpockets of the active site of M^pro (Fig. 3(d)).^20,21 The molecular docking results were validated by performing the redocking of screened hits with M^pro using Glide.⁶⁸ The key binding interactions of top hits with M^pro residues using Glide⁶⁸ were found to be approximately similar, as depicted in docked poses generated by AutoDock Vina (Table S2).

All nine top-hit compounds screened from molecular docking have an MW less than 500 Da and TPSA < 140 Ų in comparison to boceprevir (MW: 521.69 Da; TPSA: 153.86 Ų) (Table S4). All nine compounds and boceprevir possess HBA and HBD in the acceptable range as defined by Lipinski's rule (HBA ≤ 10 and HBD ≤ 5) and high lipophilicity as the partition coefficient lies below five. The low partition coefficient or high lipophilicity indicates the better membrane permeability of the selected compounds as compared to boceprevir. The results depicted no Lipinski violations in the top-hit compounds from the molecular docking, which were further screened using binding free energy analysis.

3.3. Estimation of binding free energy of top-hit compounds with M^pro

To obtain detailed insights into the key interactions between M^pro and the top nine compounds, the binding free energies of M^pro complexes were assessed using MM-PBSA (Table 2). The conformational entropy was excluded from the MM-PBSA calculations following the literature.⁷¹ The van der Waals (ΔE_vdW) interaction energy significantly contributed to the binding free energy of top hits with M^pro (Table 2), which is consistent with Khamto et al.⁸⁹ and Sanachai et al.,⁹⁰ highlighting that van der Waals interactions played an important role in the stability of M^pro–ligand complexes. The electrostatic (ΔE_elec) and non-polar solvation energy (ΔG_nps) components favoured the binding to top hits with M^pro, however, polar solvation energy (ΔG_ps) was noted to be unfavourable. The compounds C1, C3, C4, C5, C6, C8, and C9 displayed highly negative binding free energy in comparison to boceprevir (Table 2). Notably, C3, C5, and C9 exhibited significantly higher binding free energies of −65.2 ± 6.5, −66.1 ± 7.1, and −67.3 ± 5.8 kcal mol⁻¹, respectively, with M^pro as compared to M^pro–boceprevir (−51.8 ± 6.5 kcal mol⁻¹) (Table 2).

Table 2 Binding free energy of top-hit compounds with M^pro evaluated using MM-PBSA

Compounds	Binding free energy (kcal mol^−1)
	van der Waals (ΔEvdW)	Electrostatic (ΔEelec)	Molecular mechanics (ΔEMM^a)	Polar solvation (ΔGps)	Non-polar solvation (ΔGnps)	Solvation (ΔGsolv^b)	Binding energy (ΔGbinding^c)
a ΔEMM = ΔEvdW + ΔEelec. b ΔGsolv = ΔGps + ΔGnps. c ΔGbinding = ΔEMM + ΔGsolv.
C1	−33.3 ± 2.2	−8.3 ± 2.1	−41.6 ± 4.3	31.2 ± 3.5	−42.1 ± 5.2	−10.9 ± 1.7	−52.5 ± 6.0
C2	−30.0 ± 3.4	−7.3 ± 2.5	−37.3 ± 5.9	28.4 ± 5.4	−39.9 ± 6.5	−11.5 ± 1.1	−48.8 ± 7.0
C3	−42.6 ± 3.1	−14.7 ± 2.3	−57.3 ± 5.4	44.4 ± 4.4	−52.3 ± 5.5	−7.9 ± 1.1	−65.2 ± 6.5
C4	−41.5 ± 3.3	−10.2 ± 3.7	−51.7 ± 7.0	36.7 ± 5.2	−45.9 ± 5.0	−9.2 ± 0.2	−60.9 ± 7.2
C5	−45.3 ± 3.5	−11.5 ± 3.0	−56.8 ± 6.5	39.4 ± 4.6	−48.7 ± 5.2	−9.3 ± 0.6	−66.1 ± 7.1
C6	−39.9 ± 3.4	−11.1 ± 2.7	−51.0 ± 6.1	39.6 ± 8.8	−46.9 ± 7.8	−7.3 ± 1.0	−58.3 ± 7.1
C7	−33.5 ± 3.3	−5.7 ± 2.7	−39.2 ± 6.0	29.8 ± 3.9	−38.9 ± 4.9	−9.1 ± 1.0	−48.3 ± 7.0
C8	−41.4 ± 3.1	−9.0 ± 2.0	−50.4 ± 5.1	35.9 ± 3.0	−45.9 ± 5.3	−10.0 ± 2.3	−60.4 ± 7.4
C9	−46.9 ± 2.9	−13.2 ± 2.5	−60.1 ± 5.4	43.3 ± 4.3	−50.5 ± 4.7	−7.2 ± 0.4	−67.3 ± 5.8
Boceprevir	−32.1 ± 2.8	−4.1 ± 2.2	−36.2 ± 5.0	26.7 ± 4.6	−42.3 ± 6.1	−15.6 ± 1.5	−51.8 ± 6.5

---

The per-residue decomposition of binding free energy analysis depicted that a large number of residues of the binding pocket of M^pro were involved in the binding of C3, C5, and C9 with M^pro as compared to other compounds and boceprevir (Fig. 4). In M^pro–boceprevir, His41, Asn142, Gly143, Cys145, Met165, Glu166, and Pro168 of M^pro are preferably involved in the interactions with boceprevir (Fig. 4(a)). The M^pro residues Thr26, Thr27, His41 (catalytic dyad residue), Ser46, Asn142, Gly143, Ser144, Cys145 (catalytic dyad residue), Met165, Glu166, and Gln189 formed interactions with C3 (Fig. 4(b)). In the M^pro–C5 complex, C5 has interactions with Thr26, His41, Ser46, Asn142, Cys145, Met165, Glu166, and Pro168 of M^pro (Fig. 4(c)). The M^pro residues Thr26, His41, Ser46, Lys141, Asn142, Cys145, Met165, Glu166, and His172 had interactions with C9 (Fig. 4(d)). The residue-wise decomposition binding free energy analysis is consistent with that of Shao et al.,⁹¹ who determined notable contributions of Leu141, Gly143, Cys145, Met165, Glu166, Pro168, and Glu189 of M^pro in the M^pro–11a and M^pro–PF-07311332 complexes. The per-residue decomposition binding free energy analysis depicted interactions of C1, C2, C4, C6, and C8 with the active site residues of M^pro (Fig. S4). Overall, a large number of residues of the M^pro binding site displayed interactions with the top three shortlisted compounds (C3, C5, and C9).

Fig. 4 Per-residue decomposition of the binding free energies of M^pro with boceprevir, C3, C5, and C9 are shown in panels (a)–(d), respectively.

3.4. Impact of top-hit compounds on the structural stability of M^pro

The molecular docking analysis depicted that boceprevir and top hit compounds interacted with monomer A of the M^pro dimer [Fig. 3(b), (c) and Fig. S2, S3]. Furthermore, MD simulations were performed on the full dimeric form of M^pro in the presence of boceprevir and the top hit compounds to evaluate the binding stability of M^pro complexes. Notably, boceprevir and the top hit compounds remained bound to monomer A of M^pro during the simulation (Fig. 5). Thus, all the structural and dynamic analyses, including backbone RMSD, RMSF, R_g, and hydrogen bond assessments, were specifically focused on monomer A, where the ligands were bound, to interpret the ligand-induced conformational changes more precisely.

Fig. 5 Superimposed conformations of the snapshots at 10 ns intervals depicting the identical binding positions of boceprevir, C3, C5, and C9 to monomer A of M^pro are shown in panels (a–d), respectively.

The RMSD for backbone atoms of M^pro in the absence and presence of ligands has been evaluated (Fig. 6(a)). The apo–M^pro displayed a higher deviation with an average of 0.197 ± 0.007 nm during simulation, which coincides well with Handa et al.,⁹² where the average RMSD of 0.193 ± 0.038 nm was noted for the backbone atoms of apo–M^pro during the 100 ns simulation. However, the incorporation of ligands leads to significantly reduced RMSDs for M^pro complexes. M^pro–C3, M^pro–C5, and M^pro–C9 exhibit average RMSDs of 0.149 ± 0.003, 0.164 ± 0.008, and 0.158 ± 0.002 nm, respectively, which are comparatively lower than the average RMSD of M^pro–boceprevir (0.174 ± 0.009 nm), indicating the higher structural stability of M^pro in the presence of C3, C5, and C9. Furthermore, the simulations have been extended to 200 ns to assess the stability of M^pro complexes. The RMSD plot indicates that all the ligand-bound M^pro complexes (M^pro–C3, M^pro–C5, and M^pro–C9) reached equilibrium and maintained stable RMSD values during simulation (Fig. S5a), with notably reduced fluctuations compared to apo–M^pro. Similarly, the R_g profiles remained stable, indicating the compactness of the M^pro in the simulated systems over 200 ns (Fig. S5b). The ligand-bound complexes exhibited a slightly more compact and stable conformation as compared to apo–M^pro. To ensure the reproducibility of MD simulations, triplicate simulations of apo–M^pro and duplicates of M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 were performed using different initial velocities [Fig. S6(a)–(e)]. Remarkably, all simulations exhibited nearly identical RMSD profiles with similar RMSD values (Fig. S6f), which demonstrates the robust reproducibility of all MD simulations.

Fig. 6 RMSD (panel (a)) and R_g (panel (b)) of M^pro in the absence and presence of boceprevir, C3, C5, and C9. The domain regions of monomer A of M^pro are displayed in panel (c), where the loops L1 (green), L2 (orange), L3 (purple), L4 (magenta), and β-sheet regions β1 (red) and β2 (hot pink) are shown. Residue-wise RMSF for apo–M^pro and M^pro complexes are displayed in panel (d).

Furthermore, the R_g was analyzed to predict the compactness of M^pro with and without ligands. Initially, a relatively high R_g was noted in apo–M^pro up to 40 ns and then attained a stable profile for the last 60 ns. The M^pro complexes showed a stable R_g profile during the whole simulation as compared to apo–M^pro (Fig. 6(b)). The average R_g for apo–M^pro, M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 was found to be 2.241 ± 0.002, 2.214 ± 0.003, 2.203 ± 0.002, 2.213 ± 0.001, and 2.220 ± 0.000 nm, respectively. The R_g noted for apo–M^pro aligns well with Sen et al.,⁹³ where the average R_g for apo–M^pro was 2.251 nm during a 100 ns MD simulation. Samanta et al.³⁸ reported the average R_g for complexes of remdesivir analogues (58059494, 134133102, and 76314404 screened from the PubChem database) with M^pro varied from 2.200 to 2.300 nm, which corroborates well with the R_g of M^pro complexes. Also, the R_g results are consistent with the average R_g of apo–M^pro (2.213 nm), and M^pro in the presence of Hc1, Hc2, Hc3, Hc4, Cg1, and Cg2 (2.233, 2.245, 2.228, 2.240, 2.256, 2.204 nm, respectively) reported by Han et al.⁹⁴ It is noted that the average R_g of all complexes is marginally lower than that of apo–M^pro, however, only M^pro–C3 and M^pro–C5 exhibit a lower R_g compared to M^pro–boceprevir. The R_g results indicated that M^pro with C3, C5, and C9 exhibits higher compactness as compared to apo–M^pro, which depicted the stability of M^pro complexes.

The residue-wise RMSF for Cα atoms of M^pro and M^pro complexes was assessed to investigate the flexibility of M^pro (Fig. 6(d)). The RMSF plot displayed the lower fluctuations in M^pro complexes as compared to apo–M^pro in loops, L1 (Arg40–Phe66), L2 (Gly138–Ser147), L3 (Gly183–Asp197), and L4 (Arg279–Leu286), and β-sheet regions β1 (Val20–Leu30) and β2 (Cys160–Thr175) of M^pro. The loops L1, L2, L3, and L4, and β-sheet regions β1 and β2 are color coded in Fig. 6(c). The reduced flexibility in M^pro complexes is consistent with results of Li et al., indicating the notably reduced flexibility of the L1, L2, L3, L4, β1, and β2 regions of M^pro on the binding of YTV, YSP, and YU4 as compared to apo–M^pro.⁹⁵ The average RMSFs for apo–M^pro and in the presence of boceprevir, C3, C5, and C9 were noted to be 0.106 ± 0.005, 0.092 ± 0.003, 0.075 ± 0.001, 0.081 ± 0.003, and 0.087 ± 0.002 nm, respectively (Fig. 6(d)). M^pro–C3, M^pro–C5, and M^pro–C9 have reduced flexibility in 93.75%, 87.82%, and 82.23% residues of M^pro as compared M^pro–boceprevir (74.34%), which indicates the enhanced conformational stability of M^pro on the incorporation of C3, C5, and C9.

3.5. SASA and hydrogen bond analyses of apo–M^pro and M^pro complexes

The SASA analysis offers valuable insights into the stability of protein–ligand complexes, as a consistent SASA value shows that the ligand and protein remain bound to each other during simulation. The average SASA for M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 was noted to be 184.868 ± 0.638, 180.144 ± 0.188, 182.661 ± 0.665, and 184.213 ± 0.194 nm², respectively (Fig. 7(a)), which is lower than the average SASA for apo–M^pro (186.335 ± 0.162 nm²). The SASA plot displayed the marginally reduced SASA for M^pro–C3 following M^pro–C5 and M^pro–C9 complexes as compared to M^pro–boceprevir. The stable SASA for M^pro complexes suggests that ligands remained inside the subpockets of the M^pro active site during the entire simulation.

Fig. 7 SASA and the number of intramolecular hydrogen bonds of apo–M^pro and M^pro complexes are shown in panels (a) and (b), respectively. The per-residue hydrogen bond probability of M^pro with boceprevir, C3, C5, and C9 is displayed in panel (c).

The intramolecular hydrogen bonds have paramount importance and play a crucial role in maintaining the structural integrity and stability of a protein's native structure.⁹⁶ The number of intramolecular hydrogen bonds in M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 were noted to be 206.945 ± 0.690, 202.714 ± 0.675, 205.944 ± 1.401, and 204.555 ± 0.823, respectively, which are lower than that of apo–M^pro (213.714 ± 1.559) (Fig. 7(b)). The binding of C3, C5, and C9 with M^pro resulted in slightly reduced intramolecular hydrogen bonds as compared to M^pro–boceprevir and apo–M^pro. Furthermore, residue-wise hydrogen bond probability analysis indicated a significant contribution of Thr26, Asn142, and Gly143 of M^pro in the hydrogen bond interactions with the compounds, followed by Ser46, Met49, Cys145, Glu166, and Gln189 (Fig. 7(c)). The average number of intermolecular hydrogen bonds was noted to be 3.582 ± 0.352, 4.807 ± 0.225, 2.602 ± 0.222, and 2.701 ± 0.141 for M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9, respectively (Fig. 8), which depicts strong binding of boceprevir, C3, C5, and C9 with M^pro.

Fig. 8 Intermolecular hydrogen bonds of boceprevir, C3, C5, and C9 with M^pro are shown in panels (a), (c), (e), and (g), respectively. The representative members of the most-populated conformational clusters of M^pro complexes highlight the binding of boceprevir, C3, C5, and C9 with binding site residues of M^pro in panels (b), (d), (f), and (h), respectively.

Additionally, the hydrogen bonds between M^pro residues and top hits including boceprevir are noted in the representative conformers of most-populated conformational clusters (Fig. 8). Boceprevir displayed hydrogen bond interactions with oxyanion hole residues (Asn142, Gly143) and catalytic dyad residue, Cys145, with a distance of 3.0, 2.2, and 2.8 Å, respectively [Fig. 8(a) and (b)]. C3 interacted with Thr26 (2.5 Å), Asn142 (2.4 Å), Gly143 (1.9 Å), Cys145 (2.3 Å), and Glu166 (2.3 Å) of M^pro through hydrogen bonds [Fig. 8(c) and (d)] and blocked the active site of M^pro. This is consistent with the hydrogen bond occupancy analysis reported by Li et al.,³⁹ where Gly143 (<20%) and Glu166 of M^pro with 60% occupancy were involved in hydrogen bond interactions with PF-07321332. C5 formed hydrogen bonds with Thr26 (2.5 Å) and Asn142 (2.7 Å, 2.8 Å), whereas C9 displayed hydrogen bonds with Thr26 (2.1 Å, 3.0 Å) and Ser46 (2.3 Å) [Fig. 8(e)–(h)]. The hydrogen bond analysis highlighted that all compounds interacted with the active site residues of M^pro, however, C3 displayed a large number of hydrogen bond interactions in the active site as compared to boceprevir, C5, and C9.

The top hit compounds displayed hydrophobic contacts with the active site residues of M^pro as noted in the representative conformers of the most-populated conformational clusters (Fig. S7). Boceprevir displayed hydrophobic contacts with Thr25, Thr26, Leu27, His41, Ser46, Met49, Ser144, His164, Met165, Gln189, Thr190, and Ala191 of M^pro (Fig. S7a), whereas C3 showed hydrophobic contacts with Thr25, Leu27, His41, Cys44, Ser46, Met49, His164, Met165, and Gln189 (Fig. S7b); C5 with Thr25, Leu27, His41, Val42, Cys44, Thr45, Met49, Cys145, and Glu166 (Fig. S7c); and C9 with Thr25, His41, Thr45, Met49, Phe140–Asn142, Ser144, Cys145, Glu166, and Gln189 of M^pro (Fig. S7d).

Furthermore, the COM distances between the catalytic residues (His41 and Cys145) of M^pro and top hit compounds were evaluated during simulation to assess the binding interactions as well as the dynamic stability of the identified hits within the active site of M^pro (Fig. S8). Notably, the top hit compounds consistently maintained a distance of ∼7.0 Å from His41 and Cys145 during simulation, indicating favourable contacts of top hit compounds with M^pro. This is consistent with the 2D interaction maps of representative members of the most-populated conformational clusters of M^pro complexes, depicting hydrophobic contacts of top hit compounds with His41, along with notable contacts of C5 and C9 with Cys145 (Fig. S7). Furthermore, the average distances of boceprevir and C3 with Cys145 were noted to be 2.28 ± 0.06 and 2.37 ± 0.02 Å, respectively, during simulation (Fig. S8), consistent with representative members of the most-populated conformational clusters of M^pro complexes, depicting hydrogen bond interactions of boceprevir and C3 with Cys145 of M^pro (Fig. 8). The COM distances remained consistent during simulation, indicating that boceprevir and top hit compounds maintain a close proximity to His41 and Cys145, suggesting their stable interactions at the active site of M^pro.

3.6. Conformational snapshots depict the binding of top-hit compounds in the active site of M^pro

To gain further insights into the stability of the M^pro complexes, a visual inspection of all MD simulation trajectories was performed. Subsequently, several conformational snapshots from each M^pro complex were extracted at various time intervals during simulation (0, 25, 50, 75, and 100 ns) (Fig. S9). The conformational snapshots of M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 complexes demonstrate that the ligands remained tightly and consistently bound to the active site of M^pro during simulation [Fig. S9(a)–(d)]. In addition, conformations of boceprevir, C3, C5, and C9 at every 10 ns time interval were superimposed (Fig. 5), depicting the almost identical binding position of ligands during simulation.

3.7. PCA and FEL analyses depict reduced conformational motions and higher conformational stability of M^pro in the presence of top-hit compounds

PCA was conducted to distinguish the collective motions of M^pro with and without boceprevir, C3, C5, and C9 during MD simulations. It is noteworthy that the initial two eigenvectors represented higher collective motions in M^pro (42.892%) in the absence and incorporation of boceprevir (41.582%), C3 (43.057%), C5 (37.298%), and C9 (36.217%) (Fig. S10a). Thus, the first two eigenvectors (PC1 and PC2) accounting for significant amounts of overall motions were chosen to examine the conformational dynamics of apo–M^pro, M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 (Fig. S10a). The projection of MD trajectories on the first two principal components (PC1 and PC2) highlighted the lower conformational phase space in M^pro complexes as compared to apo–M^pro [Fig. S10(b)–(f)]. Also, M^pro–C3, M^pro–C5, and M^pro–C9 displayed lower conformational phase space as compared to M^pro–boceprevir [Fig. S10(c)–(f)], depicting the significantly reduced conformational motions in M^pro upon the incorporation of top hits.

Furthermore, the trace value of the covariance matrix was calculated to analyze the conformational flexibility across all systems. The M^pro complexes exhibited lower trace values of 3.659 nm² (M^pro–boceprevir), 3.108 nm² (M^pro–C3), 2.404 nm² (M^pro–C5), and 2.806 nm² (M^pro–C9) as compared to apo–M^pro (4.418 nm²). A lower trace value indicates the reduced conformational flexibility of M^pro on incorporating C3, C5, and C9 as compared to apo–M^pro and M^pro–boceprevir.

The displacements in Cα atoms of M^pro residues with and without boceprevir and three top hits along eigenvectors 1 and 2 were investigated to scrutinize the conformational flexibility of M^pro (Fig. 9). The displacements of M^pro residues in the presence of boceprevir and top three hits along eigenvector 1 exhibited lower fluctuations in loops, L1 (Arg40–Phe66), L2 (Gly138–Ser147), L3 (Gly183–Asp197), and L4 (Arg279–Leu286) and β-sheet regions β1 (Val20–Leu30) and β2 (Cys160–Thr175) of M^pro (Fig. 9(a)). Notably, lower fluctuations in L1, L2, L3, and L4 loops and β1 and β2 regions were also noted along eigenvector 2 in M^pro complexes as compared to apo–M^pro (Fig. 9(b)). The displacement analysis indicated that M^pro–C3, M^pro–C5, and M^pro–C9 have reduced flexibility in 85.52%, 79.28%, and 81.25%, residues of M^pro, respectively, along eigenvector 1, whereas 75.99%, 63.82%, and 71.38% residues of M^pro, respectively, displayed lower fluctuations along eigenvector 2 as compared to M^pro–boceprevir (lower fluctuations in 78.95% residues along eigenvector 1, and 67.11% along eigenvector 2) and apo–M^pro. The PCA depicts enhanced conformational stability of M^pro on the incorporation of C3, C5, and C9.

Fig. 9 Residue-wise displacements of Cα atoms of M^pro in the absence and presence of boceprevir, C3, C5, and C9, along eigenvectors 1 and 2 are shown in panels (a) and (b), respectively.

To examine the influence of boceprevir and the top hits on the structural robustness of M^pro, the first two principal components, PC1 and PC2, were used to generate the FELs (Fig. 10 and 11). The FELs highlighted that the energy varies between 0 and 3.2 kcal mol⁻¹ for apo–M^pro (Fig. 10(a)) and 0 and 3.0 kcal mol⁻¹ for M^pro–boceprevir (Fig. 10(b)), which is reduced on the binding of C3, C5, and C9 to M^pro [Fig. 11(a)–(c)]. The FEL plots indicated the three minimum energy basins with higher sampling of conformations in apo–M^pro, whereas two minimum energy basins with higher sampling in M^pro–boceprevir, M^pro–C9, and only one minimum energy basin with highly stable conformations in M^pro–C3 and M^pro–C5. The FEL analysis highlighted that the binding of C3 and C5 significantly modified the conformational space of M^pro and stabilized its structure.

Fig. 10 FEL of M^pro in the absence and presence of boceprevir with corresponding free energy conformations shown in panels (a) and (b), respectively. The inset views (S4, S5) portray the hydrogen bond interactions of boceprevir with M^pro.

Fig. 11 FEL plots of M^pro on the incorporation of C3, C5, and C9 with corresponding free energy conformations are shown in panels (a)–(c), respectively. The inset views of extracted conformations represent the hydrogen bond interactions of top hits with M^pro.

Additionally, the minimum-energy conformations from each free energy basin were extracted for M^pro, M^pro–boceprevir, M^pro–C3, M^pro–C5, and M^pro–C9 (Fig. 10 and 11). In apo–M^pro, three energy basins consisting of three distinct conformations (S1, S2, S3) illustrate that the conformational landscape of apo–M^pro is predominantly partitioned into three distinct subspaces, reflecting its structural heterogeneity (Fig. 10(a)). The two free energy conformations falling into energy basins S4 and S5 were extracted for M^pro–boceprevir (Fig. 10(b)), illustrating the hydrogen bond interactions of boceprevir with oxyanion hole residues, Asn142 (2.5 Å), Gly143 (3.0 Å), and catalytic dyad residue, Cys145 (2.2, 2.4 Å) of M^pro. However, only a single free energy basin with a higher sampling of conformations in M^pro–C3 (S1′) and M^pro–C5 (S2′) was observed [Fig. 11(a) and (b)], where C3 displayed hydrogen bonds with Asn142 (2.0 Å), Gly143 (2.6 Å), Cys145 (2.6 Å), and Glu166 (2.1 Å) and C5 formed hydrogen bonds with Thr26 (1.9, 3.0 Å) and Asn142 (1.7 Å) of M^pro. Three free energy conformations (S3′, S4′, S5′) corresponding to three distinct energy basins were retrieved from the FEL of M^pro–C9 (Fig. 11(c)). C9 participated in hydrogen bonds with Thr26 (1.8 Å) in S3′, Thr26 (1.8, 2.8 Å) and Asn142 (2.9 Å) in S4′, Thr26 (1.7, 2.6 Å) and Ser46 (2.4 Å) of M^pro in S5′ conformations. Overall, FEL analysis suggested that boceprevir, C3, C5, and C9 bind to the oxyanion hole and catalytic dyad residues of M^pro, which aligns well with the conformational clustering analysis. In addition to non-covalent interactions such as van der Waals interactions, hydrophobic contacts, and hydrogen bonding, the presence of electrophilic moieties (e.g., amides and carbonyls) in compounds C3, C5, and C9 suggests the potential of these compounds to bind covalently with the catalytic cysteine (Cys145) of M^pro. This mechanism is analogous to that of boceprevir, a covalent reversible inhibitor that forms a covalent bond with Cys145 through its electrophilic ketoamide warhead.^32,33c,34

The top hit compound C3 discovered as a potent inhibitor of M^pro activity in this work has been previously reported as an inhibitor (namely HC-toxin) of histone deacetylases (HDACs) in plant, yeast, and mammalian cells.⁹⁷ HC-toxin was derived from Cochliobolus (Helminthosporium) carbonum by Pringle and coworkers in 1971;⁹⁸ however, its possible structure was elucidated by Gross et al. in 1982.⁹⁹ Gross et al. reported that HC-toxin consists of cyclo(^DPro–^LAla–^DAla–^LAeo), where Aeo: 2-amino-9,10-epoxy-8-oxodecanoic acid.⁹⁹ Loidl and coworkers demonstrated that HC-toxin displayed >50% and >95% inhibitory activity against the maize HDAC at 2 and 20 μM, respectively.^97a In 2004, Sheen and coworkers reported that HC-toxin inhibited cell proliferation, triggers cell cycle arrest and apoptosis of human breast cancer cells (MCF-7 and MDA-MB-468).^46a The same group demonstrated potent antiproliferative efficacy, apoptosis, and cell cycle arrest activity of HC-toxin in T47D human breast cancer cells.^46b Later, Debuzer et al. discovered that HC-toxin at <20 nM effectively suppresses the malignant phenotype in established neuroblastoma (NB) cell lines as well as primary NB cells, triggers cell cycle arrest and apoptosis, promotes neuronal differentiation, and notably reduces colony formation and invasive cell growth.⁴⁷

Interestingly, C3 interacted with oxyanion hole residues (Asn142, Gly143) and catalytic dyad residue, Cys145, of SARS-CoV-2 M^pro, exhibiting a binding free energy of −65.2 ± 6.5 kcal mol⁻¹. Importantly, C3 reduced the residual fluctuations in M^pro leading to higher structural stability. Hydrogen bond analysis depicted that C3 displayed interactions with Thr26, oxyanion hole residues (Asn142 and Gly143), the catalytic residue (Cys145), and Glu166 of M^pro. Notably, C3 displayed hydrophobic contacts with active site residues (Thr25, Leu27, His41, Cys44, Ser46, Met49, His164, Met165, and Gln189) of M^pro. Overall, C3 demonstrated a strong potential as a therapeutic candidate against M^pro activity and has the ability to be repurposed as an antiviral agent against SARS-CoV M^pro.

4. Conclusions

In this work, the LBVS approach using boceprevir as a reference compound and MD simulations has been employed to discover new M^pro inhibitors from various small molecule databases. Notably, molecular docking followed by binding free energy evaluation using MM-PBSA identified ChEMBL144205 (C3), ZINC000091755358 (C5), and ZINC000092066113 (C9) as high-affinity binders of M^pro. Remarkably, the drug-likeness profile of top hit compounds displayed no Lipinski violations. Importantly, MD simulations depicted the high structural stability and lower conformational fluctuations in M^pro on the incorporation of C3, C5, and C9. The residue-specific binding free energy analysis depicted the binding interactions of top hit compounds (C3, C5, and C9) with the oxyanion hole residues (Asn142, Gly143), one of the catalytic dyad residues (i.e., Cys145), and Glu166 of M^pro that potentially inhibits the catalytic activity of M^pro, leading to an inactive state of M^pro. PCA displayed the lower flexibility in M^pro residues on the binding of C3, C5, and C9, and SASA depicted the stable binding of top hit compounds in the active site of M^pro during simulation. Notably, C3 depicted a large number of hydrogen bond interactions in the active site of M^pro as compared to boceprevir, C5, and C9, which, in turn, highlights its potential as a repurposed drug candidate against M^pro activity. C3 previously displayed inhibition against cell proliferation, triggers cell cycle arrest, and induces apoptosis of human breast cancer cells. Notably, previous studies highlighted that C3 at <20 nM effectively suppresses the malignant phenotype in established neuroblastoma (NB) cell lines as well as primary NB cells. Thus, the experimental validation and further optimization of the newly identified top hit compounds (C3, C5, and C9) in this work will yield more potent analogs of boceprevir as antiviral agents against SARS-CoV-2. In conclusion, the in silico data presented in this work will provide key mechanistic insights to drive medicinal chemistry projects for the discovery of more promising boceprevir analogs as M^pro inhibitors.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

The data related to this work are available within the article and the SI. Supplementary information: Chemical structures of the top hit compounds identified using LBVS, molecular docking and MM-PBSA from various small molecule databases, docked poses of the top hit compounds displaying hydrogen bonds with M^pro, 2D interaction maps of docked poses depicting hydrophobic contacts of the top hit compounds with M^pro, per-residue decomposition binding free energy plots of M^pro on the incorporation of C1, C2, C4, C6, C7, and C8, RMSD and R_g of M^pro in the absence and presence of boceprevir, C3, C5, and C9 during 200 ns MD simulations, RMSD of repeat simulations of apo–M^pro and top three M^pro complexes, 2D interaction maps of representative members of most-populated conformational clusters of M^pro complexes depicting hydrophobic contacts, variations in the COM distances of the active site residues (His41 and Cys145) of M^pro with boceprevir, C3, C5, and C9 during simulation, snapshots of M^pro complexes at various time intervals during simulation, and PCA plots depicting the conformational motions of M^pro in the absence and presence of boceprevir, C3, C5, and C9 are displayed in Fig. S1–S10. SMILES format of the top hit compounds and boceprevir, molecular docking analysis of top hit compounds with M^pro using AutoDock Vina and Glide, protonation states of M^pro residues at pH 7.4 predicted using PROPKA, and drug-likeness parameters of top hit compounds and boceprevir are listed in Tables S1–S4. The SI consists of PDBQT parameters of boceprevir and top hit compounds. See DOI: https://doi.org/10.1039/d5cp01814e

The library of small molecules was virtually screened by SwissSimilarity (http://www.swisssimilarity.ch/). 2D chemical structures of small molecules were generated using ChemDraw Professional 16.0 (https://perkinelmer-chemdraw-professional.software.informer.com/16.0/). The pdbqt coordinates of small molecules were obtained using Open Babel (https://openbabel.org). The M^pro structure was retrieved from the RCSB Protein Data Bank (PDB ID: 6Y84). Molecular docking was performed using AutoDock Vina 1.1 (https://vina.scripps.edu), available in the public domain, and Glide, available at (https://www.schrodinger.com/platform/products/glide/). The MD simulations were performed using the open-source package GROMACS 2022.4 (https://www.gromacs.org). Parameters for OPLS–AA/L force field for reference and screened compounds were generated using LigParGen Server (https://zarbi.chem.yale.edu/ligpargen/). The binding free energy of top hit compounds, screened from various small molecule databases using ligand-based virtual screening and molecular docking, with M^pro was evaluated using g_mmpbsa script (https://rashmikumari.github.io/g_mmpbsa/). Drug likeness parameters of the top hit compounds were evaluated using SwissADME, freely available at http://www.swissadme.ch/.. Origin 9.1 was used for data plotting. Binding interactions between M^pro and the top hit compounds were visualized using Ligplot+ (https://www.ebi.ac.uk/thornton-srv/software/LigPlus/) and PyMOL (https://www.pymol.org/), available for academic use. Microsoft PowerPoint was used to generate the figures (https://www.microsoft.com/en-in/microsoft-365/powerpoint). The input data, protein–ligand complexes, parameter and topology files, MD output files, and binding free energy files are available for download at https://github.com/Gurmeet-kaur06/SARS-CoV-2-Mpro-Screening.git.

Acknowledgements

BG acknowledges SERB, New Delhi, Government of India (CRG/2023/000088 and CRG/2022/008244) and CSIR, Government of India [02(0451)/21/EMR-II] for the research funding. GK thanks the CSIR, Government of India, for the Direct Senior Research Fellowship (09/0677(18151)/2024-EMR-I). The authors acknowledge the Department of Chemistry & Biochemistry, TIET Patiala, India for the research infrastructure.

References

1. https://covid19.who.int/ (Accessed on August 04, 2025).
2. B. Hu, H. Guo, P. Zhou and Z.-L. Shi, Characteristics of SARS-CoV-2 and COVID-19, Nat. Rev. Microbiol., 2021, 19, 141–154.
3. X. Zhang, Y. Zhang, W. Qiao, J. Zhang and Z. Qi, Baricitinib, a drug with potential effect to prevent SARS-CoV-2 from entering target cells and control cytokine storm induced by COVID-19, Int. Immunopharmacol., 2020, 86, 106749.
4. S. Joshi, J. Parkar, A. Ansari, A. Vora, D. Talwar, M. Tiwaskar, S. Patil and H. Barkate, Role of favipiravir in the treatment of COVID-19, Int. J. Infect. Dis., 2021, 102, 501–508.
5. P. Gautret, J.-C. Lagier, P. Parola, V. T. Hoang, L. Meddeb, M. Mailhe, B. Doudier, J. Courjon, V. Giordanengo, V. E. Vieira, H. T. Dupont, S. Honoré, P. Colson, E. Chabrière, B. La Scola, J.-M. Rolain, P. Brouqui and D. Raoult, Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial, Int. J. Antimicrob. Agents, 2020, 56, 105949.
6. A. J. Bernal, M. M. G. da Silva, D. B. Musungaie, E. Kovalchuk, A. Gonzalez, V. D. Reyes, A. Martín-Quirós, Y. Caraco, A. Williams-Diaz, M. L. Brown, J. Du, A. Pedley, C. Assaid, J. Strizki, J. A. Grobler, H. H. Shamsuddin, R. Tipping, H. Wan, A. Paschke, J. R. Butterton, M. G. Johnson and C. De Anda, Molnupiravir for oral treatment of COVID-19 in non-hospitalized patients, N. Engl. J. Med., 2022, 386, 509–520.
7. E. Mahase, COVID-19: Pfizer's paxlovid is 89% effective in patients at risk of serious illness, company reports, BMJ, 2021, 375, n2713.
8. A. Frediansyah, F. Nainu, K. Dhama, M. Mudatsir and H. Harapan, Remdesivir and its antiviral activity against COVID-19: a systematic review, Clin. Epidemiol. Global Health, 2021, 9, 123–127.
9. İ. Aygün, M. Kaya and R. Alhajj, Identifying side effects of commonly used drugs in the treatment of Covid 19, Sci. Rep., 2020, 10, 21508.
10. S. Su, G. Wong, W. Shi, J. Liu, A. C. K. Lai, J. Zhou, W. Liu, Y. Bi and G. F. Gao, Epidemiology, genetic recombination, and pathogenesis of coronaviruses, Trends Microbiol., 2016, 24, 490–502.
11. C.-C. Lai, T.-P. Shih, W.-C. Ko, H.-J. Tang and P.-R. Hsueh, Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): the epidemic and the challenges, Int. J. Antimicrob. Agents, 2020, 55, 105924.
12. R. Cannalire, C. Cerchia, A. R. Beccari, F. S. Di Leva and V. Summa, Targeting SARS-CoV-2 proteases and polymerase for COVID-19 treatment: state of the art and future opportunities, J. Med. Chem., 2022, 65, 2716–2746.
13. C. Gil, T. Ginex, I. Maestro, V. Nozal, L. Barrado-Gil, M. Á. Cuesta-Geijo, J. Urquiza, D. Ramírez, C. Alonso, N. E. Campillo and A. Martinez, COVID-19: drug targets and potential treatments, J. Med. Chem., 2020, 63, 12359–12386.
14. S. Ullrich and C. Nitsche, The SARS-CoV-2 main protease as drug target, Bioorg. Med. Chem. Lett., 2020, 30, 127377.
15. (a) J. S. Morse, T. Lalonde, S. Xu and W. R. Liu, Learning from the past: possible urgent prevention and treatment options for severe acute respiratory infections caused by 2019-nCoV, ChemBioChem, 2020, 21, 730–738; (b) Z. Jin, X. Du, Y. Xu, Y. Deng, M. Liu, Y. Zhao, B. Zhang, X. Li, L. Zhang, C. Peng, Y. Duan, J. Yu, L. Wang, K. Yang, F. Liu, R. Jiang, X. Yang, T. You, X. Liu, X. Yang, F. Bai, H. Liu, X. Liu, L. W. Guddat, W. Xu, G. Xiao, C. Qin, Z. Shi, H. Jiang, Z. Rao and H. Yang, Structure of M^pro from SARS-CoV-2 and discovery of its inhibitors, Nature, 2020, 582, 289–293.
16. F. A. Ambrosio, G. Costa, I. Romeo, F. Esposito, M. Alkhatib, R. Salpini, V. Svicher, A. Corona, P. Malune, E. Tramontano, F. Ceccherini-Silberstein, S. Alcaro and A. Artese, Targeting SARS-CoV-2 main protease: a successful story guided by an in silico drug repurposing approach, J. Chem. Inf. Model., 2023, 63, 3601–3613.
17. T. Pillaiyar, M. Manickam, V. Namasivayam, Y. Hayashi and S.-H. Jung, An overview of severe acute respiratory syndrome–coronavirus (SARS-CoV) 3CL protease inhibitors: peptidomimetics and small molecule chemotherapy, J. Med. Chem., 2016, 59, 6595–6628.
18. (a) L. Zhang, D. Lin, X. Sun, U. Curth, C. Drosten, L. Sauerhering, S. Becker, K. Rox and R. Hilgenfeld, Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α–ketoamide inhibitors, Science, 2020, 368, 409–412; (b) B. Goyal and D. Goyal, Targeting the dimerization of the main protease of coronaviruses: a potential broad–spectrum therapeutic strategy, ACS Comb. Sci., 2020, 22, 297–305.
19. J. C. Ferreira, S. Fadl, A. J. Villanueva and W. M. Rabeh, Catalytic dyad residues His41 and Cys145 impact the catalytic activity and overall conformational fold of the main SARS-CoV-2 protease 3-chymotrypsin-like protease, Front. Chem., 2021, 9, 692168.
20. A. Gahlawat, N. Kumar, R. Kumar, H. Sandhu, I. Singh, S. Singh, A. Sjöstedt and P. Garg, Structure-based virtual screening to discover potential lead molecules for the SARS-CoV-2 main protease, J. Chem. Inf. Model., 2020, 60, 5781–5793.
21. G. Kaur and B. Goyal, Insights into the interaction mechanism of boceprevir with SARS-CoV-2 main protease, ChemistrySelect, 2023, 8, e202301415.
22. H. Chen, P. Wei, C. Huang, L. Tan, Y. Liu and L. Lai, Only one protomer is active in the dimer of SARS 3C-like proteinase, J. Biol. Chem., 2006, 281, 13894–13898.
23. G. G. Chang, Quaternary structure of the SARS coronavirus main protease, in Molecular biology of the SARS-coronavirus, ed. S. Lal, Springer, Berlin, Heidelberg, 2010, pp. 115–128.
24. Q. Hu, Y. Xiong, G.-H. Zhu, Y.-N. Zhang, Y.-W. Zhang, P. Huang and G.-B. Ge, The SARS-CoV-2 main protease (M^pro): structure, function, and emerging therapies for COVID-19, MedComm, 2022, 3, e151.
25. Y. Zhang, L. Zheng, Y. Yang, Y. Qu, Y.-Q. Li, M. Zhao, Y. Mu and W. Li, Structural and energetic features of the dimerization of the main proteinase of SARS-CoV-2 using molecular dynamics simulations, Phys. Chem. Chem. Phys., 2022, 24, 4324–4333.
26. S. Lida and Y. Fukunishi, Asymmetric dynamics of dimeric SARS-CoV-2 and SARS-CoV main proteases in an apo form: molecular dynamics study on fluctuations of active site, catalytic dyad, and hydration water, BBA Adv., 2021, 1, 100016.
27. Y. Yang, Y.-D. Luo, C.-B. Zhang, Y. Xiang, X.-Y. Bai, D. Zhang, Z.-Y. Fu, R.-B. Hao and X.-L. Liu, Progress in research on the inhibitors targeting SARS-CoV-2 main protease (M^pro), ACS Omega, 2024, 9, 34196–34219.
28. A. C. Puhl, A. S. Godoy, G. D. Noske, A. M. Nakamura, V. O. Gawriljuk, R. S. Fernandes, G. Oliva and S. Ekins, Discovery of PL^pro and M^pro inhibitors for SARS-CoV-2, ACS Omega, 2023, 8, 22603–22612.
29. F. B. Omage, A. Madabeni, A. R. Tucci, P. A. Nogara, M. Bortoli, A. D. S. Rosa, V. N. D. S. Ferreira, J. B. T. Rocha, M. D. Miranda and L. Orian, Diphenyl diselenide and SARS-CoV-2: in silico exploration of the mechanisms of inhibition of main protease (M^pro) and papain-like protease (PL^pro), J. Chem. Inf. Model., 2023, 63, 2226–2239.
30. H. Liu, S. Iketani, A. Zask, N. Khanizeman, E. Bednarova, F. Forouhar, B. Fowler, S. J. Hong, H. Mohri, M. S. Nair, Y. Huang, N. E. S. Tay, S. Lee, C. Karan, S. J. Resnick, C. Quinn, W. Li, H. Shion, X. Xia, J. D. Daniels, M. Bartolo-Cruz, M. Farina, P. Rajbhandari, C. Jurtschenko, M. A. Lauber, T. McDonald, M. E. Stokes, B. L. Hurst, T. Rovis, A. Chavez, D. D. Ho and B. R. Stockwell, Development of optimized drug-like small molecule inhibitors of the SARS-CoV-2 3CL protease for treatment of COVID-19, Nat. Commun., 2022, 13, 1891.
31. M. Elsbaey, M. A. A. Ibrahim, A. M. Shawky and T. Miyamoto, Eryngium creticum L.: chemical characterization, SARS-CoV-2 inhibitory activity, and in silico study, ACS Omega, 2022, 7, 22725–22734.
32. C. Ma, M. D. Sacco, B. Hurst, J. A. Townsend, Y. Hu, T. Szeto, X. Zhang, B. Tarbet, M. T. Marty, Y. Chen and J. Wang, Boceprevir, GC-376, and calpain inhibitors II, XII inhibit SARS-CoV-2 viral replication by targeting the viral main protease, Cell Res., 2020, 30, 678–692.
33. (a) V. Minicozzi, A. Giuliani, G. Mei, L. Domenichelli, M. Parise, A. D. Venere and L. D. Paola, The dynamical asymmetry in SARS-CoV2 protease reveals the exchange between catalytic activity and stability in homodimers, Molecules, 2025, 30, 1412; (b) M. Arooj, I. Shehadi, C. N. Nassab and A. A. Mohamed, Computational insights into binding mechanism of drugs as potential inhibitors against SARS-CoV-2 targets, Chem. Pap., 2022, 76, 111–121; (c) S. Das, A. Singh, S. K. Samanta and A. S. Roy, Naturally occurring anthraquinones as potential inhibitors of SARS-CoV-2 main protease: an integrated computational study, Biologia, 2022, 77, 1121–1134; (d) A. Manandhar, B. E. Blass, D. J. Colussi, I. Almi, M. Abou-Gharbia, M. L. Klein and K. M. Elokely, Targeting SARS-CoV-2 M3CL^pro by HCV NS3/4a inhibitors: in silico modeling and in vitro screening, J. Chem. Inf. Model., 2021, 61, 1020–1032; (e) S. Borkotoky, M. Banerjee, G. P. Modi and V. K. Dubey, Identification of high affinity and low molecular alternatives of boceprevir against SARS-CoV-2 main protease: a virtual screening approach, Chem. Phys. Lett., 2021, 770, 138446.
34. (a) B. Tan, R. Joyce, H. Tan, Y. Hu and J. Wang, SARS-CoV-2 main protease drug design, assay development, and drug resistance studies, Acc. Chem. Res., 2023, 56, 157–168; (b) D. W. Kneller, H. Li, G. Phillips, K. L. Weiss, Q. Zhang, M. A. Arnould, C. B. Jonsson, S. Surendranathan, J. Parvathareddy, M. P. Blakeley, L. Coates, J. M. Louis, P. V. Bonnesen and A. Kovalevsky, Covalent narlaprevir– and boceprevir–derived hybrid inhibitors of SARS-CoV-2 main protease, Nat. Commun., 2022, 13, 2268; (c) M. Göhl, L. Zhang, H. E. Kilani, X. Sun, K. Zhang, M. Brönstrup and R. Hilgenfeld, From repurposing to redesign: optimization of boceprevir to highly potent inhibitors of the SARS-CoV-2 main protease, Molecules, 2022, 27, 4292; (d) Y. R. Alugubelli, Z. Z. Geng, K. S. Yang, N. Shaabani, K. Khatua, X. R. Ma, E. C. Vatansever, C.-C. Cho, Y. Ma, J. Xiao, L. R. Blankenship, G. Yu, B. Sankaran, P. Li, R. Allen, H. Ji, S. Xu and W. R. Liu, A systematic exploration of boceprevir-based main protease inhibitors as SARS-CoV-2 antivirals, Eur. J. Med. Chem., 2022, 240, 114596; (e) R. Oerleman, A. J. Ruiz-Moreno, Y. Cong, N. D. Kumar, M. A. Velasco-Velazquez, C. G. Neochoritis, J. Smith, F. Reggiori, M. R. Groves and A. Dömling, Repurposing the HCV NS3–4A protease drug boceprevir as COVID-19 therapeutics, RSC Med. Chem., 2021, 12, 370–379; (f) L. Fu, F. Ye, Y. Feng, F. Yu, Q. Wang, Y. Wu, C. Zhao, H. Sun, B. Huang, P. Niu, H. Song, Y. Shi, X. Li, W. Tan, J. Qi and G. F. Gao, Both boceprevir and GC376 efficaciously inhibit SARS-CoV-2 by targeting its main protease, Nat. Commun., 2020, 11, 4417.
35. A. Singh, K. Jangid, S. Nehul, P. Dhaka, R. Rani, A. Pareek, G. K. Sharma, P. Kumar and S. Tomar, Structural and mechanistic insights into the main protease (Mpro) dimer interface destabilization inhibitor: unveiling new therapeutic avenues against SARS-CoV-2, Biochemistry, 2025, 64, 1589–1605.
36. N. Khamto, K. Utama, P. Boontawee, A. Janthong, S. Tatieng, S. Arthan, V. Choommongkol, P. Sangthong, C. Yenjai, N. Suree and P. Meepowpan, Inhibitory activity of flavonoid scaffolds on SARS-CoV-2 3CL^pro: insights from the computational and experimental investigations, J. Chem. Inf. Model., 2024, 64, 874–891.
37. X. Jin, M. Zhang, B. Fu, M. Li, J. Yang, Z. Zhang, C. Li, H. Zhang, H. Wu, W. Xue and Y. Liu, Structure-based discovery of the SARS-CoV-2 main protease noncovalent inhibitors from traditional Chinese medicine, J. Chem. Inf. Model., 2024, 64, 1319–1330.
38. P. N. Samanta, D. Majumdar and J. Leszczynski, Elucidating atomistic insight into the dynamical responses of the SARS-CoV-2 main protease for the binding of remdesivir analogues: leveraging molecular mechanics to decode the inhibition mechanism, J. Chem. Inf. Model., 2023, 63, 3404–3422.
39. X. Li, Z. Fang, D. Li and Z. Li, Binding kinetics study of SARS-CoV-2 main protease and potential inhibitors via molecular dynamics simulations, Phys. Chem. Chem. Phys., 2023, 25, 15135–15145.
40. J. Hammond, H. Leister-Tebbe, A. Gardner, P. Abreu, W. Bao, W. Wisemandle, M. Baniecki, V. M. Hendrick, B. Damle, A. Simon-Campos, R. Pypstra and J. M. Rusnak, Oral nirmatrelvir for high-risk, nonhospitalized adults with Covid-19, N. Engl. J. Med., 2022, 386, 1397–1408.
41. D. R. Owen, C. M. N. Allerton, A. S. Anderson, L. Aschenbrenner, M. Avery, S. Berritt, B. Boras, R. D. Cardin, A. Carlo, K. J. Coffman, A. Dantonio, L. Di, H. Eng, R. Ferre, K. S. Gajiwala, S. A. Gibson, S. E. Greasley, B. L. Hurst, E. P. Kadar, A. S. Kalgutkar, J. C. Lee, J. Lee, W. Liu, S. W. Mason, S. Noell, J. J. Novak, R. S. Obach, K. Ogilvie, N. C. Patel, M. Pettersson, D. K. Rai, M. R. Reese, M. F. Sammons, J. G. Sathish, R. S. P. Singh, C. M. Steppan, A. E. Stewart, J. B. Tuttle, L. Updyke, P. R. Verhoest, L. Wei, Q. Yang and Y. Zhu, An oral SARS-CoV-2 M^pro inhibitor clinical candidate for the treatment of COVID-19, Science, 2021, 374, 1586–1593.
42. I. Polatoğlu, T. Oncu-Oner, I. Dalman and S. Ozdogan, COVID-19 in early 2023: structure, replication mechanism, variants of SARS-CoV-2, diagnostic tests, and vaccine & drug development studies, MedComm, 2023, 4, e228.
43. A. A. Al-Karmalawy, R. Soltane, A. A. Elmaaty, M. A. Tantawy, S. A. Antar, G. Yahya, A. Chrouda, R. A. Pashameah, M. Mustafa, M. A. Mraheil and A. Mostafa, Coronavirus disease (COVID-19) control between drug repurposing and vaccination: a comprehensive overview, Vaccines, 2021, 9, 1317.
44. N. A. Ashour, A. A. Elmaaty, A. A. Sarhan, E. B. Elkaeed, A. M. Moussa, I. A. Erfan and A. A. Al-Karmalawy, A systematic review of the global intervention for SARS-CoV-2 combating: from drugs repurposing to molnupiravir approval, Drug Des., Dev. Ther., 2022, 16, 685–715.
45. (a) Z. Tanoli, A. Fernández-Torras, U. O. Özcan, A. Kushnir, K. M. Nader, Y. Gadiya, L. Fiorenza, A. Ianevski, M. Vähä-Koskela, M. Miihkinen, U. Seemab, H. Leinonen, B. Seashore-Ludlow, M. Tampere, A. Kalman, F. Ballante, E. Benfenati, G. Saunders, S. Potdar, I. G. García, R. García-Serna, C. Talarico, A. R. Beccari, W. Schaal, A. Polo, S. Costantini, E. Cabri, M. Jacobs, J. Saarela, A. Budillon, O. Spjuth, P. Östling, H. Xhaard, J. Quintana, J. Mestres, P. Gribbon, A. E. Ussi, D. C. Lo, M. D. Kort, K. Wennerberg, M. Fratelli, J. Carreras-Puigvert and T. Aittokallio, Computational drug repurposing: approaches, evaluation of in silico resources and case studies, Nat. Rev. Drug Discovery, 2025, 24, 521–542; (b) J. L. Cummings, Y. Zhou, A. V. Stone, D. Cammann, R. Tonegawa-Kuji, J. Fonseca and F. Cheng, Drug repurposing for Alzheimer's disease and other neurodegenerative disorders, Nat. Commun., 2025, 16, 1755; (c) M. Salvadè, M. DiLuca and F. Gardoni, An update on drug repurposing in Parkinson's disease: preclinical and clinical considerations, Biomed. Pharmacother., 2025, 183, 117862; (d) N. Wankhade, A. Sharma, M. A. Wani, A. Banerjee and P. Garg, Predictive modeling and drug repurposing for type-II diabetes, ACS Med. Chem. Lett., 2024, 15, 1907–1917.
46. (a) K. N. Min, M. J. Cho, D.-K. Kim and Y. Y. Sheen, Estrogen receptor enhances the antiproliferative effects of Trichostatin A and HC-toxin in human breast cancer cells, Arch. Pharm. Res., 2004, 27, 554–561; (b) K. E. Joung, D.-K. Kim and Y. Y. Sheen, Antiproliferative effect of Trichostatin A and HC-toxin in T47D human breast cancer cells, Arch. Pharm. Res., 2004, 27, 640–645.
47. H. E. Deubzer, V. Ehemann, F. Westermann, R. Heinrich, G. Mechtersheimer, A. E. Kulozik, M. Schwab and O. Witt, Histone deacetylase inhibitor Helminthosporium carbonum (HC)-toxin suppresses the malignant phenotype of neuroblastoma cells, Int. J. Cancer, 2008, 122, 1891–1900.
48. (a) M. E. Bragina, A. Daina, M. A. S. Perez, O. Michielin and V. Zoete, SwissSimilarity 2021 web tool: novel chemical libraries and additional methods for an enhanced ligand-based virtual screening experience, Int. J. Mol. Sci., 2022, 23, 811; (b) V. Zoete, A. Daina, C. Bovigny and O. Michielin, SwissSimilarity: a web tool for low to ultra high throughput ligand-based virtual screening, J. Chem. Inf. Model., 2016, 56, 1399–1404.
49. https://www.chembridge.com .
50. A. Gaulton, L. J. Bellis, A. P. Bento, J. Chambers, M. Davies, A. Hersey, Y. Light, S. McGlinchey, D. Michalovich, B. Al-Lazikani and J. P. Overington, ChEMBL: a large-scale bioactivity database for drug discovery, Nucleic Acids Res., 2012, 40, D1100–D1107.
51. K. Degtyarenko, P. de Matos, M. Ennis, J. Hastings, M. Zbinden, A. McNaught, R. Alcántara, M. Darsow, M. Guedj and M. Ashburner, ChEBI: a database and ontology for chemical entities of biological interest, Nucleic Acids Res., 2008, 36, D344–D350.
52. D. S. Wishart, C. Knox, A. C. Guo, S. Shrivastava, M. Hassanali, P. Stothard, Z. Chang and J. Woolsey, DrugBank: a comprehensive resource for in silico drug discovery and exploration, Nucleic Acids Res., 2006, 34, D668–D672.
53. Z. Feng, L. Chen, H. Maddula, O. Akcan, R. Oughtred, H. M. Berman and J. Westbrook, Ligand Depot: a data warehouse for ligands bound to macromolecules, Bioinformatics, 2004, 20, 2153–2155.
54. W. Chan, H. Zhang, J. Yang, J. Brender, J. Hur, A. Ozgur and Y. Zhang, GLASS: a comprehensive database for experimentally-validated GPCR-ligand associations, Bioinformatics, 2015, 31, 3035–3042.
55. D. S. Wishart, D. Tzur, C. Knox, R. Eisner, A. C. Guo, N. Young, D. Cheng, K. Jewell, D. Arndt, S. Sawhney, C. Fung, L. Nikolai, M. Lewis, M. A. Coutouly, I. Forsythe, P. Tang, S. Shrivastava, K. Jeroncic, P. Stothard, G. Amegbey, D. Block, D. D. Hau, J. Wagner, J. Miniaci, M. Clements, M. Gebremedhin, N. Guo, Y. Zhang, G. E. Duggan, G. D. Macinnis, A. M. Weljie, R. Dowlatabadi, F. Bamforth, D. Clive, R. Greiner, L. Li, T. Marrie, B. D. Sykes, H. J. Vogel and L. Querengesser, HMDB: the Human metabolome database, Nucleic Acids Res., 2007, 35, D521–D526.
56. J. J. Irwin and B. K. Shoichet, ZINC–a free database of commercially available compounds for virtual screening, J. Chem. Inf. Model., 2005, 45, 177–182.
57. A. Daina, O. Michielin and V. Zoete, SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules, Sci. Rep., 2017, 7, 42717.
58. H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov and P. E. Bourne, The protein data bank, Nucleic Acids Res., 2000, 28, 235–242.
59. N. O'Boyle, M. Banck, C. A. James, C. Morley, T. Vandermeersch and G. R. Hutchison, Open Babel: an open chemical toolbox, J. Cheminf., 2011, 3, 33.
60. (a) T. A. Halgren, Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94, J. Comput. Chem., 1996, 17, 490–519; (b) T. A. Halgren, Merck molecular force field. II. MMFF94 van der Waals and electrostatic parameters for intermolecular interactions, J. Comput. Chem., 1996, 17, 520–552; (c) T. A. Halgren, Merck molecular force field. III. Molecular geometries and vibrational frequencies for MMFF94, J. Comput. Chem., 1996, 17, 553–586; (d) T. A. Halgren and R. Nachbar, Merck molecular force field. IV. Conformational energies and geometries for MMFF94, J. Comput. Chem., 1996, 17, 587–615; (e) T. A. Halgren, Merck molecular force field. V. Extension of MMFF94 using experimental data, additional computational data, and empirical rules, J. Comput. Chem., 1996, 17, 616–641.
61. O. Trott and A. J. Olson, AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading, J. Comput. Chem., 2010, 31, 455–461.
62. G. M. Morris, D. S. Goodsell, R. S. Halliday, R. Huey, W. E. Hart, R. K. Belew and A. J. Olson, Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function, J. Comput. Chem., 1998, 19, 1639–1662.
63. R. Huey, G. M. Morris, A. J. Olson and D. S. A. Goodsell, Semiempirical free energy force field with charge-based desolvation, J. Comput. Chem., 2007, 28, 1145–1152.
64. F. J. Solis and R. J.-B. Wets, Minimization by random search techniques, Math. Oper. Res., 1981, 6, 19–30.
65. G. M. Morris, R. Huey, W. Lindstrom, M. F. Sanner, R. K. Belew, D. S. Goodsell and A. J. Olson, AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility, J. Comput. Chem., 2009, 16, 2785–2791.
66. W. L. DeLano, The PyMOL molecular graphics system, DeLano Scientific, San Carlos, CA, 2002, vol. 571.
67. R. A. Laskowski and M. B. Swindells, LigPlot+: multiple ligand-protein interaction diagrams for drug discovery, J. Chem. Inf. Model., 2011, 51, 2778–2786.
68. R. A. Friesner, R. B. Murphy, M. P. Repasky, L. L. Frye, J. R. Greenwood, T. A. Halgren, P. C. Sanschagrin and D. T. Mainz, Extra precision Glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes, J. Med. Chem., 2006, 49, 6177–6196.
69. R. Kumari and R. Kumar, Open Source Drug Discovery Consortium and A. Lynn, g_mmpbsa—A GROMACS tool for high-throughput MM-PBSA calculations, J. Chem. Inf. Model., 2014, 54, 1951–1962.
70. (a) N. Forouzesh, F. Ghafouri, I. S. Tolokh and A. V. Onufriev, Optimal dielectric boundary for binding free energy estimates in the implicit solvent, J. Chem. Inf. Model., 2024, 64, 9433–9448; (b) N. Forouzesh and N. Mishra, An effective MM/GBSA protocol for absolute binding free energy calculations: a case study on SARS-CoV-2 spike protein and the human ACE2 receptor, Molecules, 2021, 26, 2383; (c) P.-C. Su, C.-C. Tsai, S. Mehboob, K. E. Hevener and M. E. Johnson, Comparison of radii sets, entropy, QM methods, and sampling on MM-PBSA, MM-GBSA, and QM/MM-GBSA ligand binding energies of F. tularensis enoyl-ACP reductase (FabI), J. Comput. Chem., 2015, 36, 1859–1873.
71. (a) T. Zhang, W. Xu, Y. Mu and P. Derreumaux, Atomic and dynamic insights into the beneficial effect of the 1,4–naphthoquinon–2–yl–l–tryptophan inhibitor on Alzheimer's Aβ1–42 dimer in terms of aggregation and toxicity, ACS Chem. Neurosci., 2014, 5, 148–159; (b) R. B. Best, N.-V. Buchete and G. Hummer, Are current molecular dynamics force fields too helical?, Biophys. J., 2008, 95, L07–L09.
72. (a) M. J. Abraham, T. Murtola, R. Schulz, S. Pall, J. C. Smith, B. Hess and E. Lindahl, GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25; (b) D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, GROMACS: fast, flexible, and free, J. Comput. Chem., 2005, 26, 1701–1718.
73. (a) G. A. Kaminski, R. A. Friesner, J. Tirado-Rives and W. L. Jorgensen, Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides, J. Phys. Chem. B, 2001, 105, 6474–6487; (b) W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids, J. Am. Chem. Soc., 1996, 118, 11225–11236.
74. (a) L. S. Dodda, J. Z. Vilseck, J. Tirado-Rives and W. L. Jorgensen, 1.14*CM1A-LBCC: localized bond-charge corrected CM1A charges for condensed-phase simulations, J. Phys. Chem. B, 2017, 121, 3864–3870; (b) L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands, Nucleic Acids Res., 2017, 45, W331–W336; (c) W. L. Jorgensen and J. Tirado-Rives, Potential energy functions for atomic-level simulations of water and organic and biomolecular systems, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6665–6670.
75. P. Mark and L. Nilsson, Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K, J. Phys. Chem. A, 2001, 105, 9954–9960.
76. (a) C. R. Søndergaard, M. H. M. Olsson, M. Rostkowski and J. H. Jensen, Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pKa Values, J. Chem. Theory Comput., 2011, 7, 2284–2295; (b) H. Li, A. D. Robertson and J. H. Jensen, Very fast empirical prediction and rationalization of protein pKa values, Proteins: Struct., Funct., Bioinf., 2005, 61, 704–721.
77. B. Nutho, P. Mahalapbutr, K. Hengphasatporn, N. C. Pattaranggoon, N. Simanon, Y. Shigeta, S. Hannongbua and T. Rungrotmongkol, Why are lopinavir and ritonavir effective against the newly emerged coronavirus 2019? atomistic insights into the inhibitory mechanisms, Biochemistry, 2020, 59, 1769–1779.
78. G. Bussi, D. Donadio and M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys., 2007, 126, 014101.
79. M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys., 1981, 52, 7182–7190.
80. (a) D. J. Evans and B. L. Holian, The Nose–Hoover thermostat, J. Chem. Phys., 1985, 83, 4069–4074; (b) S. A. Nosé, Unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519.
81. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, LINCS: a linear constraint solver for molecular simulations, J. Comput. Chem., 1997, 18, 1463–1472.
82. S. Miyamoto and P. A. Kollman, Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models, J. Comput. Chem., 1992, 13, 952–962.
83. (a) U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. A. Pedersen, Smooth particle mesh Ewald method, J. Chem. Phys., 1995, 103, 8577–8593; (b) T. Darden, D. York and L. Pedersen, Particle Mesh Ewald: an N·log (N) method for Ewald sums in large systems, J. Chem. Phys., 1993, 98, 10089–10092.
84. (a) L. J. Smith, X. Daura and W. F. van Gunsteren, Assessing equilibration and convergence in biomolecular simulations, Proteins: Struct., Funct., Genet., 2002, 48, 487–496; (b) X. Daura, K. Gademann, B. Jaun, D. Seebach, W. F. van Gunsteren and A. E. Mark, Peptide folding: when simulation meets experiment, Angew. Chem., Int. Ed., 1999, 38, 236–240.
85. (a) A. Amadei, A. B. M. Linssen and H. J. C. Berendsen, Essential dynamics of proteins, Proteins: Struct., Funct., Genet., 1993, 17, 412–425; (b) T. Ichiye and M. Karplus, Collective motions in proteins: a covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations, Proteins: Struct., Funct., Genet., 1991, 11, 205–217.
86. (a) G. Kaur and B. Goyal, Deciphering the molecular mechanism of inhibition of β-secretase (BACE1) activity by a 2-amino-imidazol-4-one derivative, ChemistrySelect, 2022, 7, e202202561; (b) S.-Q. Liu, Z.-H. Meng, Y.-X. Fu and K.-Q. Zhang, Insights derived from molecular dynamics simulation into the molecular motions of serine protease proteinase K, J. Mol. Model., 2010, 16, 17–28.
87. (a) E. Papaleo, P. Mereghetti, P. Fantucci, R. Grandori and L. De Gioia, Free-energy landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: the myoglobin case, J. Mol. Graphics Modell., 2009, 27, 889–899; (b) A. Altis, M. Otten, P. H. Nguyen, R. Hegger and G. Stock, Construction of the free energy landscape of biomolecules via dihedral angle principal component analysis, J. Chem. Phys., 2008, 128, 245102.
88. E. Seifert, OriginPro 9.1: scientific data analysis and graphing software–software review, J. Chem. Inf. Model., 2014, 54, 1552.
89. N. Khamto, K. Utama, S. Tateing, P. Sangthong, P. Rithchumpon, N. Cheechana, A. Saiai, N. Semakul, W. Punyodom and P. Meepowpan, Discovery of natural bisbenzylisoquinoline analogs from the library of Thai traditional plants as SARS-CoV-2 3CL^pro inhibitors: in silico molecular docking, molecular dynamics, and in vitro enzymatic activity, J. Chem. Inf. Model., 2023, 63, 2104–2121.
90. K. Sanachai, T. Somboon, P. Wilasluck, P. Deetanya, P. Wolschann, T. Langer, V. S. Lee, K. Wangkanont, T. Rungrotmongkol and S. Hannongbua, Identification of repurposing therapeutics toward SARS-CoV-2 main protease by virtual screening, PLoS One, 2022, 17, e0269563.
91. H. P. Shao, T. H. Wang, H. L. Zhai, K. X. Bi and B. Q. Zhao, Discovery of inhibitors against SARS-CoV-2 main protease using fragment-base drug design, Chem.-Biol. Interact., 2023, 371, 110352.
92. Y. Handa, K. Okuwaki, Y. Kawashima, R. Hatada, Y. Mochizuki, Y. Komeiji, S. Tanaka, T. Furuishi, E. Yonemochi, T. Honma and K. Fukuzawa, Prediction of binding pose and affinity of Nelfinavir, a SARS-CoV-2 main protease repositioned drug, by combining docking, molecular dynamics, and fragment molecular orbital calculations, J. Phys. Chem. B, 2024, 128, 2249–2265.
93. D. Sen, B. Debnath, P. Debnath, S. Debnath, M. E. A. Zaki and V. H. Masand, Identification of potential edible mushroom as SARS-CoV-2 main protease inhibitor using rational drug designing approach, Sci. Rep., 2022, 12, 1503.
94. D. Han, H. Yang, Y. Gao, Y. Xue, F. Liu, M. Wang, J. Lu, T. Liu and Y. Xu, Virtual screening and molecular simulation uncover potent traditional Chinese medicine small-molecules against SARS-CoV-2 M^pro, ChemistrySelect, 2025, 10, e202405037.
95. M. Li, X. Liu, S. Zhang, S. Liang, Q. Zhang and J. Chen, Deciphering the binding mechanism of inhibitors of the SARS-CoV-2 main protease through multiple replica accelerated molecular dynamics simulations and free energy landscapes, Phys. Chem. Chem. Phys., 2022, 24, 22129–22143.
96. (a) C. N. Pace, H. Fu, K. L. Fryar, J. Landua, S. R. Trevino, D. Schell, R. L. Thurlkill, S. Imura, J. M. Scholtz, K. Gajiwala, J. Sevcik, L. Urbanikova, J. K. Myers, K. Takano, E. J. Hebert, B. A. Shirley and G. R. Grimsley, Contribution of hydrogen bonds to protein stability, Protein Sci., 2014, 23, 652–661; (b) N. H. Joh, A. Min, S. Fahman, J. P. Whitelegge, D. Yang, V. L. Woods and J. U. Bowie, Modest stabilization by most hydrogen-bonded side-chain interactions in membrane proteins, Nature, 2008, 453, 1266–1270.
97. (a) G. Brosch, R. Ransom, T. Lechner, J. D. Walton and P. Loidl, Inhibition of maize histone deacetylases by HC toxin, the host-selective toxin of Cochliobolus carbonum, Plant Cell, 1995, 7, 1941–1950; (b) J. D. Walton, HC–toxin, Phytochemistry, 2006, 67, 1406–1413; (c) C. Hildmann, D. Wegener, D. Riester, R. Hempel, A. Schober, J. Merana, L. Giurato, S. Guccione, T. K. Nielsen, R. Ficner and A. Schwienhorst, Substrate and inhibitor specificity of class 1 and class 2 histone deacetylases, J. Biotechnol., 2006, 124, 258–270.
98. R. B. Pringle, Amino acid composition of the host-specific toxin of Helminthosporium carbonum, Plant Physiol., 1971, 48, 756–759.
99. M. L. Gross, D. McCrery, F. Crow, K. B. Tomer, M. R. Pope, L. M. Ciuffetti, H. W. Knoche, J. M. Daly and L. D. Dunkle, The structure of the toxin from Helminthosporium carbonum, Tetrahedron Lett., 1982, 23, 5381–5384.

---