Molecular Docking and Dynamics Studies to Explore Effective Inhibitory Peptides Against the Spike Receptor Binding Domain of SARS-CoV-2

The spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has become a pandemic due to the high transmission and mortality rate of this virus. The world health and economic sectors have been severely affected by this deadly virus, exacerbated by the lack of sufficient efficient vaccines. The design of effective drug candidates and their rapid development is necessary to combat this virus. In this study, we selected 23 antimicrobial peptides from the literature and predicted their structure using PEP-FOLD 3.5. In addition, we docked them to the SARS-CoV-2 spike protein receptor-binding domain (RBD) to study their capability to inhibit the RBD, which plays a significant role in virus binding, fusion and entry into the host cell. We used several docking programs including HDOCK, HPEPDOCK, ClusPro, and HawkDock to calculate the binding energy of the protein-peptide complexes. We identified four peptides with high binding free energy and docking scores. The docking results were further verified by molecular dynamics (MD) simulations to characterize the protein-peptide complexes in terms of their root-mean-square fluctuation (RMSF), root-mean-square deviation (RMSD), radius of gyration (Rg), solvent-accessible surface area (SASA), and hydrogen bond formation. Allergenicity and toxicity predictions suggested that the peptides we identified were non-allergenic and non-toxic. This study suggests that these four antimicrobial peptides could inhibit the RBD of SARS-CoV-2. Future in vitro and in vivo studies are necessary to confirm this.

The S proteins of coronavirus consist of spike monomers with two subunits, S1 and S2 (Gui et al., 2017;Yuan et al., 2017;Kirchdoerfer et al., 2018;Song et al., 2018;Lan et al., 2020). The S1 subunit contains the receptor-binding domain (RBD) and N-terminal domain (NTD) which are responsible for virus binding and entry. The RBD is located in the middle part of the S1 subunit and is used as an antigen to raise antibodies that interrupt virus-host binding (Xiao et al., 2003;Babcock et al., 2004;He et al., 2004;Wong et al., 2004;Lan et al., 2020;Dejnirattisai et al., 2021). The S2 domain has a proposed fusion peptide and two heptad repeats (HR1 and HR2) that facilitate cell membrane fusion between viral and target cells following proteolytic activation (Wild et al., 1994;He et al., 2004;Shang et al., 2020;Kim et al., 2021). The RBD-containing S1 domain also contains the SD1 and SD2 subdomains at the C-terminus. Although both the NTD and the RBD are immunogenic, the RBD is contains the interaction surface for ACE2. Due to the fact that the receptor-binding site (RBS) is incompletely driven into the down state, the RBD solely engages with the up state of ACE2 (Lan et al., 2020;Premkumar et al., 2020;Dejnirattisai et al., 2021;Yuan et al., 2021).
Usually, many years of research are required before vaccines enter clinical trial. However, in a record period, scientists and researchers have made great efforts to develop secure, efficient, and active SARS-CoV-2 vaccines. Currently, 13 vaccines have been approved for early or limited use, and 8 vaccines have been approved for complete use (https://www.nytimes.com/ interactive/2020/science/coronavirus-vaccine-tracker).
In this study, the selected peptides were docked to the RBD of SARS-CoV-2, leading to the identification of four peptides with high binding free energy. These peptide-RBD complexes were subsequently subjected to molecular dynamics study. Structural attributes and conformations of the docked complexes were obtained from the MD simulations, and suggested stiff and inflexible interactions between the RBD active site and the hit peptides. In comparison to earlier studies, we utilized multiple docking programs in combination to identify four peptides with high binding affinity to the active site of the RBD. Several previous studies have suggested antiviral effects of small molecules and peptides against SARS-CoV-2 through binding to the RBD. These studies identified molecules that, although they were predicted to bind to the RBD, did not interact directly with the RBD active site (Rathod et al., 2020;Padhi et al., 2021;Priya et al., 2021). In contrast, all of the peptides we identified were predicted to bind directly to the RBD active site. Although a few previous studies identified peptides that formed a single nonbonded interaction with the RBD active site, the predicted binding energies were lower than ours using the same docking software (Chowdhury et al., 2020;Hossain et al., 2021) and the complexes were less stable in MD simulations.

Peptide Screening and Preparation
In this study, we started from 27 peptide molecules that were previously identified in the venom of the wild bee Hylaeus signatus and which were screened for antimicrobial activity (Nešuta et al., 2016). Three peptides were excluded as they contained D-amino acids. Additionally, one further peptide was excluded as its amino acid sequence was incompletely characterized. The PEP-FOLD 3.5 webserver was used to predict the peptide structures from the amino acid sequences of the peptides (Lamiable et al., 2016). This webserver uses a Hidden Markov Model suboptimal sampling algorithm to predict the peptide structures. The resulting peptide structures were used as the starting point for 20 ns molecular dynamics simulations, and the root mean square deviations of the alpha carbon atoms were calculated. The final frames of these molecular dynamics simulations were used for further studies.

Protein Preparation
The three-dimensional structure of the spike receptor-binding domain (RBD) of SARS-CoV-2 at 2.43 Å resolution (PDB ID: 6M0J) was retrieved from the Protein Data Bank. The protein structure was prepared by removing heteroatoms and water molecules using Discovery Studio (Discovery Studio, 2009). Additionally, energy minimization of the protein structure was performed using the AMBER14 (Case et al., 2014) force field in YASARA software (Krieger et al., 2013). Molecular docking and dynamics studies used this energy minimized protein structure.

Molecular Docking
The peptides and RBD protein were uploaded as ligand and receptor molecules respectively to the HDOCK, HPEPDOCK, and ClusPro web servers. HDOCK uses a combined templatebased and template-free algorithm in an automatic manner (Yan et al., 2017), while HPEPDOCK uses a hierarchical algorithm (Zhou et al., 2018). After docking the peptides to the protein using HDOCK and HPEPDOCK, ClusPro was used to calculate binding energies (Comeau et al., 2004). The top ten peptides with the highest docking scores were selected for further evaluation. These ten peptides were docked to the RBD using the HawkDock web server. For further analysis, the four highest-scoring peptides were chosen based on their binding free energy and docking scores from HawkDock. On the HawkDock server, the HawkRank scoring system, the ATTRACT docking algorithm improved in groups, and MM/GBSA free energy decomposition analysis are implemented on a multipurpose platform (Weng et al., 2019). PyMOL and Discovery Studio (Discovery Studio, 2009) were used for structural analysis of the top four protein-peptide complexes.

Molecular Dynamics Simulation
Molecular dynamics simulations was performed in YASARA dynamics (Land and Humble, 2018) using the AMBER14 force field (Wang et al., 2004). The docked peptide-protein complexes were initially cleaned, optimized and the hydrogen bond network was oriented. A cubic simulation cell was created with periodic boundary conditions and the TIP3P water model was used (Harrach and Drossel, 2014). The simulation cell was extended by 20 Å in each direction beyond the protein-peptide complexes. The physical conditions of the simulation cells were set at 298 K, pH 7.4, and 0.9% NaCl (Krieger and Vriend, 2015). The initial energy minimization of the simulation cells were conducted by steepest gradient approaches with simulated annealing methods (5,000 cycles). The time step of the simulation was 2.0 fs. Long-range electrostatic interactions were calculated by the Particle Mesh Ewald (PME) method with a cut-off radius of 8.0 Å (Essmann et al., 1995;Krieger et al., 2006;Harvey and De Fabritiis, 2009). The simulation trajectories were saved every 100 ps. The simulations were run for 100 ns at constant pressure and temperature, using a Berendsen thermostat. The simulation trajectories were used to calculate root mean square deviations, root mean square fluctuations, solvent accessible surface areas, radii of gyration, and hydrogen bonds (Adji et al., 2021;Dutta et al., 2021;Obaidullah et al., 2021).
The per residue energy contribution of the peptide-protein structures was calculated using the pyDockEneRes webserver. The first, last, and average structure was extracted from simulations trajectories and utilized as input entry, and the average higher energy from the hotspot residues were tabulated. This tool can be utilized for the identification of the hotspot residues (Romero-Durana et al., 2020).

Allergenicity and Toxicity Prediction
AllerTOP (Dimitrov et al., 2014a) and AllergenFP (Dimitrov et al., 2014b) webservers were used to predict the allergenicity of the peptides. AllergenFP implements five E-descriptor-based fingerprinting, whilst AllerTOP uses both k-nearest neighbor (kNN) and amino acid E-descriptors to predict the allergenicity of peptides. The ToxinPred (Gupta et al., 2013) webserver was utilized to predict the toxicity of the peptides.

Molecular Docking
Venom extracted from hymenopteran insects, including the solitary bee Hylaeus signatus, is a prominent source of antimicrobial peptides (AMPs). Several α-helical amphipathic AMPs, referred to as HYL, have been identified in the solitary bee, with antimicrobial activity against distinct strains of pathogenic bacteria and fungi as well as the ability to lyse cancer cells Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 8 | Article 791642 (Slaninová et al., 2011;Slaninová et al., 2012;Nešuta et al., 2016). Additionally, HYL had low hemolytic activity, suggesting that it may be safe for use in humans (Nešuta et al., 2016). Notably, the earlier-discovered antimicrobial peptides, including cecropin-A from Hyalophora cecropia, exhibited antiviral activity (Both, 1980). Cecropin-A was shown to have antiviral activity against HSV-1, HIV-1, and JV, whereas melittin from A. mellifera had antiviral activity against influenza A viruses, HSV-1, HIV-1, JV, TMV, RSV, VSV, enterovirus, and coxsackievirus (Albiol Matanic and Castilla, 2004). In addition,  Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 8 | Article 791642 4 alloferon 1 and alloferon 2 from C. vicina were effective against influenza viruses (Chernysh et al., 2002), TnGlv1 and TnGlv2 from Trichoplusia ni were effective against AcMNPV (Moreno-Habel et al., 2012), attC and dptB from Drosophila melanogaster were effective against SINV (Huang et al., 2013), and a myristoylated peptide from H. virescens was effective against HIV-1 and HSV-1 (Ourth, 2004;Feng et al., 2020). Therefore, we hypothesized that the HYL antimicrobial peptides from the solitary bee could have antiviral activity and consequently devised an in silico study to test this hypothesis. We found that several HYL peptides bind the RBD of SARS-CoV-2 in silico with favorable binding energy and stable conformations through 100 ns of molecular dynamics simulation.
Molecular docking can be utilized to rapidly assess the binding affinities and modes between a target substrate, such as a protein, and diverse ligands, including small peptides, to assist antiviral drug design (Willett and Glen, 1995). Virtual screening can identify the preferred binding orientation, optimal conformation, and binding sites of protein and peptide molecules (Alonso et al., 2006;Liu et al., 2018). Our computational approach confers improved structural precision and rapid prediction of both the existence and strength of binding through the use of multiple diverse docking algorithms aimed at achieving accuracy (Bartuzi et al., 2017).
The predicted peptide structures from PEP-FOLD 3.5 are shown in Supplementary Figure 1. The peptide structures were optimized by molecular dynamics simulations in which the root mean square deviations of the peptides were found to be small (Supplementary Figure 2). After predicting the peptide structures, all 23 peptides were docked to the RBD of SARS-CoV-2 by using the HDOCK, ClusPro, and HPEPDOCK webservers (Supplementary Table 1). Binding energies and/or docking scores were recorded for the ten highest-scoring peptides (Supplementary Table 2). Moreover, the binding free energies and docking scores from HawkDock (Supplementary Table 3) were used to select the four highest-scoring peptides (Figure 1). The highest docking score was −2,848.16 kcal/mol for peptide P20, followed by −2,815.13, −2,735.49, and −2,660.34 kcal/mol for P18, P21, and P2, respectively. The maximum binding free energy was −40.38 kcal/mol for the P20 peptide, followed by −37.2, −31.17, and −28.15 for P18, P21, and P2, respectively ( Table 1).
The interactions between the RBD protein and the remaining six of the top ten peptides are shown in supplementary figure 3. The P3:RBD complex formed two hydrogen bonds at ALA344 and ASN343, five hydrophobic interactions at LEU441, LEU335, PHE374, VAL362, and VAL367, and one unfavorable bond at ARG509. The P6:RBD complex had five hydrophobic interactions at LEU441, PHE374, LEU368, VAL367, and LEU335, and six hydrogen bonds at ASN440, ARG509, TRP436, ASN343, PHE342, and GLY339 position. The P16:RBD complex had four hydrophobic interactions at CYS391, CYS525, ALA522, and VAL362, two hydrogen bonds at ASP389 and THR523, one salt bridge at ASP389 and one unfavorable bond at CYS361. The P17: RBD complex had five hydrophobic interactions at LEU335, VAL367, LEU368, PHE342, and TRP436, and three hydrogen bonds at SER371, ASN343, and ASN440. The P22:RBD complex formed five hydrogen bonds at ASN440, ASN343, SER371, GLU340, and ASP364, three hydrophobic interactions

Molecular Dynamics
Molecular dynamics simulations were conducted to explore the structural stability of the docked peptide-protein complexes. The root mean square deviations (RMSD) from the peptide-protein systems were calculated from the simulation trajectories. Figure 3A indicates that the peptide-protein complexes involving the P2, P18, P20, and P21 peptides had an initial upwards RMSD trend, suggesting flexibility of the complexes. The peptide-protein systems subsequently stabilized after 30 ns and maintained their integrity for the remaining 70 ns of the simulations. The P21 complex had a comparatively higher RMSD than the other three peptide complexes, which indicates the more flexible nature of this complex. However, all four peptideprotein complexes had RMSD lower than 2.5 Å, which indicates the stable nature of the complexes over the whole simulations. The solvent-accessible surface area (SASA) of the complexes was also explored. Figure 3B indicates that the P20 and P21 complexes increased in SASA upon binding with the protein target, and after 30 ns the P21 complex decreased in SASA,  maintaining a similar SASA trend to the rest of the complexes. The P20 complex stabilized after 50 ns and maintained its integrity. The radius of gyration (Rg) of the complexes was also explored to understand the compactness of the complexes. High Rg indicates an extended nature whereas a low Rg indicates a more compact structure. Figure 3C indicates that P20 had slightly higher Rg than other complexes which indicates that this complex is less compact. The other complexes had smaller deviations, which indicate the rigid nature of these complexes. Hydrogen bonding can play a crucial role in determining the stability of protein complexes. Figure 3D indicates that all complexes had a stable hydrogen bonding profile over the 100 ns of simulation. The root mean square fluctuations (RMSF) of the complexes were explored to understand the flexibility on a residue-byresidue basis. Figure 4 indicates that almost every residue of the complexes had an RMSF lower than 1.4 Å, indicating the stable nature of the protein-peptide systems.
The active sites of the RBD protein include amino acid residues 340, 374, 375, 378, 403, 420, 477, 478, 499, 543, and 546 (Lan et al., 2020). At the start of the simulation, the P2 peptide bound to the active site of RBD via a hydrophobic interaction with the conserved residue PHE374 (bond distance 3.31360), whereas the P18 peptide interacted with the RBD active site at both PHE374 (Bond distance 3.41730) and GLU340 (bond distance 2.70389) with a hydrophobic interaction and a hydrogen bond respectively. The P20 peptide bound to the active site of RBD at PHE374 (bond distance 3.44294) via a hydrophobic interaction, and the P21 peptide interacted with the RBD active site at SER375 (bond distance 2.49462) and LYS378 (bond distance 3.83247) position with a hydrogen bond and a hydrophobic interaction respectively (Table 2; Figure 2).
After 100 ns of simulation time, the P18 peptide was bound to the active site of RBD at two conserved residues, PHE374 (bond distance 3.44205) and GLU340 (bond distance 2.73681), via a hydrophobic interaction and a hydrogen bond respectively. The P21 peptide interacted with the active site of RBD via LYS378 (bond distance 2.81702) and PHE374 (bond distance 2.87497) via two hydrogen bonds (Supplementary table 5, Supplementary figure 4). These interactions of the peptides with the active site of the RBD indicate that these peptides may inhibit the RBD with high binding affinity.
Previously suggested multifunctional peptides, including ALPEEVIQHTFNLKSQ and DIENLIKSQ from Bacillusfermented soybean, failed to bind to the active site of the RBD. However, our top four peptides bound to the active site. In addition, P18 and P20 had higher binding energy than DIENLIKSQ when employing the same HawkDock server (Padhi et al., 2021). The computationally designed peptide APASMFLGKGDHEILM made no interactions with the active site when docked using the same HPEPDOCK server that we used for our top four peptides (Priya et al., 2021). Another peptide modeling and screening study suggested that AVP0671 can bind to the RBD, albeit not at the active site, and the HDOCK scores were not reported. Meanwhile, our top four suggested peptides not only exhibited active site binding but also had a more favorable HDOCK docking score (Rathod et al., 2020). Antiviral peptides S2P25 and S2P26 were predicted to form one bond at the active site, but had lower docking interaction energies than our P21 and P18 peptides when comparing the Cluspro scores (Chowdhury et al., 2020). Additionally, AVP1795, identified from the computational screening of 645 antiviral peptides, exhibited only one interaction (ARG403) at the active site, which was occupied for ≥90% of the time in MD simulations, whereas our peptides were stably bound over the whole simulation period (Hossain et al., 2021).
Multiple methodology developed to understand the hotspot residues in protein-protein interactions which are based on snapshot from MD simulation sampling (Gohlke et al., 2003). In addition to pyDockEneRes, mm_pbsa.pl tools successfully implemented for hotspot resdiues identifications from cytokines and receptor interface  as well as RBD from SARS-CoV-2 and antibodies (Yang et al., 2021). Moreover, the pyDockEneRes tools enables the hotspot residues identifications from the protein-protein interactions which is key to understand the biological process at molecular level (Romero-Durana et al., 2020). In our study, the P2 peptide and RBD complex, four important residues contribute to overall energy contributions: Leu335, Glu340, Val367, and Phe374; in the P18-RBD complex, hotspot residues include Asp164, Ser373, Ser375, and Trp380 (Table 3). Leu335, Asn343, Asp364, Val367, Ser373, Phe374, and Glu484 were major energy contributing residues in the P20-RBD complex, whereas Leu335, Phe342, Asn343, Asp364, Val367, Ser373, Phe374, and Glu484 were key energy contributing residues in the P22-RBD complex.

Allergenicity and Toxicity Prediction
An allergenic antigen can activate Th2 cells resulting in stimulation of B cells to generate immunoglobulin E (IgE) that binds to FcεRI and activates eosinophils leading to inflammation and tissue shrinkage (Dimitrov et al., 2013). The online AllerTOP tool, which assesses allergens using E-descriptors affined with amino acid attributes, was utilized to predict the allergenicity of our peptides (Dimitrov et al., 2014a). Three distinct web-based approaches were utilized to assess toxicity ( Table 4). According to the AllergenFP v1.0 webserver, P2 is a probable allergen, whereas P18, P20, and P21 are probably non-allergenic. The same result was obtained using the AllerTOP v2.0 webserver. In addition, all four peptides (P2, P18, P20, and P21) were categorized as non-toxic by the ToxinPred web server.

CONCLUSION
To develop novel therapeutics against SARS-CoV-2, targeting the spike protein RBD by designing antiviral peptides could identify promising leads. In this study, 23 peptides were docked to the RBD of SARS-CoV-2, leading to the identification of four peptides with high binding affinity. Molecular dynamics studies demonstrated that docked peptides were not ductile in nature, but were instead rigidly bound. Moreover, allergenicity and toxicity profiling of the peptides suggest that they have no allergenic or toxic properties. Finally, our study may facilitate the development of efficient drugs against SARS-CoV-2 through further in vitro studies.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
SM: Conceptualization, investigation, data curation, validation, software, formal analysis, writing, editing, visualization, and supervision; SB: investigation, data curation, validation, software, formal analysis, writing, editing, and visualization; MM: formal analysis, writing, and visualization; SA: formal analysis and writing; MH: formal analysis and writing; SS: formal analysis and writing, MS: editing and supervision; GM-H: editing and writing, MA: editing and writing; AO: editing and writing, GB: editing and writing. All authors have read and agreed to the current version of the manuscript.

FUNDING
The authors would like to extend their sincere appreciation to Taif University Researchers supporting Project number (TURSP 2020/309) Taif University, Taif, Saudi Arabia.