Interactive Mechanism of Potential Inhibitors with Glycosyl for SARS-CoV-2 by Molecular Dynamics Simulation

: Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is a type of Ribonucleic Acid (RNA) coronavirus and it has infected and killed many people around the world. It is reported that the receptor binding domain of the spike protein (S_RBD) of the SARS-CoV-2 virus is responsible for attachment to human angiotensin converting enzyme II (ACE2). Many researchers are attempting to search potential inhibitors for ﬁghting SARS-CoV-2 infection using theoretical or experimental methods. In terms of experimental and theoretical research, Cefuroxime, Erythromycin, Lincomycin and Oﬂoxacin are the potential inhibitors of SARS-CoV-2. However, the interactive mechanism of the protein SARS-CoV-2 and the inhibitors are still elusive. Here, we investigated the interactions between S_RBD and the inhibitors using molecular dynamics (MD) simulations. Interestingly, we found that there are two binding sites of S_RBD for the four small molecules. In addition, our analysis also illustrated that hydrophobic and π - π stacking interactions play crucial roles in the interactions between S_RBD and the small molecules. In our work, we also found that small molecules with glycosyl group have more effect on the conformation of S_RBD than other inhibitors, and they are also potential inhibitors for the genetic variants of SARS-CoV-2. This study provides in silico-derived mechanistic insights into the interactions of S_RBD and inhibitors, which may provide new clues for ﬁghting SARS-CoV-2 infection.

One experiment [13] determined the crystal structure of SARS-CoV-2 and described that the receptor binding domain of the spike protein (S_RBD) of the SARS-CoV-2 virus is responsible for attachment to human angiotensin converting enzyme II (ACE2). ACE2, as a cell receptor for SARS-CoV-2, will assist SARS-CoV-2 entry into cells and cause the virus to infect humans [13]. According to the experimental report, S_RBD has a greater binding affinity for ACE2 [14]. A number of recent articles have studied monoclonal inhibitors [15,16] of the SARS-CoV-2 virus. On the basis of the trial seventh edition of diagnosis and treatment guideline of COVID-19 issued by National Health Commission of

Molecular Docking Strategy
To screen the potential inhibitors of SARS-CoV-2, we employed semi-flexible docking in AutoDock software [31][32][33]. The grid spacing was kept at 0.375 Å and the size of grid was 92 Å × 120 Å × 42 Å. The center coordinates of the box were (−36.767, 28 28.064, 1.259) in x, y and z dimensions, for the S_RBD with Cefuroxime system, S_RBD with Erythromycin system, S_RBD with Lincomycin system, and S_RBD with Ofloxacin system, respectively.

Molecular Docking Strategy
To screen the potential inhibitors of SARS-CoV-2, we employed semi-flexible docking in AutoDock software [31][32][33]. The grid spacing was kept at 0.375 Å and the size of grid was 92 Å * 120 Å * 42 Å. The center coordinates of the box were (−36.767, 28 28.064, 1.259) in x, y and z dimensions, for the S_RBD with Cefuroxime system, S_RBD with Erythromycin system, S_RBD with Lincomycin system, and S_RBD with Ofloxacin system, respectively.
For predicted effective small molecules, the grid spacing was kept at 0.308 Å and the size of grid was 56 Å * 46 Å * 56 Å in domain 1. The coordinates of the center grid box were (−32.221, 15.724, and 28.386) in x, y, and z dimensions. However, for m_S_RBD protein and small molecules, all parameters are the same as the S_RBD with Lincomycin system.
Here we took the structure of complex with the lowest binding energy as the initial structure, as shown in Tables 1, S2 and S3.

MD Simulation
The simulation time was 400 ns and the temperature was 310 K in each MD trajectory for all simulation systems with AMBER99SB-ILDN force field [34]. To imitate neutral pH conditions, the residues of Asp, Glu, Arg, and Lys were charged (Asp −1 , Glu −1 , Arg +1 , and Lys +1 ). The N-and C-termini were also charged (NH 3+ and COO -) and counterions (Cl − , Na + ) were added to neutralize these systems. Eight Na + ions were added in the S_RBD with ACE2 system and the other systems were neutralized with two Cl − ions. All MD simulations were performed in the GROMACS-2018.1 software package [35]. The details of all the simulations are summarized in Table S1. The visual inspection of four systems was carefully performed using visual molecular dynamics (VMD) [36] and PyMol software [37].

Analysis Methods
First, we discarded the first 300 ns MD trajectory with the purpose of removing the bias of the initial states. We performed the analysis using our in-house developed codes and the GROMACS facilities [35]. The structural analysis of the conformational ensemble was performed by evaluating root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), the solvent accessible surface area (SASA) and the contact surface area (CSA). Here, we identified that the residues were in contact state when the minimum distance between residue-residue of S_RBD was within 0.54 nm. We identified that two rings formed a π-π interaction if the centroid distance between two aromatic residues was detected to be within 0.65 nm [38]. We also calculated the 2D free energy landscape, applying the formula, RTlnH (x, y), where H (x,y) is the histogram of two selected reaction coordinates, which is the angle and distance.

Molecular Docking Strategy
To screen the potential inhibitors of SARS-CoV-2, we employed semi-flexible docking in AutoDock software [31][32][33] 1.259) in x, y and z dimensions, for the S_RBD with Cefuroxime system, S_RBD with Erythromycin system, S_RBD with Lincomycin system, and S_RBD with Ofloxacin system, respectively.
For predicted effective small molecules, the grid spacing was kept at 0.308 Å and the size of grid was 56 Å * 46 Å * 56 Å in domain 1. The coordinates of the center grid box were (−32.221, 15.724, and 28.386) in x, y, and z dimensions. However, for m_S_RBD protein and small molecules, all parameters are the same as the S_RBD with Lincomycin system.
Here we took the structure of complex with the lowest binding energy as the initial structure, as shown in Tables 1, S2 and S3.

MD Simulation
The simulation time was 400 ns and the temperature was 310 K in each MD trajectory for all simulation systems with AMBER99SB-ILDN force field [34]. To imitate neutral pH conditions, the residues of Asp, Glu, Arg, and Lys were charged (Asp −1 , Glu −1 , Arg +1 , and Lys +1 ). The N-and C-termini were also charged (NH 3+ and COO -) and counterions (Cl − , Na + ) were added to neutralize these systems. Eight Na + ions were added in the S_RBD with ACE2 system and the other systems were neutralized with two Cl − ions. All MD simulations were performed in the GROMACS-2018.1 software package [35]. The details of all the simulations are summarized in Table S1. The visual inspection of four systems was carefully performed using visual molecular dynamics (VMD) [36] and PyMol software [37].

Analysis Methods
First, we discarded the first 300 ns MD trajectory with the purpose of removing the bias of the initial states. We performed the analysis using our in-house developed codes and the GROMACS facilities [35]. The structural analysis of the conformational ensemble was performed by evaluating root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), the solvent accessible surface area (SASA) and the contact surface area (CSA). Here, we identified that the residues were in contact state when the minimum distance between residue-residue of S_RBD was within 0.54 nm. We identified that two rings formed a π-π interaction if the centroid distance between two aromatic residues was detected to be within 0.65 nm [38]. We also calculated the 2D free energy landscape, applying the formula, RTlnH (x, y), where H (x,y) is the histogram of two selected reaction coordinates, which is the angle and distance.

Molecular Docking Strategy
To screen the potential inhibitors of SARS-CoV-2, we employed semi-flexible docking in AutoDock software [31][32][33]. The grid spacing was kept at 0.375 Å and the size of grid was 92 Å * 120 Å * 42 Å. The center coordinates of the box were (−36.767, 28 28.064, 1.259) in x, y and z dimensions, for the S_RBD with Cefuroxime system, S_RBD with Erythromycin system, S_RBD with Lincomycin system, and S_RBD with Ofloxacin system, respectively.
For predicted effective small molecules, the grid spacing was kept at 0.308 Å and the size of grid was 56 Å * 46 Å * 56 Å in domain 1. The coordinates of the center grid box were (−32.221, 15.724, and 28.386) in x, y, and z dimensions. However, for m_S_RBD protein and small molecules, all parameters are the same as the S_RBD with Lincomycin system.
Here we took the structure of complex with the lowest binding energy as the initial structure, as shown in Tables 1, S2 and S3.

MD Simulation
The simulation time was 400 ns and the temperature was 310 K in each MD trajectory for all simulation systems with AMBER99SB-ILDN force field [34]. To imitate neutral pH conditions, the residues of Asp, Glu, Arg, and Lys were charged (Asp −1 , Glu −1 , Arg +1 , and Lys +1 ). The N-and C-termini were also charged (NH 3+ and COO -) and counterions (Cl − , Na + ) were added to neutralize these systems. Eight Na + ions were added in the S_RBD with ACE2 system and the other systems were neutralized with two Cl − ions. All MD simulations were performed in the GROMACS-2018.1 software package [35]. The details of all the simulations are summarized in Table S1. The visual inspection of four systems was carefully performed using visual molecular dynamics (VMD) [36] and PyMol software [37].

Analysis Methods
First, we discarded the first 300 ns MD trajectory with the purpose of removing the bias of the initial states. We performed the analysis using our in-house developed codes and the GROMACS facilities [35]. The structural analysis of the conformational ensemble was performed by evaluating root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), the solvent accessible surface area (SASA) and the contact surface area (CSA). Here, we identified that the residues were in contact state when the minimum distance between residue-residue of S_RBD was within 0.54 nm. We identified that two rings formed a π-π interaction if the centroid distance between two aromatic residues was detected to be within 0.65 nm [38]. We also calculated the 2D free energy landscape, applying the formula, RTlnH (x, y), where H (x,y) is the histogram of two selected reaction coordinates, which is the angle and distance.

Molecular Docking Strategy
To screen the potential inhibitors of SARS-CoV-2, we employed semi-flexible docking in AutoDock software [31][32][33]. The grid spacing was kept at 0.375 Å and the size of grid was 92 Å * 120 Å * 42 Å. The center coordinates of the box were (−36.767, 28 28.064, 1.259) in x, y and z dimensions, for the S_RBD with Cefuroxime system, S_RBD with Erythromycin system, S_RBD with Lincomycin system, and S_RBD with Ofloxacin system, respectively.
For predicted effective small molecules, the grid spacing was kept at 0.308 Å and the size of grid was 56 Å * 46 Å * 56 Å in domain 1. The coordinates of the center grid box were (−32.221, 15.724, and 28.386) in x, y, and z dimensions. However, for m_S_RBD protein and small molecules, all parameters are the same as the S_RBD with Lincomycin system.
Here we took the structure of complex with the lowest binding energy as the initial structure, as shown in Tables 1, S2 and S3.

MD Simulation
The simulation time was 400 ns and the temperature was 310 K in each MD trajectory for all simulation systems with AMBER99SB-ILDN force field [34]. To imitate neutral pH conditions, the residues of Asp, Glu, Arg, and Lys were charged (Asp −1 , Glu −1 , Arg +1 , and Lys +1 ). The N-and C-termini were also charged (NH 3+ and COO -) and counterions (Cl − , Na + ) were added to neutralize these systems. Eight Na + ions were added in the S_RBD with ACE2 system and the other systems were neutralized with two Cl − ions. All MD simulations were performed in the GROMACS-2018.1 software package [35]. The details of all the simulations are summarized in Table S1. The visual inspection of four systems was carefully performed using visual molecular dynamics (VMD) [36] and PyMol software [37].

Analysis Methods
First, we discarded the first 300 ns MD trajectory with the purpose of removing the bias of the initial states. We performed the analysis using our in-house developed codes and the GROMACS facilities [35]. The structural analysis of the conformational ensemble was performed by evaluating root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), the solvent accessible surface area (SASA) and the contact surface area (CSA). Here, we identified that the residues were in contact state when the minimum distance between residue-residue of S_RBD was within 0.54 nm. We identified that two rings formed a π-π interaction if the centroid distance between two aromatic residues was detected to be within 0.65 nm [38]. We also calculated the 2D free energy landscape, applying the formula, RTlnH (x, y), where H (x,y) is the histogram of two selected reaction coordinates, which is the angle and distance.
For predicted effective small molecules, the grid spacing was kept at 0.308 Å and the size of grid was 56 Å × 46 Å × 56 Å in domain 1. The coordinates of the center grid box were (−32.221, 15.724, and 28.386) in x, y, and z dimensions. However, for m_S_RBD protein and small molecules, all parameters are the same as the S_RBD with Lincomycin system.
Here we took the structure of complex with the lowest binding energy as the initial structure, as shown in Table 1, Tables S2 and S3.

MD Simulation
The simulation time was 400 ns and the temperature was 310 K in each MD trajectory for all simulation systems with AMBER99SB-ILDN force field [34]. To imitate neutral pH conditions, the residues of Asp, Glu, Arg, and Lys were charged (Asp −1 , Glu −1 , Arg +1 , and Lys +1 ). The N-and C-termini were also charged (NH 3+ and COO -) and counterions (Cl − , Na + ) were added to neutralize these systems. Eight Na + ions were added in the S_RBD with ACE2 system and the other systems were neutralized with two Cl − ions. All MD simulations were performed in the GROMACS-2018.1 software package [35]. The details of all the simulations are summarized in Table S1. The visual inspection of four systems was carefully performed using visual molecular dynamics (VMD) [36] and PyMol software [37].

Analysis Methods
First, we discarded the first 300 ns MD trajectory with the purpose of removing the bias of the initial states. We performed the analysis using our in-house developed codes and the GROMACS facilities [35]. The structural analysis of the conformational ensemble was performed by evaluating root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), the solvent accessible surface area (SASA) and the contact surface area (CSA). Here, we identified that the residues were in contact state when the minimum distance between residue-residue of S_RBD was within 0.54 nm. We identified that two rings formed a π-π interaction if the centroid distance between two aromatic residues was detected to be within 0.65 nm [38]. We also calculated the 2D free energy landscape, applying the formula, RTlnH (x, y), where H (x, y) is the histogram of two selected reaction coordinates, which is the angle and distance.
Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) method is employed to calculate the binding energies between receptor and ligand, implemented in the GROMACS package [39][40][41]. In this method, the binding free energy (∆G binding ) satisfies the formula ∆G binding = ∆E MM + ∆G solv − T∆S, where ∆E MM is the gas phase energy. ∆EMM satisfies the formula ∆E MM = ∆E bonded + ∆E non-bonded , where ∆E bonded is bonded interactions taken as 0 and ∆E non-bonded consists of electrostatic (∆E elec ) and van der Waals (∆E vdW ) terms [31,42], hence ∆E MM ought to fulfil the formula ∆E MM = ∆E elec + ∆E vdW . ∆G solv represents the solvation energy, consisting of the polar solvation energy, ∆G polar , and the nonpolar solvation component ∆G surf . On the basis of numerous previous computational studies [43][44][45][46], the contribution of conformational entropy of the peptides was ignored. Therefore, in this study, we use the formula ∆G binding = ∆E MM + ∆G solv [47] to calculate the binding free energy.

Results and Discussion
According to the relevant literature and drug properties, we selected four inhibitors with antibacterial or bactericidal effects: Cefuroxime, Erythromycin, Lincomycin, and Ofloxacin from Cephalosporins, Macrolides, Lincosamides, and Quinolone, respectively [48][49][50][51][52]. To understand the interactional mechanisms of the inhibitors with S_RBD, we obtained the structures of S_RBD with small molecules/ACE2 by AutoDock method and performed molecular dynamics simulations for the six systems, the S_RBD system, the S_RBD in the presence of ACE2 system, the S_RBD in the presence of Cefuroxime system, the S_RBD in the presence of Erythromycin system, the S_RBD in the presence of Lincomycin system and the S_RBD in the presence of Ofloxacin system. Table 1 illustrates molecular formulas, binding energies, and chemical structures of the inhibitors. All of the initial states of S_RBD and small molecules/ACE2 are presented using PyMol, as shown in Figures S1 and S2.

The Convergence and Stability of S_RBD in the Research Systems
In order to check the convergence and compare the structural stability of S_RBD, we conducted the time evolution of backbone RMSD analysis in each research system. Figure 1 showed that the RMSD values of S_RBD protein increased rapidly from the initial value (0 nm) and then fluctuates in the six systems. It is obvious that the RMSD values of S_RBD fluctuate around 0.2/0.3 nm during the last 100 ns trajectories. It is shown that our simulation trajectories are converged in all of the six systems. Overall, the analysis of the RMSD shows that our MD simulations have reasonably converged and the conformation of S_RBD is stable in the absence or presence of the small molecules/ACE2. solvation energy, consisting of the polar solvation energy, ΔGpolar, and the nonpolar solvation component ΔGsurf. On the basis of numerous previous computational studies [43][44][45][46], the contribution of conformational entropy of the peptides was ignored. Therefore, in this study, we use the formula ΔGbinding = ΔEMM + ΔGsolv [47] to calculate the binding free energy.

Results and Discussion
According to the relevant literature and drug properties, we selected four inhibitors with antibacterial or bactericidal effects: Cefuroxime, Erythromycin, Lincomycin, and Ofloxacin from Cephalosporins, Macrolides, Lincosamides, and Quinolone, respectively [48][49][50][51][52]. To understand the interactional mechanisms of the inhibitors with S_RBD, we obtained the structures of S_RBD with small molecules/ACE2 by AutoDock method and performed molecular dynamics simulations for the six systems, the S_RBD system, the S_RBD in the presence of ACE2 system, the S_RBD in the presence of Cefuroxime system, the S_RBD in the presence of Erythromycin system, the S_RBD in the presence of Lincomycin system and the S_RBD in the presence of Ofloxacin system. Table 1 illustrates molecular formulas, binding energies, and chemical structures of the inhibitors. All of the initial states of S_RBD and small molecules/ACE2 are presented using PyMol, as shown in Figures S1 and S2.

The Convergence and Stability of S_RBD in the Research Systems
In order to check the convergence and compare the structural stability of S_RBD, we conducted the time evolution of backbone RMSD analysis in each research system. Figure  1 showed that the RMSD values of S_RBD protein increased rapidly from the initial value (0 nm) and then fluctuates in the six systems. It is obvious that the RMSD values of S_RBD fluctuate around 0.2/0.3 nm during the last 100 ns trajectories. It is shown that our simulation trajectories are converged in all of the six systems. Overall, the analysis of the RMSD shows that our MD simulations have reasonably converged and the conformation of S_RBD is stable in the absence or presence of the small molecules/ACE2.
To describe the conformational characteristics of S_RBD in the absence or presence of the inhibitors/ACE2, we presented the sequence of S_RBD in Figure 2A. The core domain of S_RBD consists of four α-helixes (α1: residues F338-N343, α2: residues V366-N370, α3: residues V407-Q409, α4: residues K417-N422) and seven β-sheets (β1: residues N354-I358, β2: residues S375-Y380, β3: residues T393-R403, β4: residues G431-N437, β5: residues L452-R454, β6: residues P492-S494, β7: residues P507-E516), as shown in Figure  2A,C. We also define domain 1 (the residues P337-F374 and A435-D442) and domain 2 (the residues between I472 and F490) in Figure 2. The top view and side view of binding region are shown in Figure S3 for the S_RBD protein. To observe the binding sites of the S_RBD protein and inhibitors, we plotted the snapshots of the four complex systems (S_RBD with Cefuroxime, Erythromycin, Lincomycin system, and Ofloxacin systems) in 0 ns, 100 ns, 300 ns, and 400 ns by PyMol, as shown in Figure 2D-G. We found that the position of small molecules in the complex systems tend to be stable after 300 ns. It is obvious that the small molecules of Cefuroxime and Erythromycin share the same binding region of S_RBD, near the random turn in domain 2, as shown in Figure 2D,E. However, Lincomycin and Ofloxacin share the same binding region of S_RBD, among two α-helices in domain 1, as shown in Figure 2F,G. Our data show that there are two binding regions of S_RBD in our research systems.
In order to detect the flexibility of each residue of S_RBD, we calculated the value of RMSF using the last 100 ns trajectories of all the simulations for each system in Figure 2B. Dramatically, we found that the residues in the S_RBD system without inhibitors/ACE2 have greater flexibility than the S_RBD protein systems in the presence of small molecules and ACE2. Interestingly, we also found that the residues in domain 1 are more flexible in To observe the binding sites of the S_RBD protein and inhibitors, we plotted the snapshots of the four complex systems (S_RBD with Cefuroxime, Erythromycin, Lincomycin system, and Ofloxacin systems) in 0 ns, 100 ns, 300 ns, and 400 ns by PyMol, as shown in Figure 2D-G. We found that the position of small molecules in the complex systems tend to be stable after 300 ns. It is obvious that the small molecules of Cefuroxime and Erythromycin share the same binding region of S_RBD, near the random turn in domain 2, as shown in Figure 2D,E. However, Lincomycin and Ofloxacin share the same binding region of S_RBD, among two α-helices in domain 1, as shown in Figure 2F,G. Our data show that there are two binding regions of S_RBD in our research systems.
In order to detect the flexibility of each residue of S_RBD, we calculated the value of RMSF using the last 100 ns trajectories of all the simulations for each system in Figure 2B. Dramatically, we found that the residues in the S_RBD system without inhibitors/ACE2 have greater flexibility than the S_RBD protein systems in the presence of small molecules and ACE2. Interestingly, we also found that the residues in domain 1 are more flexible in the S_RBD in the presence of Cefuroxime and Erythromycin ACE2 systems than in the S_RBD in the presence of Lincomycin and Ofloxacin systems, such as the residues F342, K356, S371, and N440 (Table S4).
For the purpose of further researching the conformational difference of S_RBD in six systems, we presented the alignment of the initial conformation (grey) and the representative snapshot at 400 ns in each simulation with the new cartoon method. The details are shown in Figure S4A-F. Interestingly, the residues of S_RBD in domain 1 and domain 2 are Processes 2021, 9, 1749 6 of 15 larger structural differences in the S_RBD system and the S_RBD with ACE2 system than other research systems.
From the above data, we can summarize that the residues of S_RBD in the absence of small molecules/ACE2 are greater flexible than those in the presence of the small molecules/ACE2. The residues of S_RBD are the most stable in the presence of Lincomycin among all the research systems.

Analysis of the Interactions among the Residue-Residue of S_RBD
To further understand the residue-residue interactions of S_RBD among different systems, we plotted the contact probability between each pair of residues in the absence and the presence of the inhibitors/ACE2. Figure 3 reflects the interaction of the protein itself in the absence and the presence of the inhibitors/ACE2, because different small molecules can have different effects on the interactions between residues of the S_RBD protein. Therefore, the difference of interaction can be compared and the differences between these systems are important. The residue-residue interaction patterns of the six systems are roughly the same, but slightly different, in Figure 3A-F, and the difference of the interactions has been marked. In the S_RBD system, a peptide-peptide interaction occurs between residues T333-W353 and residues W353-C379 around domain 1. However, the contacts between these residues become strong in the presence of the inhibitors and ACE2 in Figure 3B-F. In addition, we found most of these residue pairs display increased contact probability in the four systems by comparing the S_RBD protein system and S_RBD in the presence of ACE system. Interestingly, for the complex system and the S_RBD in the presence of ACE2 system, contact probabilities between residue E516-C525 and residue T478-F486 (around domain 2) notably increase in the presence of small molecules of Cefuroxime/Erythromycin/Lincomycin/Ofloxacin. This phenomenon shows that the potential inhibitors enhance the interaction of S_RBD. In short, our data shows that the residues in domain 1 and domain 2 have stronger interactions than other regions of the S_RBD, and this phenomenon is also consistent with the previous analysis of RMSF. We can propose that the two binding sites, domain 1 and domain 2, play critical roles in molecular progress.
the S_RBD in the presence of Cefuroxime and Erythromycin ACE2 systems than in the S_RBD in the presence of Lincomycin and Ofloxacin systems, such as the residues F342, K356, S371, and N440 (Table S4).
For the purpose of further researching the conformational difference of S_RBD in six systems, we presented the alignment of the initial conformation (grey) and the representative snapshot at 400 ns in each simulation with the new cartoon method. The details are shown in Figure S4A-F. Interestingly, the residues of S_RBD in domain 1 and domain 2 are larger structural differences in the S_RBD system and the S_RBD with ACE2 system than other research systems.
From the above data, we can summarize that the residues of S_RBD in the absence of small molecules/ACE2 are greater flexible than those in the presence of the small molecules/ACE2. The residues of S_RBD are the most stable in the presence of Lincomycin among all the research systems.

Analysis of the Interactions among the Residue-Residue of S_RBD
To further understand the residue-residue interactions of S_RBD among different systems, we plotted the contact probability between each pair of residues in the absence and the presence of the inhibitors/ACE2. Figure 3 reflects the interaction of the protein itself in the absence and the presence of the inhibitors/ACE2, because different small molecules can have different effects on the interactions between residues of the S_RBD protein. Therefore, the difference of interaction can be compared and the differences between these systems are important. The residue-residue interaction patterns of the six systems are roughly the same, but slightly different, in Figure 3A-F, and the difference of the interactions has been marked. In the S_RBD system, a peptide-peptide interaction occurs between residues T333-W353 and residues W353-C379 around domain 1. However, the contacts between these residues become strong in the presence of the inhibitors and ACE2 in Figure 3B-F. In addition, we found most of these residue pairs display increased contact probability in the four systems by comparing the S_RBD protein system and S_RBD in the presence of ACE system. Interestingly, for the complex system and the S_RBD in the presence of ACE2 system, contact probabilities between residue E516-C525 and residue T478-F486 (around domain 2) notably increase in the presence of small molecules of Cefuroxime/Erythromycin/Lincomycin/Ofloxacin. This phenomenon shows that the potential inhibitors enhance the interaction of S_RBD. In short, our data shows that the residues in domain 1 and domain 2 have stronger interactions than other regions of the S_RBD, and this phenomenon is also consistent with the previous analysis of RMSF. We can propose that the two binding sites, domain 1 and domain 2, play critical roles in molecular progress.

The Conformational Difference of S_RBD in Different Research Systems
To further understand the conformational difference of S_RBD with different small molecules and ACE2, we plotted the probability distribution function (PDF) of the Rg, SASA and CSA, as shown in Figure 4A-C. Figure 4A shows the probability of the Rg. According to Figure 4A, in the S_RBD with Lincomycin system, the radius of gyration is 1.820 nm, which has an extremely large peak. The peaks are extremely reduced in S_RBD with Cefuroxime system (1.824 nm), S_RBD with Erythromycin system (1.846 nm), S_RBD with Ofloxacin system (1.840 nm), S_RBD with ACE2 system (1.839 nm), and S_RBD system (1.833 nm). In addition, the SASA distribution curves of the S_RBD show that there exists a sharp peak at a SASA value of 100 nm 2 in the S_RBD system, S_RBD in the presence of Cefuroxime system and S_RBD in the presence of Lincomycin system, while in the S_RBD with Erythromycin, Ofloxacin system and ACE2 system, this peak of SASA value is shifted to approximately 106 nm 2 ( Figure 4B). In a word, the values of Rg and SASA indicate that conformation of S_RBD in the S_RBD in the presence of Lincomycin system is more compact than that with/without the small molecules Cefuroxime/Erythromycin/Ofloxacin and protein ACE2. S_RBD with Erythromycin, Ofloxacin system and ACE2 system, this peak of SASA value is shifted to approximately 106 nm 2 ( Figure 4B). In a word, the values of Rg and SASA indicate that conformation of S_RBD in the S_RBD in the presence of Lincomycin system is more compact than that with/without the small molecules Cefuroxime/Erythromycin/Ofloxacin and protein ACE2.
In order to further explore the influence of the inhibitors on the conformation of S_RBD, we calculated the distribution of the CSA between S_RBD and the small molecules. Our analysis shows that the average CSA values are 489.48 Å 2 in the S_RBD in the presence of Cefuroxime system, 617.884 Å 2 in the S_RBD in the presence of Erythromycin system, 560.749 Å 2 in the S_RBD in the presence of Lincomycin system, and 546.891 Å 2 in the S_RBD with Ofloxacin system, respectively, as described in Figure 4C. These demonstrate a comparatively stronger binding affinity between S_RBD and Erythromycin/Lincomycin than the other systems [53]. The CSA value of Erythromycin and S_RBD is even larger, because the molecular weight and the surface of Erythromycin is larger than other molecules. Except Erythromycin, Lincomycin has larger CSA value than other two molecules. Our data in Figure 4 agree well with the above analysis in Figure 2. (C) the contact surface area (CSA) in each MD system. The S_RBD with Cefuroxime system (black), S_RBD with Erythromycin system (red), S_RBD with Lincomycin system (blue), S_RBD with Ofloxacin system (green), S_RBD with ACE2 system (purple) and S_RBD system (orange).

Identify the Major Interaction Residues of S_RBD by Calculating the Binding Energy
After determining the effect of the inhibitors on the conformation of the S_RBD protein, we examined the binding free energy of the complex with the MM/PBSA method using the last 100 ns of each MD trajectory. In Table S3, we listed the binding energy (less than 1.0 KJ/mol) between each residue of S_RBD and the small molecules Cefuroxime/Erythromycin/Lincomycin/Ofloxacin. Interestingly, we found the binding energies between S_RBD and small molecules were −36.304 KJ/mol (Cefuroxime), −63.106 KJ/mol (Erythromycin), −62.862 KJ/mol (Lincomycin) and −69.961 KJ/mol (Ofloxacin), respectively, as shown in Table S5. To further explore the underlying interactive mechanism between the protein and the inhibitors, we first calculated the binding free energy of the inhibitors with each amino acid residue of S_RBD by discarding the first 300 ns of each MD trajectory. Figure 5A and Table S6 shows that the small molecule Cefuroxime has a Figure 4. The physical and conformational properties of the six systems. The probability distribution function (PDF) of (A) the radius of gyration (Rg), (B) the solvent accessible surface area (SASA) and (C) the contact surface area (CSA) in each MD system. The S_RBD with Cefuroxime system (black), S_RBD with Erythromycin system (red), S_RBD with Lincomycin system (blue), S_RBD with Ofloxacin system (green), S_RBD with ACE2 system (purple) and S_RBD system (orange).
In order to further explore the influence of the inhibitors on the conformation of S_RBD, we calculated the distribution of the CSA between S_RBD and the small molecules. Our analysis shows that the average CSA values are 489.48 Å 2 in the S_RBD in the presence of Cefuroxime system, 617.884 Å 2 in the S_RBD in the presence of Erythromycin system, 560.749 Å 2 in the S_RBD in the presence of Lincomycin system, and 546.891 Å 2 in the S_RBD with Ofloxacin system, respectively, as described in Figure 4C. These demonstrate a comparatively stronger binding affinity between S_RBD and Erythromycin/Lincomycin than the other systems [53]. The CSA value of Erythromycin and S_RBD is even larger, because the molecular weight and the surface of Erythromycin is larger than other molecules. Except Erythromycin, Lincomycin has larger CSA value than other two molecules. Our data in Figure 4 agree well with the above analysis in Figure 2.

Identify the Major Interaction Residues of S_RBD by Calculating the Binding Energy
After determining the effect of the inhibitors on the conformation of the S_RBD protein, we examined the binding free energy of the complex with the MM/PBSA method using the last 100 ns of each MD trajectory. In Table S3, we listed the binding energy (less than 1.0 KJ/mol) between each residue of S_RBD and the small molecules Cefuroxime/Erythromycin/Lincomycin/Ofloxacin. Interestingly, we found the binding energies between S_RBD and small molecules were −36.304 KJ/mol (Cefuroxime), −63.106 KJ/mol (Erythromycin), −62.862 KJ/mol (Lincomycin) and −69.961 KJ/mol (Ofloxacin), respectively, as shown in Table S5. To further explore the underlying interactive mechanism between the protein and the inhibitors, we first calculated the binding free energy of the inhibitors with each amino acid residue of S_RBD by discarding the first 300 ns of each MD trajectory. Figure 5A and Table S6 shows that the small molecule Cefuroxime has a relatively low binding energy with the hydrophobic residue (P479) and the aromatic residues (Y473 and Y489) in domain 2 for the S_RBD in the presence of Cefuroxime system. We also found that the small molecule Erythromycin has a relatively low binding energy with the hydrophobic residues (L455, F456, and A475) and the aromatic residues (Y473 and Y489) in domain 2 for the S_RBD with Erythromycin system. Some researchers reported the residues (L455, F456, and Y489) that S_RBD interact with human ACE2 protein in  [14,54]. In our research, the binding sites in domain 2 consists with the reports well. From Figure 5E,G and Table S3, it can be seen that major interacting residues of Lincomycin to the S_RBD protein (residues: F342, Y367, L368, F374, and W436) are exceedingly similar to those of Ofloxacin to the protein of S_RBD (residues: F342, F374, and W436) in domain 1. The binding residues between S_RBD protein and small molecules, in domain 1, have also been reported in some simulation calculations [55,56].

Analysis of Physical Interactions between the Residues and the Small Molecules
To better understand the interaction of stacking patterns between benzene-like rings of small molecules and the aromatic rings of S_RBD residues, we analyzed the potential To further have a detailed view of the characteristic interaction of the small molecules and S_RBD, we plotted the most representative binding model of the small molecules to S_RBD in Figure 5B,D,F,H using Discovery Studio. The types of interactions between the small molecule Cefuroxime and the residues belong to conventional hydrogen bond (Gln474 and Gly476), and carbon hydrogen bond (Tyr473 and Pro479) in the S_RBD in the presence of Cefuroxime system. The interactions between the small molecule Erythromycin and the residues belong to conventional hydrogen bond (Leu455), alkyl interaction (Ala475) and pi-alkyl interactions (Phe456, Tyr489, and Tyr473) in the S_RBD in the presence of Erythromycin system. The interactions between the small molecule Lincomycin and the residues belong to alkyl interaction (Val367 and Leu368) and pi-alkyl interactions (Phe342, Phe374, and Trp436) in the S_RBD with Lincomycin system. In addition, the interactions between the small molecule Ofloxacin and the residues also belong to pi-alkyl interactions (Phe342 and Trp436) and π-π stacked interaction (Phe374) in the S_RBD in the presence of Ofloxacin system. From Figure 5, we can also conclude that the four inhibitors have the two different binding regions with S_RBD. There are strong interactions between small molecules Cefuroxime/Erythromycin and S_RBD in domain 2 for S_RBD with Cefuroxime system and S_RBD with Erythromycin system. There are also strong interactions between small molecules Lincomycin/Ofloxacin and S_RBD in domain1 for the S_RBD with Lincomycin system and S_RBD with Ofloxacin system. Hydrophobic interactions and π-π interactions play major roles between the inhibitors and the residues of S_RBD. In previous computer simulation studies, hydrophobic interactions and π-π interactions between small molecules and protein have also been reported [57][58][59][60][61].

Analysis of Physical Interactions between the Residues and the Small Molecules
To better understand the interaction of stacking patterns between benzene-like rings of small molecules and the aromatic rings of S_RBD residues, we analyzed the potential mean force (PMF) (in kcal/mol) as a function of the distance between the mass center of the two rings and the angle among the two rings ( Figure 6A-D). The stacking patterns were mainly classified into three categories: parallel (0 • -30 • ), herringbone (30 • -60 • ), and perpendicular or T-shaped (60 • -90 • ) [31]. As shown in Figure 6A, it is obvious that the function image of free energy surface has two minimum-energy basins centered at (0.51 nm < distance < 0.815 nm and 15.2 • < angle < 80.2 • ). The representative snapshot is shown in Figure 6E and the specific values of angle as well as distance are 70.202 • and 0.515 nm, respectively, which correspond to perpendicular-aligned stacking interaction between aromatic residue TYR 473 and the benzene-like ring of Cefuroxime. There are three minimum-energy basins that the distances range from 0.42 nm to 1.12 nm and the angles range from 5.143 • to 80.143 • in the S_RBD with Erythromycin system. The typical detail image is shown in Figure 6F, which corresponds to herringbone π-π stacking patterns between aromatic residue TYR 473 and the aromatic heterocycle of Erythromycin. From Figure 6C,D, we can clearly see that two minimum-energy basins are plotted on the illustration of the free energy surface, which illustrates that the herringbone and Tshaped π−π stacking patterns are dominant in the S_RBD with Lincomycin system, and the parallel and herringbone π−π stacking patterns are populated in the S_RBD with Ofloxacin system. As shown in the Figure 6G,H, there are the representative snapshots, which corresponds to T-shaped/parallel π-π stacking patterns between the aromatic ring of residue PHE 374 and the benzene-like ring of Lincomycin/Ofloxacin. As a result, we can summarize that the stacking interaction between benzene-like rings of small molecules and the aromatic rings plays a vital significant role in maintaining the stability of the complex. Our previous works also reported the three π−π stacking patterns between protein and small molecules [43,58,59,62].
shown in the Figure 6G,H, there are the representative snapshots, which corresponds to T-shaped/parallel π-π stacking patterns between the aromatic ring of residue PHE 374 and the benzene-like ring of Lincomycin/Ofloxacin. As a result, we can summarize that the stacking interaction between benzene-like rings of small molecules and the aromatic rings plays a vital significant role in maintaining the stability of the complex. Our previous works also reported the three π−π stacking patterns between protein and small molecules [43,58,59,62]. Figure 6. Analysis of aromatic stacking interactions between the aromatic heterocycle of the inhibitors and aromatic residues. The free energy landscape as a function of the centroid distance and the angle between the rings of aromatic residues of S_RBD and the benzene-like rings of inhibitors for (A) TYR473 and Cefuroxime in the S_RBD with Cefuroxime system, (B) TYR473 and Erythromycin in the S_RBD with Erythromycin system, (C) PHE374 and Lincomycin in the S_RBD with Lincomycin system and (D) PHE374 and Ofloxacin in the S_RBD with Ofloxacin system. Representative snapshots using PyMOL showing a perpendicular-aligned stacking orientation between the aromatic ring of TYR473 and the benzene-like ring of Cefuroxime (E) in the S_RBD with Cefuroxime system, a herringbone-aligned stacking orientation between the aromatic ring of TYR473 and the benzene-like ring of Erythromycin (F) in the S_RBD with Erythromycin system, a perpendicularaligned stacking orientation between the aromatic ring of PHE374 and the benzene-like ring of Lincomycin (G) in the S_RBD with Lincomycin system and a parallel-aligned stacking orientation between the aromatic ring of PHE374 and the benzene-like ring of Ofloxacin (H) in the S_RBD with Ofloxacin system.

Prediction Potential Inhibitors with Glycosyl for SARS-CoV-2 and the Genetic Variants of SARS-CoV-2 from Existing Drugs
From the above analysis, we found that Lincomycin has more effect on the conformation of S_RBD than other inhibitors. First, based on the structural character of the small molecules, the Lincomycin is divided into two parts: part I and part II, as shown in Figure  7A. Part I belongs to the pyrrolidine, and part II is the glycosyl group. In order to further explore the crucial part of Lincomycin, we calculated the contact probability of S_RBD protein with each part of Lincomycin, as shown in Figure S5. The contact probability of S_RBD protein with part II is higher than part I of small molecule in domain 1 (Table S7). Figure 6. Analysis of aromatic stacking interactions between the aromatic heterocycle of the inhibitors and aromatic residues. The free energy landscape as a function of the centroid distance and the angle between the rings of aromatic residues of S_RBD and the benzene-like rings of inhibitors for (A) TYR473 and Cefuroxime in the S_RBD with Cefuroxime system, (B) TYR473 and Erythromycin in the S_RBD with Erythromycin system, (C) PHE374 and Lincomycin in the S_RBD with Lincomycin system and (D) PHE374 and Ofloxacin in the S_RBD with Ofloxacin system. Representative snapshots using PyMOL showing a perpendicular-aligned stacking orientation between the aromatic ring of TYR473 and the benzene-like ring of Cefuroxime (E) in the S_RBD with Cefuroxime system, a herringbone-aligned stacking orientation between the aromatic ring of TYR473 and the benzene-like ring of Erythromycin (F) in the S_RBD with Erythromycin system, a perpendicular-aligned stacking orientation between the aromatic ring of PHE374 and the benzene-like ring of Lincomycin (G) in the S_RBD with Lincomycin system and a parallel-aligned stacking orientation between the aromatic ring of PHE374 and the benzene-like ring of Ofloxacin (H) in the S_RBD with Ofloxacin system.

Prediction Potential Inhibitors with Glycosyl for SARS-CoV-2 and the Genetic Variants of SARS-CoV-2 from Existing Drugs
From the above analysis, we found that Lincomycin has more effect on the conformation of S_RBD than other inhibitors. First, based on the structural character of the small molecules, the Lincomycin is divided into two parts: part I and part II, as shown in Figure 7A. Part I belongs to the pyrrolidine, and part II is the glycosyl group. In order to further explore the crucial part of Lincomycin, we calculated the contact probability of S_RBD protein with each part of Lincomycin, as shown in Figure S5. The contact probability of S_RBD protein with part II is higher than part I of small molecule in domain 1 (Table S7). This means that the glycosyl group in part II plays important roles in the interactions between Lincomycin and S_RBD.
According to the reports, the existing drugs, such as Lianhuaqingwen capsules [63], Double Coptis [64], and Ribavirin [65], are effective treatment methods for SARS-CoV-2. Interestingly, we found the main components of these drugs, Salidroside [66] Figure S2. We also calculated the free binding energy of potential inhibitors in domain 1 using AutoDock in Table S2. The binding models between the inhibitors and S_RBD in binding domain 1 are also shown in Figure S7 with the Discovery Studio method. We found that the glycosyl moieties in small molecules interact with hydrophobic residues of S_RBD protein. Among all the potential inhibitors, Phillyrin and Liquiritin have relatively lower docking energy in domain 1. The three-dimensional binding modes between the Phillyrin/Liquiritin and S_RBD are shown in Figure 7D,E. From Figure S7, we can see that hydrophobic residues are important in the interactions between S_RBD and the inhibitors with glycosyl. Here, we propose that the glycosyl functional groups in inhibitors may be instrumental in inhibiting SARS-CoV-2. ties in small molecules interact with hydrophobic residues of S_RBD protein. Among all the potential inhibitors, Phillyrin and Liquiritin have relatively lower docking energy in domain 1. The three-dimensional binding modes between the Phillyrin/Liquiritin and S_RBD are shown in Figure 7D,E. From Figure S7, we can see that hydrophobic residues are important in the interactions between S_RBD and the inhibitors with glycosyl. Here, we propose that the glycosyl functional groups in inhibitors may be instrumental in inhibiting SARS-CoV-2.  To explore the potential inhibitors for the genetic variants of SARS-CoV-2, we calculated the lowest free binding energy of Lincomycin, Phillyrin, and Liquiritin with the spike RBD of the genetic variants of SARS-CoV-2 using AutoDock, as shown in Table S3. Interestingly, we found that the value of the binding energy between Phillyrin and m_S_RBD (mutant site: E484K) is −20.585 KJ/mol, and it is −22.635 KJ/mol between Liquiritin and m_S_RBD (mutant sites: L452R, T478K). These are lower than the value between Lincomycin and S_RBD. We also plotted the binding models of the inhibitors and the m_S_RBD in the genetic variants of SARS-CoV-2, as shown in Figure 7F,G. The types of interactions between the inhibitors and m_S_RBD are the hydrophobic interactions, π-sigma interaction and pi-sulfur interaction. The glycosyl group plays a crucial role in the interactions of the protein and inhibitor.
In short, the small molecules with the glycosyl group show the potential inhibitions for S_RBD and m_S_RBD.

Conclusions
In this work, we investigated the conformational characteristics of S_RBD in the six simulation systems. Based on the 12 MD simulations, we concluded that two types of binding sites are explored in the four systems with inhibitors, and small molecules can reduce the flexibility of the S_RBD protein. Hydrophobic and π-π stacking interactions played vital important roles in the interplay between S_RBD and the small molecules. More interestingly, our data indicated that Lincomycin containing glycosyl group shows a stronger role in stabilizing S_RBD protein than the other inhibitors and the human ACE2 protein.
In short, our results explained the interactional mechanism of S_RBD and inhibitors at the atomic level, revealed two binding sites of S_RBD, and discovered the importance of the glycosyl group of inhibitors in inhibiting SARS-CoV-2 and the genetic variants of SARS-CoV-2. In silico-derived mechanistic insights into the interactions of S_RBD and inhibitors are provided in this study, which may provide new clues for fighting SARS-CoV-2 infection.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/pr9101749/s1, Figure S1: The initial structures of complex. Figure S2: The initial structure in the S_RBD with ACE2 system and the S_RBD system by cartoon representation. Figure S3: The side and top views of two binding regions. Figure S4: The alignment of the initial conformation and the representative snapshot at 400 ns in each simulation with the cartoon representation. Figure S5: The residue-based contact probability between S_RBD and part I/part II of the inhibitor Lincomycin. Figure S6: The chemical structures of potential inhibitors with the glycosyl group. Figure S7: The binding models between the inhibitors and S_RBD in binding domain 1 with Discovery Studio method. Table S1: The detail of each simulation system. Table S2: The lowest free binding energy of potential effective inhibitors and the in different domains using AutoDock. Table S3: The lowest free binding energy of Lincomycin, Phillyrin and Liquiritin with the mutated Spike RBD of the genetic variants of SARS-CoV-2 using Autodock. Table S4: The respective RMSF value of residues belonging to domain 1 and domain 2. Table S5: Free energy between S_RBD and the small molecules Cefuroxime, Erythromycin, Lincomycin and Ofloxacin. Table S6: The residues of binding energy less than −1 KJ/mol in four systems. Table S7: The specific value of residue-based contact probability more than 0.2 between S_RBD and part I/part II of Lincomycin in domain 1.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.