In this study, we focused on RBD of Omicron VOC. The crystal structure was retrieved from the PDB library (PDB ID: 7WBP) (Patel et al., 2021b). All inorganic ions and water molecules were removed. The combined library (containing 92,699 compounds in total) was obtained from Approved Drugs in Major Juridications (https://zinc.docking.org/substances/subsets/world/), Enamine Coronavirus Library (https://enamine.net/), Asinex small molecule PPI inhibitors (http://www.asinex.com/), traditional Chinese medicine natural products (http://tcm.cmu.edu.tw/) and Alinda natural products (https://www.alinda.ru) databases. The corresponding structures of these compounds were downloaded from the ZINC database (http://zinc.docking.org/). The Openbabel software was used to separate the sdf files of individual small molecules.

The compounds were processed to protonation states and pdbqt file using the MGLTools software for docking. Docking was conducted via AutoDock Vina 1.2.3 (Trott et al., 2010; Eberhardt et al., 2021), with different parameters. The resulting docked structures were ranked according to their predicted binding energies. For the initial screening, “exhaustiveness” was set to 50. The compounds representing cluster centroids were used with a docking box of 26.25 × 45.50 × 22.50 Å and center at [−35.72, 32.42, 0.22]. The top 1,000 hits from all docking procedures were selected for further screening. For the second screening, “exhaustiveness” was set to 400. Finally, the top 10 of the 1,000 highest binding energies and best-docked conformation were considered for the next MD simulation.

In order to validate the docking procedure, we have run an enrichment test. We selected three small molecules as control molecules, lifitegrast sodium, evans blue and lumacaftor, that have been shown experimentally to have detectable binding to RBD proteins and block the recognition of RBD and ACE2 (Day et al., 2021). The wild-type RBD protein was obtained from the PDB (PDB ID: 6M0J). Ten thousand small molecules were randomly selected, and the same docking parameters as virtual screening were used for docking test.

The MD Simulations were performed using the Gromacs 2019.6 package (Van Der Spoel et al., 2005). Parameterization of the protein was conducted using the AMBER14SB force field (Maier et al., 2015). Then, the ligand was parametrized with the general AMBER force field (Wang et al., 2004) obtained from the AmberTools21 program (Case et al., 2021). The binding conformation with the highest affinity in docking was used as the simulated initial complex conformation. The resultant system was solvated in a dodecahedron box with TIP3P water (Jorgensen et al., 1983) molecules extending 10 Å from any atoms of the protein in any directions. Three sodium ions were added to maintain the neutrality of the system. The final system had 28,608 atoms. All bond lengths involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert et al., 1977). The time step for integration of equations of motion was 2 fs. The particle mesh Ewald method (Darden et al., 1993) with a cutoff distance of 10 Å was used to calculate Coulomb interactions. The steepest descent method (Hratchian et al., 2010) was performed to minimize the system with a tolerance value of 1,000.0 kJ/mol/nm. Then, all systems were sequentially minimized and equilibrated in NVT and NPT ensembles. Then, the system was heated to 300 K using the v-rescale temperature coupling scheme (Bussi et al., 2007) with the NVT ensemble in 1,000 ps, followed by another 1000-ps NPT simulation via the Parrinello Rahman pressure coupling scheme (Martyna et al., 1994). Each of the ten systems was performed for MD simulation of 100 ns.

We further assessed the binding free energies between RBD and ligands with the MD results. The MMPBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) method (Homeyer and Gohlke, 2012), implemented in the gmx_MMPBSA program (Valdés-Tresanco et al., 2021) was adopted along with AmberTools21 (Case et al., 2021). The PB model estimates only the polar component of the solvation. The non-polar component is usually assumed to be proportional to the molecule’s total solvent accessible surface area, with proportionality constant derived from experimental solvation energies of small non-polar molecules. This method has been proven to balance accuracy and computational efficiencies, especially when dealing with large systems. In the 100-ns simulation, the most stable 40-ns trajectory was selected for free energy calculation. Based on the calculated results, the 5 compounds with strongest binding abilities to RBD were used for energy decomposition, MD trajectory analysis, quantum chemical calculations and ADMET prediction.

To establish the more precise molecular properties of the top 5 compounds, the DFT method was used to perform quantum chemical calculations with the Gaussian 16 program (Frisch et al., 2016). Geometries were optimized without any constraint with the B3LYP method (Becke, 1993; Stephens et al., 1994) using 6-311G (d,p) basis (Ditchfield et al., 1970; Hehre et al., 1972; Hariharan and Pople, 1973) with DFT-D3 empirical dispersion corrections (BJ damping) (Grimme et al., 2011). The polarizable continuum model implicit solvent model (Di Remigio et al., 2016) was used to study the effects of water solvents. The Multiwfn 3.8 program (Lu and Chen, 2012a) was used for Electrostatic potential (ESP) maps generation with the five selected molecules. Then, the electrostatic potential involved in analyses was evaluated by Multiwfn based on the highly effective proposed algorithm (Lu and Chen, 2012b; Zhang and Lu, 2021). The Visual Molecular Dynamics program (Humphrey et al., 1996) was used for visualization of ESP surface. Based on Density Functional Reactivity Theory (LIU, 2016), the lowest unoccupied molecular orbital energy (ELUMO), the highest occupied molecular orbital energy (EHOMO), and other molecular chemical descriptors such as chemical hardness (η), chemical softness (S), band gap (GAPE), electron affinity (EA), ionization potential (IP), electrophilicity index (ω), electronic potential (μ), and electronegativity (χ) were calculated.

Pharmacokinetics and toxicity are important considerations in drug development. Evaluating the absorption, distribution, metabolism, excretion and toxicity (ADMET) properties (Ferreira and Andricopulo, 2019) of compounds can help in determining whether they are a potential drug candidates. We predicted the ADMET properties of the top 5 compounds obtained from the above screening tests. The ADMET prediction was performed using the ADMETlab 2.0 online server (Xiong et al., 2021) (https://admetmesh.scbdd.com/). Human intestinal absorption (HIA), through the blood−brain barrier permeability, cytochrome P450 enzyme inhibition, hepatotoxicity, as well as mitochondrial toxicity and other important indicators were selected as predictors of ADMET.

In this study, the amino acid residues where RBD interacts with ACE2 at the PPI interface are referred to as recognition residues. The mutation sites unique to Omicron VOC are described as hotspot residues. In Supplementary Figure S1, the binding energy values of 244 compounds ranged from −8.5 to −10.7 kcal/mol in the second screening. Based on their structural diversity, 50 of these compounds were selected. Then, the binding conformation of these compounds were visually inspected (Supplementary Table S1). Candidates that had good interactions with SARS CoV-2 Omicron RBD protein were selected. After screening twice with different precisions, 10 compounds were selected for subsequent analyses (Figure 1). The docking results are shown in Table 1. Each molecule has varying degrees of contact with hotspot residues (1–4). The number of recognition residues that interacted with the 10 candidate molecules ranged from 1 to 9. The greater the number of contact recognition residues in the candidate molecule, the stronger the PPI inhibition. The top 5 candidate molecules with the highest binding affinities interacted with an average of 8 recognition residues. The top 5 molecules with the highest binding affinities showed high inhibitory activities. This strongly inhibits the recognition of ACE2 receptors by RBD. The results of the enrichment test showed that the binding energy ranking of the three control molecules was all within 1%. This result confirms the reliability of the docking procedure in this study, indicating that it can be used for candidate compound enrichment. Furthermore, we docked three control molecules with Omicron RBD protein, and the results showed that their binding affinity to proteins was weaker than that of the first 10 compounds screened. This also implies that the top 10 compounds are potential RBD protein binders.

The binding mechanism of the top 5 compounds was also analyzed. Interactions of candidate molecules with the RBD protein are presented in Figure 2. ZINC95919448, which had the highest binding affinity (−10.7 kcal/mol), was tightly bound to RBD by interacting with 11 residues. It interacted with ARG403 via H-bonds. It also formed hydrophobic interactions with VAL445, TYR449, TYR453, TYR495, and TYR501. ZINC85531210 interacted with ARG403 via hydrogen bonds and also interacted with VAL453, TYR455, TYR456, TYR495 and TYR501 via hydrophobic interactions. Upon close inspection, residues TYR453, TYR495 and TYR501 were found to be important amino acid residues that form hydrophobic interactions with small molecules, which also exists in ZINC95610651, ZINC000035399302, and ZINC95910594. Additionally, residues SER494 and HIS505 interacted with the ZINC95610651 hydroxy oxygen and amino hydrogen, respectively, through hydrogen bonds. Residues SER496 and HIS505 formed H-bonds with ZINC000035399302 hydroxy oxygen. ZINC95910594 formed hydrogen bonds with HIS505 to enhance the binding affinity to RBD. Candidate compounds (ZINC95919448, ZINC85531210, ZINC95610651, and ZINC95910594) were obtained from the TCMNP database. ZINC000035399302 was obtained from the Alinda natural products database. The close interaction between hot residues and the ACE2 protein increases the probability of the Omicron VOC infecting the human body significantly. The top 5 compounds have strong interactions with hotpot residues, which may prevent Omicron RBD from recognizing ACE2. Notably, their large molecular weight results in a high degree of spatial matching with PPI interfaces with large, flat features (Supplementary Table S1). ZINC000035399302, also known as Deoxybouvardin (RA-V), is a cyclopeptide of the Rubiaceae type that was isolated from Bouvardia ternifolia in 1977 (Jolad et al., 1977). Deoxybouvardin and other rubiaceae cyclic peptides are bicyclic cyclic hexapeptides with strong antitumor activity. RA-V is a specific inhibitor of Yes-associated protein (YAP) and transcriptional coactivator with PDZ-binding motif (TAZ) (Ji et al., 2018). RA-V inhibits YAP activation-induced development of liver tumors and is a potential anti-cancer drug. Interestingly, RA-V was found to have anti-inflammatory effects on neutrophil recruitment and edema in carrageenan-induced mouse foot inflammation models (Rupachandra et al., 2020). Moreover, RA-V can also inhibit oxidative stress and expressions of cox as well as TNF-α in inflammatory tissues. RA-V is a potential anti-inflammatory agent. Studies on the other four candidate compounds are not conclusive.

The 10 compounds obtained from the previous screening were subjected to final screening via Molecular mechanics/Poisson-Boltzmann surface area calculations. The binding free energies were further estimated, including four energy sub-items (ΔEvdW, ΔEele, ΔGpolar, and ΔGnon-polar) based on stable Gromacs trajectories.

The calculation of the contribution of entropy to the binding free energy is a challenging project. Most of approaches that can be used to estimate the entropy of a molecule are time-consuming, and the magnitude of the standard error value is high compared to other contributions. On the other hand, multiple studies have shown that the contribution of net entropy is usually small, and the correction for the change in the free energy of the system configuration leads to only a small improvement in the correlation with the experiment (Brown and Muchmore, 2009; Rastelli et al., 2010; Valdés-Tresanco et al., 2021). When the contribution of entropy is ignored, the value obtained is the effective free energy, which is usually sufficient for comparing relative binding free energies of related ligands, for example to compare different ligands binding to the same receptor protein (Wang et al., 2019; Case et al., 2021).

The findings are shown in Figure 3 and Table 2. For all candidate molecules, van der Waals and electrostatic interactions were the main contributors of binding energy and electrostatic interactions were responsible for differences in final binding free energy. The weaker the electrostatic effects, the higher the value of the binding free energy (such as ZINC85542617, ZINC85546719, ZINC85543430, ZINC19764220, ZINC85593889). This also leads to a discrepancy between their binding energy and that of other compounds.

Binding free energy rankings of candidate molecules showed striking agreement with docking results. ZINC95919448, ZINC85531210, ZINC95610651, ZINC000035399302, and ZINC95910594 were the top 5 molecules with the highest binding free energy. ZINC95919448 was not only the best performing candidate molecule in docking but also had the lowest binding free energy in MD simulations. The values of its van der Waals interactions (−43.11 kcal/mol) and its binding free energy (−43.18 kcal/mol) were much higher than those of other compounds. ZINC95919448 was the most potential PPI inhibitor. Interestingly, the chemical structures of ZINC95910594 and ZINC000035399302 have very low similarities (Figure 2), but their binding free energy values and each of its terms are almost the same. ZINC85531210 and ZINC95610651 also showed strong binding abilities to the RBD protein (−22.63 kcal/mol, −23.97 kcal/mol). It implies that they have a comparable degree of spatial and energy matching.

To explore the impact of simulation time on binding free energy, binding free energy fluctuations of top10 compounds during the simulation were calculated (Figure 4 and Supplementary Figure S2). The binding free energy of ZINC95919448 with protein remained in dynamic equilibrium during the simulation. It shows that it has a continuous and stable interaction with the protein during the simulation. In the simulation process, most compounds fluctuate greatly in the early stage and maintain the equilibrium state in the late stage (ZINC95610651, ZINC000035399302, ZINC95910594, ZINC85531210, ZINC85593889, ZINC85542617 and ZINC85546719). This is due to the optimization of the initial conformation of the docking in the early simulation, and when the binding conformation reached the lowest energy state, the protein and ligand maintained a stable interaction. Compounds ZINC19764220 and ZINC85543430 showed the largest fluctuation in binding free energy during the simulation, indicating that they did not have stable interaction with proteins, and that they were not ideal binders.

The combined docking results showed that ZINC95919448, ZINC85531210, ZINC95610651, ZINC000035399302, and ZINC95910594 constitute compounds that bind well to RBD proteins. They were then subjected to the MD analysis and prediction of molecular properties.

The dynamic characteristics of the top5 compounds were calculated during the 100-ns simulation (Figure 5). The root mean square deviations (RMSD) of 5 complexes were analyzed using backbone atoms (Figure 5A). During MD simulation, the RMSD of the five systems fluctuated between 1 and 2 Å. The ZINC95919448 system exhibited slight deviations during 90–100 ns (∼0.3 Å). ZINC85531210 was the system with the smallest fluctuation among the five systems, and the average value of its RMSD was 1.27 Å. The five systems also showed a constant RMSD after 20 ns, implying that they were stable over the entire simulation time and the resulting data was of high confidence. The radius of gyration (Rg) for each frame against the simulated time also showed that the RBD protein that bound with all proposed small molecules attained compactness and rigidity (Figure 5B). The range of Rg values of RBD protein that bound with proposed compounds ZINC95919448, ZINC85531210, ZINC95610651, ZINC000035399302, and ZINC95910594 were 17.9–18.6 Å, 18.0–18.7 Å, 18.0–18.7 Å, 17.9–18.5 Å and 18.0–18.6 Å, respectively. The relatively consistent Rg values also indicated stable folding properties of the protein during the entire MD simulation period.

Fluctuations of protein carbon alpha (Cα) atoms and effects of five candidate molecules binding in the RBD protein were analyzed by root mean square fluctuations plot (Figure 5C). It is reasonable for terminal residues to have high RMSF values. Similar patterns of changes in RMSF values of RBD residues that bound to ZINC95919448 and ZINC000035399302, ZINC85531210 and ZINC95610651 were identified. LEU371 showed the highest RMSF value for ZINC95919448 and ZINC000035399302 systems. LEU371, in a flexible loop and far from the PPI interface, did not affect the results. ASN477 was the most active residue for ZINC95919448 and ZINC000035399302 systems. Although ASN477 is one of the important residues at the RBD-ACE2 PPI interface, it is located at the edge of the RBD protein and did not interact with any of the proposed compounds. Furthermore, other residues (480–505) on the binding surface showed natural fluctuations (1–1.5 Å), indicating that the residues are flexible enough to interact with the ligand.

The RMSD and solvent accessible surface area (SASA) of RBD PPI residues were analyzed (Figures 5D, E). The PPI-RBD RMSD of ZINC000035399302 and ZINC95910594 systems were highly stable which signifies that the residues had the least conformational changes. The convergence of simulation revealed that there were some deviations in RMSD of ZINC95919448, ZINC85531210, and ZINC95610651 systems. The PPI-RBD SASA showed that the ZINC95610651 system was highly exposed to the solvent compared with other complexes (Figure 5E).

The stable RBD-ligand interaction mode during MD simulation is shown in Figure 6. ZINC95919448 formed hydrogen bonds with three residues and mainly maintained hydrophobic interactions with TYR449 and TYR501 (Figure 6A). A hydrogen bond was formed between the hydroxyl oxygen on the six-membered ring of ZINC95919448 and ARG403, consistent with the docking mode. Each of LYU493and SER496 formed a hydrogen bond with the hydroxyl oxygen, which is important for PPI inhibitions. ZINC95610651 formed a new binding conformation with the RBD protein distinct from docking (Figure 6B). ZINC95610651 formed two hydrogen bonds with TYR473 via a hydrogen bond donor and acceptor, respectively. In addition, it formed hydrogen bonds with ASN417 and LEU455. ZINC95910594 formed a close interaction with the PPI interface residues of RBD protein (Figure 6C). TYR449, ASN450, LEU492 and SER494 were involved in the formation of hydrogen bonds. ZINC000035399302 maintained the hydrogen bond with SER496 in the docked conformation and added a hydrogen bond with TYR453, which are necessary for PPI inhibition (Figure 6D). Due to the hydrophobic scaffold of ZINC85531210, it mainly formed hydrophobic interactions with RBD protein, including ASN450, TYR451, ILE468, THR470, and LEU492 among others (Figure 6E).

As expected, the residues mentioned above that formed interactions with candidate compounds were also found to be significant contributors to the binding (Figure 7, Supplementary Tables S2–S6). Residues TYR501 and HIS505 stabilized small molecules through van der Waals interactions in five systems. ARG403 contributed −1.02 kcal/mol to the binding energy, further demonstrating its important role in ZINC95919448 binding. Moreover, the binding contributions of TYR449, SER494, SER496, and ARG498 to ZINC95919448 and RBD protein were also non-negligible. LEU455, PHE456, TYR473, SER496 played important roles in ZINC95610651 binding to RBD, consistent with findings from the interaction analysis. For ZINC95910594, the contribution of TYR449 was very significant (−3.55 kcal/mol), which may be due to hydrophobic interaction and hydrogen bonding. Residues ARG403, TYR453, and SER496 contributed to the binding of ZINC000035399302 to RBD. For ZINC85531210, TYR449, LYS493, and SER494 had a significant destabilizing effect. Therefore, hotspot residues and other recognition residues were widely involved in the binding of candidate small molecules. These are consistent with the analysis results of stable binding conformation.

From the above interaction analyses and binding energy decomposition analysis, we concluded that (1) more stable and reasonable conformation of the interaction were obtained in MD simulation, suggesting that MD simulations play a role in the refinement of the docking results, (2) the top 5 compounds are the efficient RBD protein binding agents, and (3) all these compounds have stable interactions with the unique mutant residues of Omicron VOC, so they might prevent Omicron from recognizing ACE2.

The pharmacokinetic and side effects for the five proposed compounds are presented in Table 3. The blood-brain barrier (BBB) prevents dangerous substances from entering the brain. However, the BBB can block some drugs from entering the brain and the central nervous system. ZINC95919448, ZINC95610651, and ZINC85531210 could not penetrate the BBB, while ZINC85531210 was shown to penetrate the BBB. Apart from ZINC95910594, all other compounds were absorbed through the gastrointestinal tract. P-glycoprotein (P-gp) is an important protein of the cell membrane that pumps many foreign substances out of cells. ZINC95610651 was not shown to be a substrate for P-gp while ZINC000035399302 was a P-gp inhibitor. All candidate molecules showed very low inhibitions of cytochrome P450 enzymes (CYP). ZINC000035399302 had a good performance with regards to Caco-2 cell permeability exhibiting the permeability of molecules into the large intestines. ZINC95919448 and ZINC95910594 do not cause severe drug-induced liver injury.

The maxima and minima of candidate molecular ESPs were localized on the vdW surface by quantitative molecular surface analysis (Figures 8A–E). In addition, the distributions of different ESP intervals on the vdW surface were calculated (Manzetti and Lu, 2013; Lu and Manzetti, 2014) and plotted (Figures 8F–J). With regards to color scale in the right side of the figure, negative values (blue) display the electron-rich negatively charged part of the molecule, while positive values (red) indicate opposite characteristics. The ubiquitous presence of red and blue in the ESP map demonstrates the polar character of the molecule, with the blue part of the map more likely to interact with the positively charged residues or the red part with the negatively charged residues. The median value of ESP (white) represents the weakly polar hydrophobic part of the molecule. As shown in from the labelled text in Figures 8A–E, the ESP of the vdW surface of ZINC95919448 ranges from −45.22–71.74 kcal/mol, which is a wider range relative to the other four molecules, indicating that ZINC95919448 might form strong electrostatic interactions with amino acid residues. ZINC95919448 and ZINC95910594 tended to act as a Lewis acid to bound Lewis base species. ZINC000035399302 and ZINC85531210 tended to act more as a Lewis base to dock to the Lewis acid species. ZINC95910594 had no obvious tendency. This was based on comparisons of the magnitude of ESP at the maximum and minimum on the vdW surface as well as on comparisons of the surface area of positive and negative ESPs. The electrostatic potential is one of the fundamentals of electrostatic interaction. Because electrostatic interaction is the primary long-distance interaction between molecules. Electrostatic potential plays a very unique role in the interaction of inter-molecules and protein, the reaction site of molecules in metabolism, and molecular recognition.

The following descriptors were obtained: chemical hardness (η), chemical softness (S), band gap (GAPE), electron affinity (EA), ionization potential (IP), electrophilicity index (ω), nucleophilicity index (N), electronic potential (μ), and electronegativity (χ) (Table 4). Based on Koopman’s theorem (LIU, 2016), opposite numbers of the values of HOMOs and LUMOs correspond to the electron affinity (EA) and ionization potential (IP), respectively. Bases on the Pearson’s Hard/Soft Acid-Base (HSAB) principle, Chemicals can be described as hard or soft acids or bases (Parr and Pearson, 1983). Hard species are small in volume, difficult to polarize, highly charged, and have a large HOMO-LUMO gap, while soft species are the opposite. ZINC95919448 has lower ionization potential (IP = 5.16 eV), electron affinity (EA = −0.46 eV) and chemical hardness (η = 6.9 eV) values, which support its reactivity. Nucleophilicity is an important parameter that can reflect the reactivity of compounds. ZINC95919448 has strong nucleophilicity (N = 3.96 eV), while ZINC85531210 has soft nucleophilicity (N = 3.23 eV). Quantum chemical descriptors are of great significance for evaluating the molecular properties and guiding experiments.