In this modeling study, we used the 3CL^pro (PDB code: 7LME) from SARS-CoV-2 in complex with ML300 as the crystal structure. We selected this crystal structure according to the following criteria (Tian et al., 2022): 1) The organism of the selected crystal structure should be SARS-CoV-2 rather than other species. The organism of the 3CLpro structure is SARS-CoV-2. 2) One of the quality indexes of protein crystal structure is resolution, which represents the uncertainty of atomic position in crystal structure model. When there are many crystal structures available, we choose the one with high resolution (that is, the one with small resolution value). Generally, structures with a resolution less than 3 Å are sufficient for pharmacophore modeling. The 3CL^pro crystal structure has a high resolution of 2.10 Å 3) The selected crystal structure should include an active pocket; the crystal structure of 3CL^pro contains the active-binding pocket. This structure file was prepared by the MOE program (Zhou et al., 2019a), with the following protocol: hydrogens were added, water molecules were removed, partial charges were computed and energy minimization was carried out using the Amber10:EHT forcefield. In the Pharmacophore Query Editor, PCH pharmacophore scheme including hydrogen-bond acceptor, aromatic, and hydrophobic feature was selected. Amber10:EHT forcefield was assigned to the system. Pharmacophore Query Editor was used to manually construct a visualized 3D-pharmacophore model by analyzing protein-ligand interaction in the 3CL^pro complex (Zhou et al., 2019a; Zhou et al., 2019b). In our study, the generated pharmacophore model was exported and translated into a pharmacophore file by a script.

According to a previously reported method (Zhou et al., 2019a), the GH method was used to evaluate precision of model selectivity. In the study, we used the generated pharmacophore model to successfully screen actives from a 1,090 molecular database. This database consists of the decoy set containing 1,080 molecules (retrieved from DUD-E database) (Mysinger et al., 2012) and 10 known 3CL^pro inhibitors (Han et al., 2022). The resulted mapping data was used to evaluate pharmacophore quality by solving the following equation (The GH score with greater than 0.6 indicated a good model): G H = ( H a ( 3 A + H t ) 4 H t A ) ( 1 − ( H t − H a ) ( D − A ) )

The virtual screening was performed to screen a chemical database containing 35,000 molecules (Zhou et al., 2019a). In the course of screening, we applied the Pharmacophore Search as a screening protocol to identify active compounds that have good matching with query features of the generated pharmacophore model (Zhou et al., 2019a). Hit molecules can be ranked by their RMSDx values which indicate how well the matching ligand annotation points of the chemical structures were mapped onto the query features of the model. A low RMSDx value indicates a good matching of the ligand annotation points with query features of the model.

In this study, we used the 3CL^pro structure from SARS-CoV-2 with ML300 (PDB code: 7LME) from the pharmacophore study (Han et al., 2022). It has high resolution and relatively complete structure, so it is used as the docking model. By using default settings, the hit compounds obtained were docked into the 3CL^pro active site containing some amino acid residues (such as Phe140, Leu141, Met165, Met49, Thr25, and Thr24) through the Triangle Matcher Docking protocol of MOE (Zhou et al., 2019a; Zhou et al., 2019b). Amber10:EHT forcefield was assigned to the system. The dG docking scoring function are selected to rank the compounds.

ADMETlab web server (https://admetmesh.scbdd.com/) was used to predict the ADME properties of selected hits (Tian et al., 2022). The molecular weight (mol_MW), number of hydrogen bond acceptors (nHA), number of hydrogen bond donors (nHD), log of the octanol/water partition coefficient (logP), and log of the aqueous solubility (LogS) were evaluated.

According to a previously reported method (Tian et al., 2022), protein-hit complexes were investigated by molecular dynamics simulation using Groningen machine for chemical simulations software (GROMACS).

The purified 3CL^pro protein was purchased from Abcam (Cambridge, MA, United States). According to a previously reported method (Hang et al., 2019), the binding affinity of the compounds with SARS-CoV-2 3CL^pro was detected by MST assay. Briefly, purified SARS-CoV-2 3CL^pro was labelled with the Monolith NT Protein Labelling Kit RED (NanoTemper Technologies). Serially diluted compounds, with concentrations of 0.76 nM–25 μM, were mixed with 100 nM labelled SARS-CoV-2 3CL^pro at room temperature and loaded into Monolith standard-treated capillaries. The fluorescent signal was detected by Monolith NT.115 (NanoTemper, Munich, Germany). The K _d value was calculated by fitting a standard binding curve. Experiment was performed in triplicate.

According to a previously reported method (Han et al., 2022), 200 nM of SARS-CoV-2 3CL^pro were incubated with different concentrations of compounds (0, 0.19, 0.76, 3.05, 12.21, 48.83, 195.31, 781.25, 3125, 12,500, and 50,000 nM) for 20 min. The reaction was initiated by adding fluorophore-quencher peptide substrate HiyteFluor-488ESATLQSGLRKAK-(QXL)-NH₂, followed by 30 min incubation at 25°C. Fluorescence intensity was measured on a microplate reader (λex = 485 nm; λem = 528 nm). Experiment was performed in triplicate.

Recently, researchers have reported the structural complex of SARS-CoV-2 3CL^pro with ML300 (PDB code: 7LME) (Han et al., 2022). Therefore, we used this structure analyze the chemical features of 3CL^pro-ML300 interaction (Figure 1A). Based on their interaction between 3CL^pro and ML300 (Figure 1B), the complex structure was used to generate a 3CL^pro-pharmacophore model (Figure 1C). This model consists of four features: F1 hydrophobic feature, F2 hydrogen-bond acceptor feature, F3 aromatic feature, and F4 aromatic feature. The F1 feature mapped with its aromatic ring of ML300 described hydrophobic interactions with residues Phe140 and Leu141 in the 3CL^pro-active site (Figure 1D). The F2 feature mapped with the oxygen atom of ML300, formed an important hydrogen-bond interaction between ML300 and Glu166, while the aromatic rings of ML300 mapped on the aromatic features are involved in interactions with hydrophobic residues Met165 and Met49. The results indicates that the 3CL^pro-model features could effectively map the interactions between ligands and 3CL^pro-active residues.

The quality of 3CL^pro-model was validated using GH scoring method (Tian et al., 2022). The GH analysis were done by computing statistical parameters such as the enrichment factor (E) and goodness of hit score (GH). The 3CL^pro-model was successful in retrieving 90% of active compounds from the decoys set (Table 1). Moreover, an enrichment factor of 89 and a GH score of 0.84 indicated good quality of the model. The result reveals that the model can distinguish active inhibitors of 3CL^pro from the decoy set.

The flowchart in Figure 2 is a detailed schematic diagram of our virtual screening process in this study. The validated SARS-CoV-2-3CL^pro-model is applied as a 3D search query to virtually screen a 35,000-compound database. The obtained 922 compounds that could be mapped onto 3CL^pro-model model were successfully retrieved, including 103 compounds with RMSDx values lower than 0.65 Å (the lower the value, the better the pharmacophore mapping). Subsequently, the retrieved 103 compounds were docked into the 3CL^pro-active site to reduce the false positive molecules. The 3CL^pro inhibitor ML300 with a docking score of -10.90 kcal/mol was used as a positive control. Finally, the top four hits (hits 1–4) with lower docking scores (lower than -10.90 kcal/mol) were selected (Table 2). The analysis of the model matching results of four hit compounds showed that the oxygen atoms of four hit compounds matched the F2 hydrogen-bond acceptor feature of 3CL^pro-model while their two aromatic rings were mapped with the F3 and F4 aromatic features, respectively (Figure 3). Moreover, the hydrophobic groups of four hits matched the features of F1 hydrophobic feature. The superimposition results indicated that the four structurally similar hits match the 3CL^pro-model very well. The structural and chemical similarities of four hits suggested good pharmacophore mapping and docking scores. Furthermore, additional characterization of hits one to four such as the molecular weight (mol_MW), number of hydrogen bond acceptors (nHA), number of hydrogen bond donors (nHD), log of the octanol/water partition coefficient (logP), and log of the aqueous solubility (LogS) was evaluated and their parameter values are in the optimal range (Supplementary Table S2), suggesting their druggable and pharmacokinetics properties. Finally, four candidate hits (hits 1–4) were used as candidate molecules for performing molecular dynamics simulations.

To evaluate the binding stability of each 3CL^pro-hit complex from the docking study, we conducted 50 ns MD simulation to analyze the parameters such as root mean square deviation (RMSD) and root mean square fluctuation (RMSF). As shown in Figure 4, the movement of 3CL^pro protein and each hit-3CL^pro complex was monitored. We observed that 3CL^pro protein and each 3CL^pro-hit complex can exhibit stable internal motion throughout the simulation process. The RMSF was also calculated that provides knowledge on the flexibility of the protein residues. Figure 5 shows that compared with that of 3CL^pro protein alone, the RMSF fluctuation values of the active residues (Glu166, Phe140, Leu141, Met165, Met49, Thr25, and Thr24) from all hit-3CL^pro complexes were relatively small. In addition, we also monitor the movement of each ligand in the four complexes. As depicted in Supplementary Figure S1, the RMSD value for each individual hit (from each 3CL^pro-hit complex during MD simulation) is below 0.33 nm, indicating its stability at the binding site. Overall, these data show that four hits can stably bind to the active site of 3CL^pro.

To characterize hits one to four targeting 3CL^pro in vitro, we used the MST assay to measure the binding affinity of hits one to four to 3CL^pro. The previously reported 3CL^pro inhibitor ML300 was used as a positive control. MST assay showed that the K_d values of hits one to four to 3CL^pro were ranged from 0.012 μM to 0.27 μM (Supplementary Figure S2). We also evaluated the inhibition effect of hits one to four on 3CL^pro by enzyme inhibition assay. As shown in Table 2 and Supplementary Figure S3, all the four hits inhibited SARS-CoV-2-3CL^pro activity, with IC₅₀ values ranging from 0.017 to 0.83 μM. In particular, hit one is a promising 3CL^pro inhibitor. Its inhibition activity (IC₅₀ = 0.017 ± 0.003 µM) was about 236 times stronger than that of ML300 (IC₅₀ = 4.01 ± 0.66 µM). Through the search analysis of PubChem and SciFinder, it was found that hits one to four with the inhibition activity of 3CL^pro were reported for the first time, indicating that they are novel SARS-CoV-2-3CL^pro inhibitors. These experimental data demonstrate that our virtual screening protocol is very reliable in identifying novel and effective 3CL^pro inhibitors.

The binding modes of hit 1 with SARS-CoV-2-3CL^pro were further analyzed. The oxygen atom of hit one formed a hydrogen-bond interaction with Glu166 and hydrophobic interactions with key residues including Phe140, Leu141, Met165, Met49, Thr25, and Thr24 (Figure 6A). As shown in Figure 6B, it can be observed that the 3CL^pro-active site has multiple hydrophobic sites, and the active compound targeting SARS-CoV-2-3CL^pro should occupy each site at the same time. Our molecular docking results indicated that hit one showed excellent geometric matching with SARS-CoV-2-3CL^pro. In addition, we also measured the inter-atomic distance profile of the important interacting atom pairs (hydrogen bond: for details see Supplementary Table S3) of the 3CL^pro and hits during the MD simulation. The result indicates that the oxygen atoms of all hits can form hydrogen-bond interactions with Glu166. In order to further evaluate the stability of the 3CL^pro in complex with hits one to four, we also used DSSP algorithm to monitor the secondary structure changes of each 3CL^pro-hit complex during MD simulation (Supplementary Figure S4). The result suggests that the changes of structural elements (such as α- Helix and β-Sheet content) were not observed, indicating the stability of each hit-complex system. As shown in Supplementary Figure S5, it can be clearly observed that the 3CL^pro complexes superpose well with the 3CL^pro alone with a Cα-RMSD < 2.5 Å, indicating that all systems retained the structural integrity and maintained the same fold with minor changes of the loop region in the 3CL^pro of SARS-Cov-2.