Molecular Docking and Dynamics Simulation of a Screening Library from Life Chemicals Database for Potential Protein-Protein Interactions (PPIs) Inhibitors against SARS-CoV-2 Spike Protein

The ongoing pandemic of coronavirus 2 represents a major challenge for global public health authorities. Coronavirus disease 2019 can be fatal especially in elderly people and those with comorbidities. Currently, several vaccines against coronavirus 2 are under application in multiple countries with emergency use authorization. In the same time, many vaccine candidates are under development and assessment. It is worth noting that the design of some of these vaccines depends on the expression of receptor binding domain for viral spike protein to induce host immunity. As such, blocking the spike protein interface with antibodies, peptides or small molecular compounds Original Research Article Odhar et al.; JPRI, 33(20A): 74-84, 2021; Article no.JPRI.66801 75 can impede the ability of coronavirus 2 to invade host cells by intervention with interactions between viral spike protein and cellular angiotensin converting enzyme 2. In this virtual screening study, we have used predictive webservers, molecular docking and dynamics simulation to evaluate the ability of 3000 compounds to interact with interface residues of spike protein receptor binding domain. This library of chemicals was focused by Life Chemicals as potential proteinprotein interactions inhibitor. Here, we report that hit compound 7, with IUPAC name of 3‐cyclohexyl‐N‐(4‐{[(1R,9R) ‐6‐oxo‐7,11‐ diazatricyclo [7.3.1.0 2,7 ] trideca‐2,4‐dien‐11‐yl] sulfonyl} phenyl) propenamide, may have the capacity to interact with interface of receptor binding domain for viral spike protein and thereby reduce cellular entry of the virus. However, in vitro and in vivo assessments may be required to validate these virtual findings.


INTRODUCTION
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative pathogen responsible for the ongoing pandemic of coronavirus disease 2019 (COVID-19) [1]. SARS-CoV-2 is a novel RNA virus and has about 79% of genomic sequence identity with SARS-CoV, a previously identified beta-coronavirus [2]. SARS-CoV-2 can be transmitted mainly through exposure to respiratory droplets generated by infected persons. Airborne transmission by aerosols is also possible during aerosol generating medical procedures like intubation [3]. Upon exposure to SARS-CoV-2, the mean incubation period can range between 5.6 and 6.7 days after which symptoms of COVID-19 may develop [4]. Common clinical features of COVID-19 include fever, dry cough, shortness of breath, fatigue, myalgia, nausea, vomiting and diarrhea. The COVID-19 infection fatality rate (IFR) is estimated to be 0.68% in general with 95% confidence interval between 0.53% and 0.82%. COVID-19 fatality is usually the result of complications like pneumonia, acute respiratory distress syndrome (ARDS), cardiac injury, liver injury, kidney injury and coagulopathy [5,6]. Therapeutic options currently available to control COVID-19 complications involve oxygen supply, management of sepsis, use of corticosteroids, heparin and antiviral agents like Remdesivir [5,[7][8][9]. According to World health organization (WHO), as of 18 February 2021, seven anti-COVID-19 vaccines are under emergency use authorization in multiple countries and more than 60 vaccine candidates are under clinical development [10]. Different platforms have been utilized to produce these COVID-19 vaccines like inactivated virus, protein subunit, recombinant viral-vector and nucleic acid (DNA and mRNA). Some of these vaccines have focused on the expression of full-length SARS-CoV-2 spike protein inside human body to elicit immune response, other vaccines used only the receptor binding domain (RBD) of spike protein to induce immunity [11]. SARS-CoV-2 can infect susceptible host through interaction with angiotensin converting enzyme 2 (ACE2) on the surface of target cells. The receptor binding domain (RBD) of SARS-CoV-2 is located within S1 subunit of the viral spike protein [12]. Virtual alanine scanning approach have been used to determine the interacting residues at the interface of SARS-CoV-2 receptor binding domain and angiotensin converting enzyme 2. Identification of these interface residues is considered important to construct a protein-protein interactions (PPIs) inhibitor capable of blocking cellular entry of SARS-CoV-2 [13]. A three-dimensional illustration for the complex between RBD of SARS-CoV-2 spike protein and ACE2 can be seen in Fig. 1 (A). A close view for interaction interface of crystal complex can be seen in Fig. 1 (B). The sequence of amino acids for RBD of SARS-CoV-2 spike glycoprotein is presented in Fig. 1 (C), RBD residues involved in interaction with ACE2 are surrounded by orange line. In this virtual screening study, we have used a proteinprotein interactions (PPIs) focused library of compounds from Life Chemicals database. Our aim was to identify potential inhibitor of SARS-CoV-2 spike protein binding to ACE2 thereby blocking viral entry to target cells. The chemical compounds in this library were curated based on rule of four (RO4) principle. According to rule of four, a compound may have the potential to inhibit protein-protein interactions (PPIs) if it possesses the following four chemical criteria: a molecular weight (M.W.) ≥ 400 g/ mol, a calculated logarithm of partition coefficient (cLog P) ≥ 4, the number of rings ≥ 4 and the number of hydrogen bond acceptor (HBA) ≥ 4 [14,15].

Preparation of Receptor Binding Domain Crystal and Database
The crystal complex of receptor binding domain RBD of SARS-CoV-2 spike protein and angiotensin converting enzyme 2 ACE2 was downloaded from protein data bank website, the identification code for the downloaded crystal PDB is 6LZG [16,17]. We have used UCSF chimera version 1.15 to prepare the crystal file for docking study and dynamics simulation By using UCSF chimera software, we have removed ACE2 peptidase domain (chain A) from 6LZG crystal, we have also removed any bound ligands, water molecules and ions. The residues sequence of chain B (receptor binding domain) was visualized as one letter format by using UCSF chimera. We have used 3000 compounds downloaded from Life Chemicals website as a screening library, this library of chemicals was focused by Life Chemicals as potential inhibitors for protein-protein interactions (PPIs) by using rule of four (RO4) filtration criteria. Based on RO4 principle, any compound may have the capacity to interfere with interactions between two proteins if it has the following criteria: a molecular weight ≥ 400 g/ mol, a logarithm of dimensional view for 6LZG crystal complex between receptor binding 2 spike protein and its receptor angiotensin converting enzyme 2 A close view for interaction interface between SARS-CoV-2 spike protein RBD and 2, the amino acid residues of spike protein RBD are shown as sticks. (C) The amino acids 2 spike protein RBD is presented in one letter format, RBD residues that are involved in interaction with ACE2 are encircled by orange line

Receptor Binding and Chemicals
The crystal complex of receptor binding domain 2 spike protein and converting enzyme 2 ACE2 was downloaded from protein data bank website, the identification code for the downloaded crystal . We have used UCSF chimera version 1.15 to prepare the crystal file for docking study and dynamics simulation [18]. By using UCSF chimera software, we have removed ACE2 peptidase domain (chain A) from 6LZG crystal, we have also removed any bound ligands, water molecules and ions. The residues sequence of chain B (receptor binding domain) r format by using UCSF chimera. We have used 3000 compounds downloaded from Life Chemicals website as a screening library, this library of chemicals was focused by Life Chemicals as potential inhibitors protein interactions (PPIs) by using e of four (RO4) filtration criteria. Based on RO4 principle, any compound may have the capacity to interfere with interactions between two proteins if it has the following criteria: a ≥ 400 g/ mol, a logarithm of partition coefficient ≥ 4, the number of hydrogen bond acceptor ≥ 4 and the number of rings ≥ 4. Additionally, the compounds in this focused library have been enriched for sp3 hybridization of carbon atoms to ensure high level of three dimensional diversity and hence sufficient dru likeness score [14,15].

Structure-based Virtual Screening
We have used MCULE online drug discovery website to screen the downloaded chemicals library against the prepared crystal of SARS CoV-2 spike protein RBD (chain B of 6LZG crystal) [19]. Here, we used a virtual screening workflow steps similar to those applied in our previously published researches summary, we have used the default screening steps with the addition of REOS (Rapid Elimination of Swill) filter to reduce possibility of frequent and non-selective hits. dimensional view for 6LZG crystal complex between receptor binding 2 spike protein and its receptor angiotensin converting enzyme 2 2 spike protein RBD and The amino acids 2 spike protein RBD is presented in one letter format, RBD residues E2 are encircled by orange line the number of hydrogen ≥ 4 and the number of rings ≥ 4. Additionally, the compounds in this focused library have been enriched for sp3 hybridization of carbon atoms to ensure high level of threedimensional diversity and hence sufficient drug-

Virtual Screening
We have used MCULE online drug discovery website to screen the downloaded chemicals library against the prepared crystal of SARS-2 spike protein RBD (chain B of 6LZG we used a virtual screening workflow steps similar to those applied in our previously published researches [20,21]. In summary, we have used the default screening steps with the addition of REOS (Rapid Elimination of Swill) filter to reduce possibility of selective hits. The MCULE website employs AutoDock Vina tool to perform structure based virtual screening [22]. Also, this online drug discovery platform uses AutoDock tools to add polar hydrogen atoms and Gasteiger [23]. We have used a binding site area of (X: 30, Y:30, Z:30) Angstrom while the coordinates were (X: -33, Y: 26, Z: 7.5). The final hits were ranked according to their least binding energy to the receptor binding domain (RBD) of SARS-CoV-2 spike protein crystal (chain B). According to the minimum binding energy rank, we have selected the top ten hits for further virtual characterization and visualization. For each compound of these top ten hits, we have saved the complex of ligand and protein with least binding energy pose as PDB file. The generated PDB file was then visualized as two and three-dimensional illustrations by using PyMOL version 2.4.1 and Discovery Studio Visualizer version 21.1.0.20298 respectively [24,25]. The two-dimensional chemical structures of these top ten hits were drawn by using MarvinSketch version 20.1 [26].

Prediction of Chemical, Pharmacokinetics and Mutagenic Characteristics of Hit Compounds
MCULE online platform provides the opportunity to predict various chemical features of the generated hit compounds. Drug-likeness score for the top ten hits were predicted by using Molsoft webserver [27]. Various pharmacokinetics and mutagenic characteristics for the top hits were predicted and calculated by using SwissADME and pkCSM webservers [28,29]. These webservers use predictive regression and molecular similarity to analyze molecules under investigation [30,31].

Molecular Dynamics (MD) Study
We have employed YASARA Dynamics version 20.12.24 to perform molecular dynamics (MD) simulation of the ligand-protein complex with the least binding energy pose [32]. The steps for molecular dynamics simulation study are similar to what we had followed in our previously published researches [21,33]. In summary, the MD protocol involves hydrogen bonds network optimization and a pKa prediction to fine-tune the protonation of residues at pH value equal to 7.4 [34]. NaCl was added with a concentration of 0.9%, an excess of either sodium or chloride ions were used to neutralize ligand-protein complex. To eliminate any possible clashes, steepest descent and simulated annealing minimizations were applied. The simulation period employed was 10 nanoseconds by using AMBER14 force field for solute, TIP3P for water, AM1BCC and GAFF2 for ligand [35][36][37]. The cutoff value for van der Waals (vdW) forces was 8 Angstrom, the default parameters were employed by AMBER [38]. By applying Particle Mesh Ewald algorithm, no cutoff value was used for electrostatic forces [39]. The motions equations were used as a multiple timestep of 1.25 and 2.5 femtoseconds for bonded and non-bonded interactions respectively at a pressure of 1 atm and a temperature of 298K [40]. After assessment of root-mean-square deviation (RMSD) for the solute as a function of simulation period, the first 10 nanoseconds were regarded as the equilibration time and precluded from further analysis.

RESULTS AND DISCUSSION
A prediction of different chemical characteristics for the top ten hit compounds that were computationally screened against RBD crystal of SARS-CoV-2 spike protein can be seen in Table  1. In this table, the hit compounds were ranked according to their minimum binding energy to RBD crystal. According to Table 1, all these protein-protein interactions (PPIs) inhibitor candidates have a molecular weight greater than 400 g/ mol and a number of hydrogen bond accepters greater than 4. Regarding the prediction of partition coefficient logarithm (Log P), all the reported hits have Log P value greater than 4 with the exception of compounds 4, 9 and 10. The chemical structures for these top ten hits can be seen in Fig. 2 where all these compounds have more than 4 rings within their structures. Based on chemical data presented in Table 1 and Fig. 2, we can assess the adherence of these hit compounds to rule of four criteria.
In Table 2, various pharmacokinetics properties together with mutagenic potential are predicted and listed for the top ten hit compounds. These hits were ordered in Table 2 according to their predicted least binding energy as presented in this table. It is obvious that compounds 4, 9 and 10 have one violation for rule of four (RO4) criteria as seen in Table 2. The logarithm of partition coefficient (Log P) for these compounds is less than 4 as reported in Table 1. All the reported hits have a high predicted drug-likeness score with the exception of compounds 2, 3 and 8. Both compounds 1 and 6 were precluded from further analysis due to their anticipated poor water solubility. Additionally, compounds 4 and 10 may have a mutagenic capacity as predicted by AMES toxicity and therefore are excluded from any further examination. According to Table  2, only compounds 5 and 7 pass all the criteria for rule of four (RO4), have high drug-likeness score, possess moderate water solubility with no mutagenic capacity as indicated by predictive webservers. As such, only compounds 5 and 7 were considered for further evaluation of docking images and dynamics simulation.    Careful examination of two and three dimensional illustrations for compound 5 and 7 docking against RBD crystal with least binding energy pose, as presented in Fig. 3 both compounds may be involved in different kinds of chemical interactions with interface residues of RBD crystal. Of interest is the possible ability of compound 5 to form conventional hydrogen bond with Glycine 496 as seen in Fig. 3 (A). On the other hand, compound 7 may have the capacity to form conventional hydrogen bond with Tyrosine 449 as noticed in Careful examination of two and threeimensional illustrations for compound 5 and 7 docking against RBD crystal with least binding Fig. 3, reveal that both compounds may be involved in different kinds of chemical interactions with interface Of interest is the possible ability of compound 5 to form conventional hydrogen bond with Glycine 496 as . On the other hand, compound 7 may have the capacity to form conventional hydrogen bond with Tyrosine 449 as noticed in Fig. 3 (B). Both Tyrosine 449 and Glycine 496 are considered interface residues of RBD crystal as can be seen in Fig. 1.
Finally, the results for compounds 5 and 7 molecular dynamics (MD) simulation were reported in both Fig. 4 and superposing the receptor on its reference structure, ligand proximity to RBD interface residues can be recorded as a function of simulation period as seen in Finally, the results for compounds 5 and 7 molecular dynamics (MD) simulation were and Fig. 5. By superposing the receptor on its reference structure, ligand proximity to RBD interface residues can be recorded as a function of Fig. 4. It is wellevident that compound 5 has moved away from initial binding site by the end of simulation period as noticed in Fig. 4 (A), this may refer to weak interactions between compound 5 and RBD interface residues. On the other hand, shows that compound 7 can maintain a more constant distance from RBD interface as reported by root-mean-square deviation (RMSD) of ligand movement throughout 10 nanoseconds. This may indicate that compound 7 has stronger interactions with RBD crystal interface than does compound 5. Conformational changes RMSD of d of simulation period , this may refer to weak interactions between compound 5 and RBD interface residues. On the other hand, Fig. 4 (B) shows that compound 7 can maintain a more constant distance from RBD interface as square deviation (RMSD) of ligand movement throughout 10 nanoseconds. This may indicate that compound 7 has stronger interactions with RBD crystal interface than does compound 5. Conformational changes RMSD of compounds 5 and 7 can be seen in and (B) respectively. By superposing the ligand on its reference structure throughout simulation period, we can observe that the conformational changes RMSD for compounds 5 and 7 in are consistent with ligand movement RMSD for these two compounds in Fig. 4. Based on these results, it is obvious that compound 7 may have more capacity than do compound 5 to bind interface residues of receptor binding domain for SARS-CoV-2 spike protein.

square deviation (RMSD) of ligand movement as a function of simulation interval. Plot (A) is for compound 5 movement and plot (B) is for compound 7 movement
; Article no.JPRI.66801 compounds 5 and 7 can be seen in Fig. 5 (A) respectively. By superposing the ligand on its reference structure throughout simulation period, we can observe that the conformational changes RMSD for compounds 5 and 7 in Fig. 5 are consistent with ligand movement RMSD for . Based on these results, it is obvious that compound 7 may have more capacity than do compound 5 to bind interface residues of receptor binding domain for square deviation (RMSD) of ligand movement as a function of simulation (B) is for compound 7 movement
Prediction of chemical and with outputs of docking and molecular dynamics studies are supporting the possible capacity of compound 7 to inhibit cellular entry of SARS-CoV the outputs of this virtual screening study may need further in vitro and in vivo evaluation.

CONSENT
It is not applicable. CoV-2. However, the outputs of this virtual screening study may need further in vitro and in vivo evaluation.

ETHICAL APPROVAL
As per international standard or university standard ethical approval has been collected and preserved by the authors.