Structure-Based Discovery of Novel Nonpeptide Inhibitors Targeting SARS-CoV-2 Mpro.

The continual spread of novel coronavirus disease 2019 (COVID-19) is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), posing a severe threat to the health worldwide. The main protease (Mpro, alias 3CLpro) of SARS-CoV-2 is a crucial enzyme for the maturation of viral particles and is a very attractive target for designing drugs to treat COVID-19. Here, we propose a multiple conformation-based virtual screening strategy to discover inhibitors that can target SARS-CoV-2 Mpro. Based on this strategy, nine Mpro structures and a protein mimetics library with 8960 commercially available compounds were prepared to carry out ensemble docking for the first time. Five of the nine structures are apo forms presented in different conformations, whereas the other four structures are holo forms complexed with different ligands. The surface plasmon resonance assay revealed that 6 out of 49 compounds had the ability to bind to SARS-CoV-2 Mpro. The fluorescence resonance energy transfer experiment showed that the biochemical half-maximal inhibitory concentration (IC50) values of the six compounds could hamper Mpro activities ranged from 0.69 ± 0.05 to 2.05 ± 0.92 μM. Evaluation of antiviral activity using the cell-based assay indicated that two compounds (Z1244904919 and Z1759961356) could strongly inhibit the cytopathic effect and reduce replication of the living virus in Vero E6 cells with the half-maximal effective concentrations (EC50) of 4.98 ± 1.83 and 8.52 ± 0.92 μM, respectively. The mechanism of the action for the two inhibitors were further elucidated at the molecular level by molecular dynamics simulation and subsequent binding free energy analysis. As a result, the discovered noncovalent reversible inhibitors with novel scaffolds are promising antiviral drug candidates, which may be used to develop the treatment of COVID-19.

Introduction

Infection with severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) will cause novel coronavirus disease 2019 (COVID-19),[1] and the pandemic of the disease has rapidly become a global health concern[2] and led to 160,074,167 confirmed cases and 3,325,260 deaths worldwide as of May 13, 2021.[1] To cope with the severe crisis, great efforts have been paid to developing therapeutic approaches and vaccines against SARS-CoV-2.[3,4] Discovering inhibitors of key proteins involved in the viral life cycle is an often-used and efficient approach to disrupt the replication of virus.[5] Like SARS-CoV, the encoded 4 structural and 16 nonstructural proteins (NSPs) of SARS-CoV-2 provide multiple avenues to identify potential drug targets.[6,7] Among the encoded proteins, the main protease (Mpro, alias 3CLpro), which has no human homolog, has become an attractive therapeutical target for the drug discovery and development of anti-COVID-19.[8,9]

Mpro belongs to the 16 NSPs of coronavirus (CoV) and is a vital enzyme that has an essential role in mediating the replication and transcription of CoVs.[8] Together with papain-like proteases (PLPs), the enzyme processes the polyproteins that are translated from CoV RNA.[10] Mpro is a highly conservative protein existing in all CoVs consisting of three domains (domains I to III).[8] Crystal structures of SARS-CoV-2 Mpro (Figure )[9,11] show that they are the chymotrypsin-like domain (domain I, residues 10 to 99), picornavirus 3C protease-like domain (domain II, residues 100 to 182), and a globular cluster formed by five helices (domain III, residues 198 to 303). The substrate-binding site (active site) of Mpro composed of four subsites (S1, S2, S3, and S4) is located at the six-stranded antiparallel β barrels between domains I and II.[9]

Figure 1

(A) Workflow of ensemble docking-based virtual screening of novel nonpeptide inhibitors targeting SARS-CoV-2 Mpro. (B) Ensemble SARS-CoV-2 Mpro 3D structures shown in cartoon representation with different colors. Domain I (residues 10 to 99), Domain II (residues 100 to 182), and Domain III (residues 198 to 303) of the protease are labeled. The substrate-binding site (active site) of Mpro composed of four subsites (S1′, S1, S2, and S4) marked by the gray surface.

Based on the crystal structures of SARS-CoV or SARS-CoV-2 Mpro, computer-aided drug design techniques have been successfully used in anti-COVID-19 studies regarding the rapid discovery of potential inhibitors,[12−16] drug repurposing,[14,16−20] and making the action mechanism of the active compound against SARS-CoV-2 more understandable.[21] Though these timely research studies have led to the design of several first-in-class SARS-CoV-2 Mpro inhibitors as promising drug candidates,[8,9,11] currently no Mpro-based therapeutics have been officially approved for COVID-19.[3] The need to develop novel as well as more effective antiviral drugs to inhibit SARS-CoV-2 has become more urgent.[3] However, larger flexibility and figurability of active sites on SARS-CoV-2 Mpro proved to be a challenge for the rational design of small molecule inhibitors.[22,23] For addressing this problem, the crystal structures of Mpro could be complemented by the all-atom molecular dynamics (MD) trajectory data released publicly in the spirit of open science.[24,25]

In the present work, based on nine different conformations of the SARS-CoV-2 Mpro substrate-binding site, a multiple conformational-based virtual screening strategy in combination with experimental validation was proposed to identify the enzyme inhibitors from a protein mimetics library with 8960 commercially available compounds (Figure ). Considering the docking pose and scaffold diversity, 49 selected candidates were purchased for testing their binding profiles to SARS-CoV-2 Mpro using the surface plasmon resonance (SPR) assay. The identified six compounds were further evaluated by the fluorescence resonance energy transfer assay (enzyme kinetics study) and bioluminescence resonance energy transfer (BRET) assay. All six compounds showed inhibition activities against the cell lines of SARS-CoV-2 Mpro. The live virus assay indicated that two out of the six inhibitors had the activity to interdict the viral infection of SARS-CoV-2. In addition, computational absorption-distribution-metabolism-excretion (ADME) analysis showed that the two inhibitors had good pharmacokinetic properties and low toxicity.

Results and Discussion

Compounds Selected through Ensemble Docking

To account for protein flexibility, nine Mpro ensembles including five MD-sampled apo structures and four holo structures (one homology model and three crystal structures in complex with different ligands) were collected. Meanwhile, a protein mimetics library (with 8960 compounds) from Enamine was prepared and resulted in 28,727 conformations. As an initial step, ensemble docking of the prepared library against the active site of the defined nine Mpro ensembles was performed to identify potential lead candidates. By considering the GlideScore, binding pose, and scaffold diversity profiles, the nine sets of hits from the ensemble docking were then used for selecting 50 top-ranked compounds (with 49 purchased) for experimental testing using SPR assays. As shown in Table S1, the selected compounds as potential Mpro inhibitors have GlideScore ≤−5.651 kcal/mol, and each of the compound forms at least two hydrogen bonds with the residues located at the protease active site. The hierarchical clustering of the fingerprints using Tanimoto similarity and Ward’s cluster linkage method[26] of the selected compounds shown in Figure S1 indicated the high diversity of the scaffolds. HPLC chromatograms and mass spectrograms were applied to verify the chemical structures and purity of the 49 compounds, and the data of the six active compounds (see the next section) are provided in the Supporting Information.

Evaluation of Compounds as Inhibitors of SARS-CoV-2 Mpro In Vitro

By using dipyridamole (DIP) as a positive control, the binding toward SARS-CoV-2 Mpro of the purchased 49 compounds was tested using the SPR assay at 100 μM concentrations (Figure S2). In addition, 6 out of the 49 compounds that have the abilities of binding to Mpro (Table ) were selected to investigate whether their binding alters the biochemical function of the enzyme. The ranking of GlideScore for the six compounds Z236230776, Z1244904919, Z225729516, Z1759961356, Z108564100, and Z106460362 was 42/50, 33/50, 19/50, 25/50, 11/50, and 1/50, respectively. There is only one compound (Z106460362) that was at the top 10 of the list (Table S1). Furthermore, we determined the biochemical half-maximal inhibitory concentration (IC50) values of the six chemical compounds, ranging from 0.69 to 2.05 μM (Figure ). All compounds presented a strong inhibitory effect on Mpro activity, among which Z1759961356 (IC50 = 0.69 ± 0.05 μM) had the strongest effect (Figure ). Consistent with the IC50 results, the BRET ratio showed that all compounds had a good inhibitory effect on Mpro in HEK293T cells (Figure ). However, in this structure-based virtual screening study, although multiple conformation strategy was employed, the success rate was still very low (only 12% cases were correctly predicted by Glide docking). It is hypothesized that this is because those compounds were selected from specific conformations of the SARS-CoV-2 Mpro. However, according to the experimental test, the specific structure may not occupy the preferred conformation of the protease. Therefore, to increase the success rate of virtual screening, enhanced conformational sampling of the protease by state-of-the-art MD simulation is needed. In addition, the flexibility of the protease active site was not considered during each docking process, which was crucial for the protein–ligand recognition. Therefore, the induced fit docking method may be used to address this problem even if the calculation is time consuming.

Table 1

Information of Six Compounds That Have the Abilities of Binding to Mpro Using Dipyridamole (DIP) as a Positive Control

The resonance units (RU) of the SPR assay in the presence of each compound at a concentration of 100 μM.

The nine SARS-CoV-2 Mpro structures including five apo forms (extracted per 2 μs from 10 μs MD simulation of 6LU7(11)) and the four holo forms (one homology model using 3ATW(33) as a template and three crystal structures 6LU7,[11]6Y2F,[9] and 6Y2G(9) in complex with different ligands.

The docking scores (kcal/mol) were calculated by the Glide extra precision algorithm.[39]

Figure 2

Inhibitory activity profiles of compounds against SARS-CoV-2 Mpro. The median inhibitory concentration (IC50) values were determined by a fluorescence resonance energy transfer (FRET)-based cleavage assay.

Figure 3

Dose dependence of six inhibitors on intracellular SARS-CoV-2 Mpro activity measured using a bioluminescence resonance energy transfer (BRET) ratio.

Inhibitors Suppress SARS-CoV-2 Infection In Vitro

For examining whether these two lead candidates could prevent viral replication, further qRT-PCR and plaque-reduction assays were carried out in Vero E6 cells infected by SARS-CoV-2. As can be seen from Figure , quantitative qRT-PCR results showed that Z1244904919 and Z1759961356 exhibited a stronger effect on anti-SARS-CoV-2 (Figure A,B). The plaque-reduction assay indicated that Z1244904919 and Z1759961356 displayed inhibitory effect on SARS-CoV-2, and the individual EC50 values were 4.98 ± and 8.52 ± μM, respectively (Figure C,D). Furthermore, the SPR assay showed that Z1244904919 and Z1759961356 bound to SARS-CoV-2 Mpro with Kd values of 465 and 133 μM, respectively (Figure A,B). In conclusion, these data suggest that the inhibition of Z1244904919 and Z1759961356 on Mpro is mainly achieved through direct binding to the enzyme active site.

Figure 4

In vitro inhibition of viral main protease inhibitors against SARS-CoV-2. (A) At 72 h after infection, viral RNA (vRNA) copy numbers in Vero E6 cells monitored by qRT-PCR. (B–D) Mean percent inhibition of virus yield in the cells treated with a series concentration of DIP, Z1244904919, and Z1759961356.

Figure 5

Binding of inhibitors to SARS-CoV-2 Mpro. (A, B) Surface plasmon resonance (SPR) assay of Z1244904919 and Z1759961356 bound to the protease. (C, D) The binding modes and (E, F) energy contributions of key residues in the Mpro active site for Z1244904919 and Z1759961356 to the protease. The protein and ligand were displayed as cartoon and stick representation, respectively. The hydrogen bond is shown in green dashed lines.

MD Simulation of the Inhibitor–Mpro Complex

Though the two lead candidates were recognized by ensemble docking, we thought that their predicted binding modes in Mpro were not enough because the protease flexibility was not considered in each independent docking. To investigate inhibitor–Mpro interaction flexibility, 1 μs MD simulation was executed for sampling enough conformations of the two complexes. The root-mean-square deviation (RMSD) of the backbone atoms on protein and heavy atoms on the ligand referred to the starting structure was computed to reflect the stabilities of the studied systems during the period of simulation (Figure S3). The RMSD value variation suggested that the two complexes had small changes of conformation on the process of simulation. The average RMSD values of the binding site residues for Z1244904919 and Z1759961356 bound Mpro were 0.88 and 1.86 Å, respectively. The values for Z1244904919 and Z1759961356 were 1.09 and 1.96 Å. The trends of RMSD variation in Figure S3 indicated that the poses of ligands predicted were consistent with the active site of Mpro. In addition, we have compared the predicted poses of inhibitors Z1244904919 and Z1759961356 with positive control DIP in Mpro (Figure S4A,B). The results showed that the three ligands occupied the same binding site of the protease.

Binding Free Energy and Interaction Mode of Inhibitors in Mpro

On the account of the MD trajectories, the binding free energy of the two inhibitors bound to Mpro (ΔGcalc) was estimated using the MM/GBSA method.[27] As shown in Table , the ΔGcalc for Z1244904919 and Z1759961356 bound to Mpro was −45.72 and −48.01 kcal/mol, respectively. The variation trend of ΔGcalc values is compatible with the order of the experimental binding free energies (ΔGexp). The energy terms of ΔGcalc are listed in Table , indicating that the electrostatic (ΔEele) and hydrophobic (ΔEvdW + ΔGnonpol) interactions were of great importance for the binding of the four anticoagulants; however, polar solvent energies (ΔGpolar) were not conducive to the binding of inhibitors. In order to acquire a more particular understanding of the protein–ligand interaction, we decomposed the binding free energies into each residue. Residues with an absolute energy contribution of ≥0.5 kcal/mol would be identified as key residues, which were conducive to the binding of inhibitors to the pocket; these key residues are displayed in Table S2. Meanwhile, the recognized key residues of the two complexes suggested that there was a certain degree of similar interactions between them. As shown in Table S2, a total of 14 and 13 residues in SARS-CoV-2 Mpro were identified to play an important role in Z1244904919 and Z1759961356 binding, respectively. Compared with the characterized interactions between the protease with the substrate[28] and N3,[29] 10, 6, 7, and 5 common residues were found for Z1244904919- and Z1759961356-bound complexes (Figure S5), indicating that the key interactions between the protease pocket and ligands were maintained for the identified new nonpeptide inhibitors. Meanwhile, the superposition between SARS-CoV-2 Mpro in complex with N3 and Z1244904919 (Figure S4C) and ZN1759961356 (Figure D) indicates the overlap between the occupied pockets of these inhibitors, especially N3 and ZN1759961356.

Table 2

Biochemical Half-Maximal Inhibitory Concentration (IC50, μM) and Binding Free Energies (ΔG, kcal/mol) of Z1244904919 and Z1759961356 to Mpro

complexes	ΔEelea	ΔEvdWb	ΔGpolc	ΔGnonpold	ΔGcalce	IC50	ΔGexpf
Z1244904919-Mpro	–18.76	–46.88	–65.63	–3.81	–45.72	0.73 ± 0.05	–8.70
Z1759961356-Mpro	–6.03	–52.51	–58.54	–4.25	–48.01	0.69 ± 0.05	–8.73

Electrostatic (ΔEele) energy terms in the gas phase.

van der Waals (ΔEvdW) energy terms in the gas phase.

Polar (ΔEpol) solvent energies by solving the GB equation.

Nonpolar (ΔEnonpol) solvent energies by solving the GB equation.

Calculated binding free energy (ΔGcalc), ΔGcalc = ΔEele + ΔEvdW + ΔEpol + ΔEnonpol.

Experimental binding free energy (ΔGexp), ΔGexp ≈ RTln(IC50).

The binding modes of Z1244904919 (Figure C) and Z1759961356 (Figure D) were investigated by the representative conformations extracted from the MD trajectories. As is known, the active site of Mpro consists of four subpockets, which are S1, S2, S3, and S4.[9] Residues Leu27, His41, Met49, His164, Met165, and Gln189 identified as key residues were of great importance for both Z1244904919-Mpro (Figure E) and Z1759961356-Mpro (Figure F) complexes. All the identified key residues uniformly distributed in the four subpockets of the Mpro active site (Figure C,D). Taking Z1244904919 as an example (Figure C), the backbone atoms of Gly143, Ser144, Cys145, and Asn166 interact with the compound via hydrogen bonds. The fluorophenol moiety of Z1244904919 embedded into the S1 site consisted of residues Phe140 and Asn166, and the piperidine moiety took up the S4 site containing residues M165 and Gln189, while the indole analogue moiety and linkages in contact with residues Met49, Thr25, and His41 located at the S2 and S3 sites. The piperidine moiety acts like a linker to connect fluorophenol and indole analogue motifs. Compared to Z1244904919, the higher binding abilities of the Z1759961356 may come from the energy contributions of residues His164 and Met165 in the S1 pocket and residue Asn47 in the S2 pocket of the Mpro active site. In this study, histidine (His41, His163, and His164) and cysteine (Cys145) located at the binding site of inhibitors were treated as neutral states during docking and MD simulation. However, it is important to note that the altering protonation states of titratable groups in histidine and cysteine in SARS-CoV-2 Mpro, which can modulate protein dynamics and stability, is important in virtual screening studies. This has been well studied in the recently published work by Pavlova et al.[30]

In Silico Pharmacokinetic Analysis

The pharmacokinetic properties of new lead candidates are essential for the development of an effective druggable molecule. Herein, the ADME properties of Z1244904919 and Z1759961356 were calculated in QikProp (v. 4.5) (Table ). The QikProp method is based on 1700 known oral drugs, and the rms errors of its predictions are 0.5–0.6 log unit.[31] The ADME properties of the recently reported two SARS-CoV-2 Mpro inhibitors (13a and 13b)[9] and FDA approved drug DIP,[32] which have favorable pharmacokinetic properties, were also calculated and are included in Table . The predicted ADME values of Z1244904919 and Z1759961356 compare favorably with the drug leads 13a and 13b or FDA approved drug DIP (Table ). Moreover, some properties, such as QPPCacod and PercentHumanOralAbsorption, are better than those of 13a, 13b, and DIP. Therefore, we are optimistic about the application of the two compounds Z1244904919 and Z1759961356 as new drug leads targeting the Mpro protein.

Table 3

Calculated Pharmacokinetic Properties of Compounds from QikProp (Version 4.5)

compounds	MWa	QPlogPo/wb	QPlogSc	QPPCacod	absorptione
Z1244904919	383.465	3.408	–3.241	255.308	89.98
Z1759961356	382.505	2.899	–2.93	550.893	92.983
13af	585.699	2.782	–5.863	29.197	43.542
13bf	593.678	2.682	–3.504	67.514	49.474
DIPg	504.631	1.974	–3.468	120.031	49.804

The molecular weight of the molecule.

The predicted log of the octanol/water partition coefficient.

The predicted aqueous solubility; S in mol/L.

Predicted Caco-2 cell permeability in nm/s.

Predicted percent of human oral absorption (%).

13a and 13b are two recently reported two SARS-CoV-2 Mpro inhibitors with favorable pharmacokinetic properties.[9]

DIP is an FDA approved drug[32] S4. The range or recommand values of MW (130.0 to 725.0), QPlogPo/w (−2.0 to 6.5), QPlogS (−6.5 to 0.5), QPPCaco (>500 is great; <25 is poor), QPlogBB (−3.0 to 1.2), and PercentHumanOralAbsorption (>80% is high; <25% is poor).

Conclusions

We report that the IC50 values of the six identified inhibitors targeting Mpro ranged from 0.68 to 2.05 μM here. Among them, Z1244904919 and Z1759961356 inhibit the purified recombinant SARS-CoV-2 Mpro and their IC50 values were 0.73 ± 0.04 and 0.69 ± 0.05 μM, individually. Further experiments show that Z1244904919 and Z1759961356 display inhibition against SARS-CoV-2, and EC50 values were 4.98 ± and 8.52 ± μM, respectively. In addition to this, the recognized key residues that contributed to the binding modes between Z1244904919 and Z1759961356 with SARS-CoV-2 Mpro were elucidated by MD simulation and binding free energy analysis. The results from this study provide a new starting point for the design of antiviral drugs to treat COVID-19.

Materials and Methods

Multiple Conformation-Based Virtual Screening

Protein Preparation and Grid Generation

In this study, nine SARS-CoV-2 Mpro structures including five apo forms and four holo forms were collected from released public data. The five apo forms (no ligand bound) were extracted from 10 μs MD simulation of SARS-CoV-2 Mpro (PDB ID: 6LU7(11)) per 2 μs,[25] and the four holo forms were homolog model (the crystal structure of SARS-CoV Mpro3ATW[33] was regarded as a template) or crystal structures (PDB IDs: 6LU7,[11]6Y2F,[9] and 6Y2G(9)) of SARS-CoV-2 3CLpro with diverse ligands. The Protein Preparation Wizard[34] was used to add hydrogen atoms, assign partial charges, assign protonation states, and minimize the structure with the OPLS3 force field[35] to prepare each structure. When the RMSD value reached the maximum of 0.30 Å, this minimization would be terminated. After minimization, the Receptor Grid Generation program of Glide (Version 6.8) was used to define the docking grid for each Mpro in the monomer state.[36] For each structure, by centering on the ligand (holo form) or selecting active site residues (apo form), the docking grids were generated. By using the cocrystal structure 6LU7 as a reference, the active site residues of SARS-CoV-2 3CLpro in the apo form were chosen. The center and size of defined nine docking grids are summarized in Table S3.

Small Molecule Database Preparation

A protein mimetics library with 8960 commercially available compounds from Enamine was used for ensemble docking.[37] The library was prepared by using the LigPrep (Version 3.5) program, and then, all compounds were processed through generating tautomers, stereoisomers, and ionization states by Epik (Version 3.2).[38] All the ligand preparations were under the condition of 7.0 ± 2.0 pH value with the OPLS3 force field.[35] QikProp (Version 4.5) was used for calculating the five compounds’ ADME properties summarized in Table , and the library was prefiltered using druglike properties.

Ensemble Docking

Screening the prepared library via docking them into the generated grids using Glide (Version 6.8).[36,39] High-throughput virtual screening was first carried out for maintaining 10% top-ranked structures, and those molecules were redocked on the scoring algorithm of standard precision, reserving the 10% top-scored molecules. The resulting set was further filtered at the extra precision level, and a database involving 245 compounds was retained ultimately. From the retained sub-database, 50 compounds were selected by considering the docking scores, binding mode, and scaffold diversity. Finally, 49 compounds available commercially were purchased from TargetMol for further biological evaluation.

Biacore Assay

Performing SPR experiments in a Biacore 8K device (Cytiva, Previously GE Healthcare Life Sciences) using CM5 sensor chips (Cytiva, Previously GE Healthcare Life Sciences) on the basis of the protocol provided by the manufacturer. Briefly, recombinant SARS-CoV-2 Mpro protein was fixed in a CM5 chip. Compounds of different concentrations were injected at a flow rate of 30 μL/min lasting for 2 min. Subsequently, collecting data for a 2 min association followed by a 5 min dissociation. The chip was regenerated by injecting 1 × PBS, 0.05% Tween-20, pH 7.4, 5% DMSO for 60 s. All procedures were run in 1 × PBS, 0.05% Tween-20, pH 7.4, 5% DMSO as a running buffer. The software Biacore Insight Evaluation Software with a 1:1 Langmuir binding model was applied to analyze the binding kinetic. The Kd was calculated by the Biacore Insight Evaluation Software.

Mpro Activity and Inhibition Assay

Fluorescence Resonance Energy Transfer

Chemical compounds were dissolved in 100% DMSO. Half-maximal inhibitory concentration (IC50) was determined using the 3CL Protease, MBP-tagged (SARS-CoV-2) Assay Kit (BPS Bioscience, San Diego, CA, USA). In brief, 3CL protease (5 ng/μL) was preincubated with chemical compounds at indoor temperature for 30 min with slow shaking. Afterward, a substrate solution with a 50 μM final concentration was added to each well to initiate the reaction. The samples were incubated overnight at indoor temperature. The fluorescence intensity was surveyed at 360 nm excitation. GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA) was used for the calculation of the IC50 values.

Bioluminescence Resonance Energy Transfer

BRET was used to detect the inhibitory effect of the chemical compounds against SARS-CoV-2 Mpro in HEK293T cells. The 3CL protease recognition sequence linker (ITSAVLQSGFRK) was fused with enhanced yellow fluorescent protein (EYFP) and inserted into pRLuc-N2 plasmid. A full-length coding sequence of 3CL protease was inserted into pcDNA3.1-Flag plasmid. pEYFP-linker-Rluc and pcDNA3.1-3CL-Flag were co-transfected into HEK293T cells and treated with chemical compounds in different concentrations for 48 h. The BRET ratio was detected at 475 nm emission.

Antiviral Activity Assay

qRT-PCR Assay

The in vitro antiviral efficacy of compounds was determined in Vero E6 cells as previously described.[11] Briefly, the cells were pretreated with the chemical compound with a concentration of 10 μM for 1 h and then infected with SARS-CoV-2 with multiplicity of infection of 0.01 for 2 h. After this, the virus–drug mixture was wiped out, and the cells were placed in the medium filled with fresh drugs for further cultivation. At 72 h post infection, viral RNA (vRNA) was extracted from the culture supernatant and detected by quantitative real-time PCR (qRT-PCR).

Plaque-Reduction Assay

The anti-SARS-CoV-2 activity of selective compounds was determined through the plaque-reduction assay. Compounds in different dilution concentrations were mixed with SARS-CoV-2 (100 plaque-forming units), and 200 μL of mixtures was injected into 1 × 105 monolayer Vero E6 cells lasted for 1 h. Then, the cells were washed twice with a fresh medium; after this, the cells were incubated with 0.9% agarose containing indicated chemical compounds. After infection, at the 4th day, the cells were fixed in 4% polyoxymethylene for 30 min and finally dyed with crystal violet. The plaque-forming units were counted.

MD Simulation and Binding Free Energy Calculation

The docking poses of Z1244904919- and Z1759961356-bound Mpro complexes were used to perform MD simulation by the GPU-accelerated PMEMD module of AMBER14 software as previously described.[40] Before MD simulation, the AMBER ff14SB[41] was assigned to Mpro protein and Antechamber[42] with GAFF[43] and RESP partial charges were applied for two inhibitors to generate the force field parameters. The HF/6-31G* level of Gaussian09 suite[44] was employed for the calculations of ligand geometric optimization and the electrostatic potential. Then, the two complexes were neutralized through adding an appropriate number of counterions and immersed into a rectangular periodic box of TIP3P[45] water molecules with an edge of 10.0 Å. For each complex, two steps of 1000 cycles of energy minimization were performed, the first one was under a harmonic restraint of a 10.0 kcal·mol–1·Å–2 force constant followed by the second minimization without restraint. After this, the two complexes were heated from 0 to 100.0 K in 2500 steps and gradually to 310.0 K within 5000 steps; both of them are under a force constant of 10.0 kcal·mol–1·Å–2. Then, equilibration with 50 ps at 310.0 K was conducted by freeing all atoms. At the end, production run with 1000 ns was performed for the two systems under the NPT ensemble at 310.0 K and 1 atm by the periodic boundary condition.

The binding free energies (ΔGcalc) of Z1244904919-Mpro and Z1759961356-Mpro complexes were calculated via the end-point molecular mechanics generalized Born surface area (MM/GBSA) approach[27] using the following equation (eq ):

Furthermore, we decomposed the total binding free energy into each residue by eq to recognize the key residues responsible for the binding of ligand–Mpro complexes.