The structures of compounds were obtained from the Coconut database, ³ and the crystal structure of 3CL^pro (6W63) was downloaded from the PDB database. ⁴ Discovery Studio software (version 4.0, Biovea Inc., Omaha, NE, USA) was applied to molecule docking. The compounds of the Coconut database were prepared using the “Prepare Ligands” module. After removing crystallographic water molecules and ligands, the “Prepare Protein” module was used to remove the conformation of the target protein and to fill the incomplete amino acid residues. Then, the virtual screening was performed in the “LibDock” module, which is a virtual high-throughput algorithm for docking ligands into the receptor active sites. CDOCKER, as an implementation of a CHARMm based docking tool using a rigid receptor, was used to further molecular docking. Ligand conformations were aligned to receptor interaction sites (hotspots) and the best scoring poses were reported (Wu et al., 2003; Mohan et al., 2015).

In order to verify the time-dependent stability of ligand-receptor complexes, Desmond_maestro 2020.4 with OPLS_2005 force filed ⁵ tools was used to perform a 100-ns molecular dynamics simulation on the receptor-ligand complexes obtained by virtual screening (Okimoto et al., 2009; Gahlawat et al., 2020). Firstly, in an orthorhombic box, the system was solvated with the SPC water model and neutralized by Cl⁻ or Na⁺ ions. The system energy was minimized at T = 0 K while the force field could be used to regulate the structure of the molecules, followed by ensuring the system stability by the equilibration processes and heating. Subsequently, MD simulation was carried out under periodic boundary conditions and NPT ensemble at T = 300 K and 1 atm pressure using a Nose–Hoover thermostat and Martyna–Tuckerman–Klein barostat (Li et al., 2019; Srivastava et al., 2021). The parameters, such as root-mean-square-deviation (RMSD) for the receptor and ligand, root-mean-square-fluctuation (RMSF) of the protein and ligand, or numbers of hydrogen bonds, were calculated to evaluate the stability of receptor-ligand complexes. The change in conformation can be analyzed by the change in RMSD, showing whether the equilibrium is reached during the simulation. When the shift between the specific frame and the reference frame is 1–3 Å, it can be considered that the shift of the specific frame is not much different from the reference frame. However, it indicates that a larger conformational change has occurred during the simulation when the shift change is apparent. The backbone RMSD curve represents the RMSD of the protein backbone over time. Lig_fit_Prot curve represents the RMSD of the ligand when the protein-ligand complex is first aligned to the protein backbone and then the RMSD of the heavy atom is measured. It is likely that the ligand has diffused away from its initial binding site if a value significantly larger than the RMSD of the protein is observed. The lig_fit_lig curve represents the RMSD of the ligand when the ligand is aligned in its reference conformation which is used to measure the internal fluctuations of the ligand atoms.

AdmetSAR2 ⁶ is a website which allows user to search for ADME/T (Absorption, Distribution, Metabolism, Excretion, and Toxicity including Ames toxicity, carcinogenic properties, acute oral toxicity, and rat acute toxicity) properties profiling by name, CASRN, and similarity search (Yang et al., 2019). At the same time, Swiss-ADME online web server ⁷ was used to predict physicochemical properties, water solubility, pharmacokinetics, lipophilicity, drug-likeness drug-like nature (according to drug-likeness rules: Lipinski, Ghose, Veber, Egan, and Muegge), and medicinal chemistry friendliness (PAINS, Brenk, Lead-likeness, Synthetic accessibility) (Daina et al., 2017).

With the rapid development of various computer methods such as molecular modeling, molecular docking, and molecular dynamics simulations, modern computer tools provide new strategies to explore the scientific basis of TCM (Jiao et al., 2021). Molecular docking is an established in silico structure-based method widely used in drug discovery. Docking enables the identification of novel compounds of therapeutic interest, predicting ligand-target interactions at the atom level, or delineating structure-activity relationships (SAR), without knowing a priori the chemical structure of other target modulators (Pinzi and Rastelli, 2019).

It has been shown that structural domain III (residues 201-301) of the 3CL^pro has catalytic activity (Tahir Ul Qamar et al., 2020). X77 can bind at the active site to inhibit the catalytic activity of the 3CL^pro (Gahlawat et al., 2020). The interaction diagram of inhibitor X77 bound to 3CL^pro is shown in Figure 1A. X77 bound to 3CL^pro by two hydrogen bonds between residues GLY143 and GLU166, Pi-Sulfur interaction with residues MET49, amide-Pi stacked with residues LEU141, and other hydrophobic interactions. Very interestingly, as shown in the middle part in Figure 1A, X77 is docked into the active pocket by an “X” shape.

The top 20 virtual hits with different structural skeletons are listed in Table 1. The binding affinities of all compounds were assessed through LibDock Scores and -CDOCKER energy value. According to the LibDock score results, all compounds had apparent interactions with the target protein. The top 20 virtual hits showed stronger interactions with 3CL^pro than the original inhibitor X77. Two -CDOCKER interaction energy values of 38.6734 and 35.6722 kcal/mol, compared to that of X77 (-CDOCKER_ENERGY = 30.4747), were calculated, indicating that CNP0348829 and CNP0474002 are promising candidate compounds for the treatment of SARS-CoV-2. The ADMET property predictions of the virtual hits were demonstrated in Tables 2, 3 in order to investigate how molecules can access the target site of 3CL^pro after entering the bloodstream. This analysis is also crucial for analyzing the efficacy of molecules (Fadaka et al., 2022). All parameters were within the Lipinski’s rule of five (ROF) cut-off range for the test compounds except CNP0465303 and CNP0447420 and present no bystander toxicity effects, since managing toxicity is the main task in developing new medications. Ames toxicity, carcinogenic properties, and rat acute toxicity which are important parameters of drug safety, were predicted in the current investigation. According to the results, most of the compounds showed no Ames toxicity, no carcinogenic properties, and low toxicity, meaning they are safe when used as drugs (Fadaka et al., 2021; 2022). The ideal drug molecules fall under such criteria as molecular weight (<500 g/mol), MlogP values (≤4.5), rotatable bonds (≤10), HB acceptors (<10), and number of HB donors (<5) (Lipinski, 2004). We screened out two star-hits based on ROF, while others were eliminated because of the molecular weight or -CDOCKER_ENERGY values.

Interaction diagrams of CNP0348829 and CNP0474002 bound to 3CL^pro are shown in Figures 1B, C. As seen in Figure 1B, CNP0348829 was stabilized by three strong hydrogen bond interactions between residues SER144, CYS145, and HIS163. In addition, the ring exhibited very strong Pi–Pi stacked interactions with HIS41, Pi–sulfur interaction with CYS44. CNP0474002-3CL^pro bound to 3CL^pro by hydrogen bond interactions between C=O of amide with CYS145. The Pi–alkyl stacked interactions between residues MET165 and 18-phenyl and X–sulfur interaction with MET49 like CNP0348829 were also observed (Figure 1C; Table 4). Both CNP0348829 and CNP0474002 bound to the 3CL^pro in a spreading conformation, CNP0474002 was docked into the active pocket by an “X” shape, somewhat like X77.

To verify the results of molecular docking, molecular dynamics simulations were introduced in this study. Molecular dynamics simulations can show the microscopic evolution of the system at the atomic level and visualize the mechanisms and patterns of experimental phenomena, which is a great advantage in studying the mechanism of ligand-receptor interaction (Adcock and McCammon, 2006).

The change in conformation can be analyzed by the value of RMSD, showing whether the equilibrium is reached during the simulation. When the position-shift between the specific frame and the reference frame is 1–3 Å, it can be considered that the position-shift of the specific frame is not much different from the reference frame. It indicates that a larger conformational change occurred during the simulation when the position-shift change is apparent. The backbone RMSD curve represents the RMSD of the protein backbone over time. Lig_fit_Prot curve represents the RMSD of the ligand when the protein-ligand complex is first aligned to the protein backbone and then the RMSD of the heavy atom is measured. It is likely that the ligand has shifted away from its initial binding site if the value is significantly larger than the RMSD of the protein. The lig_fit_lig curve represents the RMSD of the ligand when the ligand is aligned in its reference conformation which is used to measure the internal fluctuations of the ligand atoms.

The results of molecular dynamics simulations are shown in Figures 2, 3. The major deviation is in RMSD-P, which shows the protein RMSD with a ligand in the active site. The compound CNP0348829 showed initial light fluctuation in RMSD-L, at around 1 Å and 3 Å. However, it showed very strong fluctuations at 52 ns, mainly resulting from the change of conformation at the 12-phenyl, while remaining in the active site of the receptor and stabilized for the remaining period of time in the range of 3.5 Å to 5.5 Å. (Figures 2A, C, D). Unlike CNP0348829, CNP0474002 was stabilized in the range of 3 Å to 5.5 Å at all times under study (Figure 3A). The major initial fluctuation in the 3CL^pro was mainly due to C-terminal fluctuations, which was confirmed by the RMSF plot for the protein C (alpha) and backbone, where except for C-terminal, all other residues showed RMSF less than 2.5 Å. Over the period of analysis, C-terminal exhibited RMSF up to 4.5 Å (Figures 2B, 3B). RMSF of the ligand with the protein and the ligand with initial conformation was very stable and was between 1.00 and 2.25 Å and 0.25 to 1.25 Å, respectively (Figure 3C). Protein interactions with the ligands were monitored throughout the simulation, and these interactions can be summarized and classified into four types: hydrogen bonds, hydrophobic interactions, ionic bonds, and water bridges. Among all interactions, hydrogen bonds play an important role in the overall dynamics simulation because of their greater influence on compound specificity, metabolic characteristics, and adsorption parameters (Kruth et al., 2005). The following amino acids, such as HIS41, CYS145, and GLU166 in CNP0474002-3CL^pro complex, and CYS44, ASN142, GLY143, SER144, and GLU166 in CNP0348829-3CL^pro complex, contribute greatly to the interaction between protein and ligand. These interactions were categorized by type and summarized as shown in Figures 4A, B. The above analysis suggests that the ligand-receptor complexes were maintained well within a relatively stable range with small fluctuations throughout the molecular dynamic simulations, indicating that the binding modes of CNP0348829-3CL^pro and CNP0474002-3CL^pro complexes are stable during the molecular dynamics simulation.