The SARS-CoV-2 3CLpro protein found in the PDB network database (https://www.rcsb.org) with the following PDB ID: 8FY7 and resolution: 1.94 Å was utilized in this research (Ngo et al., 2023). The high-resolution structure of 1.94 Å ensured precise and reliable spatial information, which was crucial for accurate molecular docking and dynamics simulations. The protein was subjected to pre-processing and optimization utilizing the Protein Preparation Wizard available in Schrödinger Release 2023-4. Initially, the spatial configuration was examined for any structural anomalies. Afterwards, hydrogen atoms were incorporated via a hydrogen bond assignment procedure. After eliminating any water or solvent molecules, the protein was optimized using the OPLS4 force field, while keeping all other settings at their default values. The energy-minimized structure was then utilized to create a grid file for ligand docking. The grid box was established by utilizing the native ligand within the 3CLpro as a point of reference, and its dimensions were adjusted to encompass an optimal calculation range during receptor grid setup.

Topscience database (https://www.tsbiochem.com/service/topscience-database) stands as a leading international chemical database dedicated to supporting the scientific community with a comprehensive collection of active compounds, facilitating drug discovery research. And it is a publicly accessible database of more than 10 million compounds. After an initial screening of over 10 million compounds from the Topscience database using the Lipinski’s rule, the effects of the remaining 5 million compounds on 3CLpro were further examined. The Create Phase Database module, provided by Schrödinger Release 2023-4, was utilized to prepare the ligands by applying specific filters and generating a customized ligand library, which included refining the 3D structures by adding hydrogen atoms, generating stereoisomers, and identifying ionization states. Further ligand refinement was achieved through energy minimization and optimization of the low-energy 3D structure, using the OPLS_2005 force field implemented in PHASE.

The stepwise molecular docking process was carried out using the Glide module in Schrödinger to identify the most promising hit compounds. The process began with high throughput virtual screening (HTVS), after which the ligands with the greatest glide scores (top 10%) were chosen for additional analysis. These selected ligands underwent standard precision (SP) docking. The highest-ranking 10% of ligands from the SP docking results were further subjected to extra precision (XP) docking to enhance accuracy and identify the most promising candidates. The Prime module was used to employ the molecular mechanics/generalized born surface area (MM/GBSA) method in order to evaluate the free binding energy between the protein and the docked ligands, which yielded an assessment of the affinities of ligands. The VSGB solvation model and the OPLS4 force field were used for the free binding energy calculations. Additionally, the strain energy of the compounds was calculated using the Strain Energy module. Finally, the most promising compounds were selected based on a combination of their docking scores, calculated binding free energies, and strain energies. The comprehensive evaluation ensured that the top candidates exhibited strong docking interactions, favourable binding affinities, and acceptable strain energies.

The absorption distribution metabolism excretion toxicity (ADMET) properties of all candidate compounds were analysed by using SwissADME (Daina et al., 2017) (http://www.swissadme. ch/) and ProTox 3.0 (Banerjee et al., 2018; Banerjee et al., 2024) (https://tox.charite.de/protox3/) web facilities. SwissADME was utilized to predict various pharmacokinetic properties such as absorption, distribution, metabolism, and excretion profiles of the compounds. ProTox 3.0 which integrated various models to predict the oral toxicity, hepatotoxicity, cytotoxicity, and potential adverse effects of compounds was employed to assess the toxicity profiles of candidates.

The 44 candidate compounds were purchased from Topscience (Topscience, Shanghai, China). DMSO (Innochem, Beijing, China) was utilized as a solvent for dissolving the compounds, and all solutions were stored at a temperature of −20°C. The screening kit for the SARS-CoV2 3CLpro inhibitor was acquired from the Beyotime Institute of Biotechnology, which was located in Shanghai, China. Fluorescence resonance energy transfer (FRET) method was used to detect the activity of 3CLpro in the novel coronavirus 3CLpro inhibitor screening kit. The 96 wells of black plate were used for the reactions that were carried out. Initially, 91 μL of 3CLpro assay reagent and 5 μL of the test substance were consecutively introduced into each sample well, with DMSO replacing compounds in the model wells. For the control wells, 5 μL of DMSO and 91 μL of assay buffer were introduced. To achieve comprehensive homogenization of the reaction mixture, the 96-well plate was agitated for a duration of 1 min. After that, a rapid addition of 4 μL of the substrate was made to every well, and it was completely mixed. The plate was placed in a dark environment and kept at a temperature of 37°C for a duration of 15–20 min. The measurements of fluorescence were carried out with the assistance of a multifunction enzyme labeling reader (SpectraMax M5, Molecular Devices) that had an excitation wavelength of 325 nm and an emission wavelength of 393 nm. The inhibition rates were determined at single points, and IC₅₀ values were plotted using Origin 2024 (n = 6). The inhibition rate was calculated as follows: I n h i b i t i o n   r a t e   100 % = RFU enzyme − RFU sample RFU enzyme − RFU control × 100 %

The GROMACS 2020.7 beta suite was utilized for all-atom MD simulations to analyse the dynamic binding behaviour and stability of protein-candidate compounds complexes. The simulations operated over a 200 ns duration at a temperature of 298.15 K by utilizing the Amber ff14SB force field (Maier et al., 2015). Parameters for the ligands were generated by using the general amber force field (GAFF) parameters (Wang et al., 2004) the fitted charges were calculated at the B3LYP/cc-pVDZ level by using Gaussian 16 (Frisch et al., 2016), Multiwfn (Lu and Chen, 2012) and Sobtop (Lu, 2022) software. For each system, solvation was performed using a 10 Å cubic box filled with TIP3P water molecules. The genion module in GROMACS was employed to introduce the necessary counter ions (Na⁺ and Cl⁻) to neutralize the system. The energy minimization process was conducted with periodic boundary conditions and utilized the steepest descent algorithm over 5,000 steps. Electrostatic interactions were computed using the Particle Mesh Ewald (PME) method. Using the V-rescale temperature coupling approach under the NVT ensemble, the system temperature was progressively raised from 0 K to 298.15 K. Afterwards, the NPT ensemble’s Parrinello-Rahman barostat was used to keep the system pressure at 1 atm. Finally, each system was subjected to a 200 ns molecular dynamics simulation to ensure that the systems reached equilibrium and stability.

The DFT is used to ascertain the electron’s density and energy properties. Gaussian 16 program was employed to perform analysis of compounds with visualization conducted using GaussView 6.0. The structural coordinates of compounds were optimized at the B3LYP/6–311++G (d,p) level of theory without imposing any symmetry constraints. The molecules have one multiplicity and zero charges. The electrostatic surface potential and HOMO-LUMO energy levels of compounds were derived from the optimized geometry utilizing Multiwfn.

Molecular docking is a crucial and essential technology that provides a valuable alternative to in vitro screening in the field of drug discovery (Stanzione et al., 2021). Schrödinger’s software, particularly the Glide module, was chosen for molecular docking due to its high accuracy, comprehensive suite of integrated tools, user-friendly interface, and advanced scoring functions. Before the molecular docking of 3CLpro, docking validation was conducted to ensure the reliability of the process, with the results available in the (Supplementary Figure S1). The molecular docking of 3CLpro was performed on 5 million compounds obtained by preliminary screening using Lipinski’s rule from the Topscience database. These candidate compounds were identified through a rigorous selection process. The initial screening considered the XP GScore, a metric that predicts the binding affinity of compounds to the target protein. Compounds with the most favourable XP GScore were then subjected to further analysis using MM/GBSA binding free energy calculations to provide a more accurate estimation of their binding affinities. Strain energy was also evaluated to ensure that the binding of these compounds did not induce significant conformational changes in the target protein, which could affect the stability of the protein-ligand complex (Supplementary Table S1). Lower strain energy was generally favorable as it suggested less conformational distortion upon binding.

Figure 4 provided a visual representation of key metrics used in the selection process. Figure 4A showed box plots of the distribution of XP GScore, MM/GBSA binding free energy, and strain energy for the top 44 compounds, with S-217622 highlighted. The compound S-217622 was chosen as a reference due to its noncovalent, nonpeptidic nature, broad-spectrum antiviral activity against SARS-CoV-2 variants, favorable pharmacokinetic profile and efficacy in animal models. These box plots illustrated that most of the top 44 compounds had more favorable XP GScore and MM/GBSA binding free energy compared to S-217622, indicating stronger binding affinities and stable binding. The XP GScore ranged from −10.13 to −4.72 kcal/mol among the compounds, with compound 15 having the lowest score (−10.13 kcal/mol), suggesting the strongest binding affinity. And S-217622 had a relatively moderate XP GScore of −6.58 kcal/mol. The ∆ G bind values ranged from −76.38 kcal/mol to −50.32 kcal/mol, with compound 11 showing the most stable binding (−76.38 kcal/mol). S-217622 had a ∆ G bind of −57.03 kcal/mol, which was relatively moderate compared to other compounds. The strain energy for most compounds was less than 10 kcal/mol, indicating minimal conformational changes. Besides, S-217622 had a strain energy of 5.300 kcal/mol, which was within an acceptable range but not the lowest (Supplementary Table S1).

The evaluation of properties is pivotal in pinpointing molecules that hold the greatest promise as effective drugs for particular diseases during the drug design and development phases (Zhou et al., 2016). These evaluations focus on determining the compounds’ capabilities for absorption, their bodily distribution and mobility, their metabolism and excretion potential, and the assessment of any related toxicological implications (Jana et al., 2018). In addition to the inhibitory efficacy, a comprehensive evaluation of the ADMET properties of the lead compounds was conducted to ensure their potential as drug candidates (Supplementary Tables S2, S3). Figure 4B summarized the ADMET properties and toxicity profiles for top 44 compounds compared to S-217622. Most compounds exhibited high absorption and moderate water solubility, with a balanced distribution regarding BBB penetration and P-gp interaction. The majority fell into toxicity class IV. In terms of toxicity, while many compounds were inactive for hepatotoxicity, respiratory toxicity, and carcinogenicity, a significant number showed activity for neurotoxicity, immunotoxicity, and mutagenicity. S-217622 displayed high absorption, moderate water solubility, no BBB penetration, P-gp interaction, and fell into toxicity class IV, with specific activities in neurotoxicity, respiratory toxicity, immunotoxicity, and mutagenicity. These findings highlighted the strengths and potential safety concerns of the top 44 compounds as 3CLpro inhibitors. After comprehensive analysis of multiple parameters, we selected and purchased 44 compounds for wet screening.

Although this study successfully identified potential SARS-CoV-2 3CLpro inhibitors using traditional computer-aided drug design (CADD) methods, the role of artificial intelligence (AI) in accelerating the drug discovery process has gained significant attention in recent years. AI technologies can rapidly process large-scale data, automate compound screening, and predict compound activity, selectivity, and ADMET properties, offering a time- and resource-efficient alternative to traditional methods. Future studies could explore how AI techniques can be integrated with CADD methods, particularly to expedite the discovery of SARS-CoV-2 3CLpro inhibitors, enhancing the efficiency of identifying potential drug candidates and reducing the cost of experimental screening.

The top 44 compounds were evaluated and experimentally analysed to confirm their inhibitory action against 3CLpro. The inhibitory activity of 3CLpro was assessed using the FRET technique. There was no difference between the sequence of amino acids found in 3CLpro and that of the natural new coronavirus 3CLpro. The detection mechanism involves linking the fluorescent donor (MCA) and fluorescent receptor (Dnp) to both ends of the natural substrate of the 3CLpro protease. Specifically, this substrate is MCA-avlqsgfr-lys (Dnp)-Lys-NH2. After 3CLpro cuts the substrate, fluorescence of MCA can be detected. The initial screens were conducted with inhibitor doses of 50 μM (Supplementary Table S4). Among the tested compounds, compound 1 exhibited a notable inhibition percentage of 85.4%, closely followed by compounds 43 (84.8%), 5 (84.0%), 36 (81.9%), and 34 (78.3%). Compound 10 also demonstrated significant inhibitory activity at 78.0%. Notably, the reference compound S-217622 displayed the highest inhibition at 99.2%, serving as a benchmark for evaluating the efficacy of the tested compounds. Overall, the screening identified several potent inhibitors of 3CLpro, particularly compounds 1, 5, 10, 34, 36, and 43 which demonstrated inhibition percentages above 70%. These results highlighted the potential of these compounds as lead candidates for further optimization and development of effective antiviral therapies against SARS-CoV-2.

To gain further insights into the potency of these compounds, we conducted concentration-response assays to determine their IC₅₀ values. Figure 5A illustrated the structure and IC₅₀ values for selected compounds against 3CLpro activity. Figure 5B presented the inhibition curves for compounds 1, 5, 10, 34, 36, 43, and the reference compound S-217622, showing their effectiveness across a range of concentrations. Among these, compound 36 demonstrated the highest potency with an IC₅₀ value of 4.47 ± 0.39 μ M, followed by compound 34 (6.12 ± 0.42 μ M). The inhibition curve for the reference compound S-217622 displayed a significantly lower IC₅₀ value of 0.047 ± 0.0032 μ M, which was consistent with the activity data reported in the literature. These results suggested that while the tested compounds exhibit promising inhibitory effects on SARS-CoV-2 3CLpro, particularly compounds 34 and 36, they are less potent than the reference inhibitor S-217622. Further optimization and structural modifications may enhance their efficacy, providing a foundation for developing effective antiviral therapies.

The all-atom MD simulations further examined the effects of explicit solvent molecules on the protein, focusing on fluctuations and conformational adjustments. Additionally, the simulations enabled the collection of time-averaged attributes of the complex system over different timescales (Wichapong et al., 2013). From the evaluation of comprehensive calculation and pharmacological results, compounds 34 and 36 were identified as outstanding candidates. Consequently, ligand-free 3CLpro (apo-3CLpro), and the docked complexes of 3CLpro with compounds 34 and 36 were selected for all-atom MD simulation studies, with an aim to evaluate the stability and dynamic behaviour of compounds 34 and 36 when complexed with the 3CLpro, providing deeper insights into their potential as effective inhibitors. To ensure the robustness and accuracy of the results, three independent MD simulations were conducted for each system. The RMSD values across all three trajectories were nearly identical, confirming the consistency and reliability of the simulations.

The energies of the HOMO and LUMO boundary orbitals play an essential part in the field of quantum chemistry as they affect the behaviour of a molecule reaction to other compounds (Guerroudj et al., 2021). Moreover, the boundary orbit makes it easy to characterize the chemical reactivity as well as the kinetic stability of analysed molecule (Hao et al., 2020). The electronic properties and molecular interactions of compounds 34 and 36 were thoroughly examined through frontier molecular orbital analysis and molecular electrostatic potential surfaces, as depicted in Figure 10. Compound 34 exhibited a HOMO energy of −6.141 eV and a LUMO energy of −2.895 eV, resulting in a gap of 3.246 eV. In contrast, compound 36 had a HOMO energy of −5.543 eV and a LUMO energy of −0.549 eV, producing a much larger gap of 4.994 eV. The smaller HOMO-LUMO gap in compound 34 suggested higher reactivity, which could correlate with a greater ability to participate in charge transfer interactions or other electron-related mechanisms with the biological target. Meanwhile, the larger gap in compound 36 indicated greater stability and potentially lower reactivity.

The molecular electrostatic potential maps provide further insights into the charge distribution across the molecular surface. Both compounds exhibit distinct patterns of positive and negative regions, marked by red (negative) and blue (positive) zones. For compound 34, the molecular electrostatic potential map indicated a significant region of negative potential localized near electronegative atoms, which suggested potential for forming strong hydrogen bonds or electrostatic interactions with positively charged residues in the protein binding site. On the other hand, compound 36 showed a more balanced distribution of positive and negative potential, with the central part of the molecule bearing considerable electron density. These analyses provided valuable insights into the reactivity, stability, and potential binding affinities of the compounds with their respective biological targets.