Molecular Docking Reveals Ivermectin and Remdesivir as Potential Repurposed Drugs Against SARS-CoV-2

SARS-CoV-2 is a newly emerged coronavirus that causes a respiratory disease with variable severity and fatal consequences. It was first reported in Wuhan and subsequently caused a global pandemic. The viral spike protein binds with the ACE-2 cell surface receptor for entry, while TMPRSS2 triggers its membrane fusion. In addition, RNA dependent RNA polymerase (RdRp), 3′–5′ exoribonuclease (nsp14), viral proteases, N, and M proteins are important in different stages of viral replication. Accordingly, they are attractive targets for different antiviral therapeutic agents. Although many antiviral agents have been used in different clinical trials and included in different treatment protocols, the mode of action against SARS-CoV-2 is still not fully understood. Different potential repurposed drugs, including, chloroquine, hydroxychloroquine, ivermectin, remdesivir, and favipiravir, were screened in the present study. Molecular docking of these drugs with different SARS-CoV-2 target proteins, including spike and membrane proteins, RdRp, nucleoproteins, viral proteases, and nsp14, was performed. Moreover, the binding affinities of the human ACE-2 receptor and TMPRSS2 to the different drugs were evaluated. Molecular dynamics simulation and MM-PBSA calculation were also conducted. Ivermectin and remdesivir were found to be the most promising drugs. Our results suggest that both these drugs utilize different mechanisms at the entry and post-entry stages and could be considered potential inhibitors of SARS-CoV-2 replication.

INTRODUCTION SARS-CoV-2 emerged in 2019 as the causative agent of a pneumonia outbreak in Wuhan, Hubei Province, China (Zhu et al., 2020). The disease outbreak spread globally, causing a pandemic with dozens of millions laboratory-confirmed cases (WHO, 2020) with more folds of the infections could be passed undetected (Hosein et al., 2020). The disease has resulted in more than 1,500,000 deaths as of 12 December 2020 (WHO, 2020).
The S protein is proteolytically cleaved by type-II transmembrane serine protease (TMPRSS2) into S1 and S2 subunits (Hoffmann et al., 2020). The former subunit binds to the host cell surface receptor, while the latter is responsible for the fusion of the viral envelope and the cell membrane. The M protein is one of the most abundant envelope proteins. It plays an important role in determining the morphology of the virus. The E protein is present in a small amount on the envelope; however, it is important for the assembly and release of the virus. The N protein binds to the viral genome and forms the nucleocapsid of the virus (Abdel-Moneim and Abdelwhab, 2020). The replicase proteins encode the papain-like protease (PLpro) and the serine-type protease or main protease (Mpro) (Ziebuhr et al., 2000;Mielech et al., 2014). In addition, many other nsps, including RNAdependent RNA polymerase (RdRp; nsp12) (Xu et al., 2003), RNA helicase (nsp13) (Ivanov et al., 2004), N7 MTase and 3 -5 exoribonuclease (nsp14) (Eckerle et al., 2010), form the replicase-transcriptase complex (RTC), which is essential for RNA replication and transcription. The accessory proteins are also involved in viral replication and pathogenesis (Zhao et al., 2012;Fehr and Perlman, 2015).
Although extensive research has been conducted on SARS-CoV-2, no approved drug is available against the virus so far. Hence, there is an urgent need to develop sensitive and effective antiviral agents against zoonotic coronaviruses, including SARS-CoV-2. The registration of new antiviral drugs takes a long time; therefore, many studies have been conducted to screen the efficacy of already approved drugs against SARS-CoV-2. Many drugs have been found to possess potential antiviral activity against SARS-CoV-2. These include old antimalarial drugs (chloroquine phosphate, chloroquine, and hydroxychloroquine) , an anthelmintic drug (ivermectin) (Caly et al., 2020), viral RNA polymerase inhibitors (remdesivir and favipiravir) (National Health Commission of the People's Republic of China, 2020;Wang M. et al., 2020) and viral protease inhibitors (Mugisha et al., 2020). Although some of such drugs possess know antiviral action (Gurung, 2020;Shannon et al., 2020;Wang M. et al., 2020;Yin et al., 2020), however, their possible antiviral effects on different viral proteins are not fully understood. Therefore, the present study aimed to perform molecular docking in order to assess the binding affinity of different viral proteins to different drugs with potential antiviral activities against SARS-CoV-2.

METHODOLOGY Protein Retrieval and Preparation
The 3D structures of recently identified SARS-CoV-2 proteins, namely the S glycoprotein (PDB ID = 6VXX) (Walls et al., 2020), RdRp (PDB ID = 67M1) (Gao et al., 2020), Mpro (PDB ID = 6Y2E) , PLpro (PDB ID = 6W9C) (Osipiuk et al., 2020), and the N protein (PDB ID = 6VYO) (Kang et al., 2020), were obtained from the RCSB Protein Data Bank 1 . On the other hand, the 3D structure of the viral M protein was not available; therefore, structural protein sequences of SARS-CoV-2 were downloaded from the NCBI Protein database (Accession No. QJA17755). Homology modeling of the viral proteins was performed using the SWISS-MODEL server 2 with default settings. The M protein sequence of SARS-CoV-2 was entered in FASTA format, and the 3D homology model was retrieved from the SWISS-MODEL server as a PDB file and used for the docking process. Similarly, because of the lack of experimental 3D structure of the non-structural protein Nsp14 the sequence of SARS-CoV-2 ExoN/nsp14 (P0DTD1) was used to build a 3D homology model in the SWISS-MODEL web server based on the x-ray structure of Nsp14 from SARS-CoV (PDB ID: 5C8S, chain B) (Gurung, 2020).
The 3D structure of human ACE-2 was downloaded from the RCSB Protein Data Bank (see text footnote 1) (PDB ID = 1R42). However, 3D x-ray crystallographic data of TMPRSS2 were not available; therefore, the sequence of human TMPRSS2 (O15393) was retrieved from UniProt (UniProt Consortium, 2019), and loaded into the SWISS-MODEL server (see text footnote 2) with default settings to create three different 3D homology models of the protein. The top-ranked homology model built using the serine protease hepsin as the template (PDB ID = 5CE1) was subjected to protein preparation and optimization using the default protein preparation protocol in the Molegro Virtual Docker (MVD) software. Finally, the verified homology model of TMPRSS2 with good quality was used for molecular docking studies.

Ligand Preparation
The 2D structures of all the drugs used in the study (chloroquine, hydroxychloroquine, ivermectin, remdesivir, favipiravir, lopinavir, and camostat), human ACE-2 cell receptor and cellular serine protease TMPRSS2 were compiled using ChemDraw, and the 3D structures of all the drug ligands were constructed using Chem3D ultra 15.0 software (Molecular Modeling and Analysis; Cambridge Soft Corporation, United States 2014). Then, the constructed structures were energetically minimized using MOPAC (semi-empirical quantum mechanics) (Job Type with 100 iterations and minimum RMS gradient of 0.01) and saved as an MDL Molfile ( * .mol).

Molecular Docking
The 3D protein structures of all the proteins being studied were loaded onto the MVD 6.0 (2013) platform for the docking process. Potential binding sites (referred to as cavities) were identified using the built-in cavity detection algorithm of MVD. For each PDB file, protein preparation was performed using the default parameters in MVD before conducting the docking experiment. Subsequently, the docking process between different ligands and active sites of different protein structures was performed using the MolDock score as the scoring function of MVD with a grid resolution of 0.30 Å. The number of runs for each docking process was 10. Moreover, the maximum iterations were 2000, with an energy threshold of 100 Kcal/mol. The best conformations for each docking process were selected on the basis of the lowest docked binding energy.

MD Simulation and MM-PBSA Calculation
To validate the accuracy of the ivermectin and remdesivir results, a short molecular dynamics (MD) simulation for each of the predicted receptor-ligand complexes was conducted. All simulations were conducted using GROMACS 2018 (Abraham et al., 2015), employing AMBER 99SB-ILDN force field (Maier et al., 2015). For every ligand-receptor complex, all receptor nonterminal missing loops/atoms were modeled using MODELLER (Sali and Blundell, 1993). Then missing hydrogen atoms were added using chimera (Pettersen et al., 2004) built-in method, putting special consideration on histidine residues to add their hydrogen atoms based on each histidine microenvironment. The complex was placed in an isotonic box (Main Protease and M-Protein were placed in a triclinic box, while the rest were placed in a dodecahedron cell). The box contains a neutralizing number of Na + and Cl − ions according to the protein charge. Topology parameters for the ligands were built using ACPYPE (Sousa da Silva and Vranken, 2012) and ANTECHAMBER (Wang et al., 2006) to generate general amber force field (GAFF) parameters. Then the regular steps of energy minimization, equilibration (NVT and NPT), then MD simulation with 2 femtosecond integration steps for 30 ns were conducted. The output trajectory was PBC corrected and the system was fitted to its start position based on the receptor's backbone before further analysis was performed. The analysis included: (i) plotting the RMSD (root-mean-square deviation) of the ligand after fitting, (ii) plotting the number of hydrogen bonds between the ligand and the receptor, and (iii) MM-PBSA energy calculation. The first two measures were conducted using GROMACS suite, and the MM-PBSA energy was calculated using the g_MM-PBSA package (Baker et al., 2001;Kumari et al., 2014). For MM-PBSA calculation, the SAV model for calculation of the non-polar solvation energy was adopted, taking a sample every 1 ns. A ligand that could maintain its RMSD within 1 nm was considered fixed, less than 2 nm to be stable, but a ligand that possessed more than 2 nm RMSD to be non-stable.

Binding Interactions of Selected Drugs With SARS-CoV-2 S Glycoprotein
The binding energies of the selected drugs (Supplementary Figure 1) with the S glycoprotein (PDB ID = 6VXX) were studied. The docking results are given in terms of the MolDock score, interaction energy, H-bond energy and interacting amino acid residues present at the predicted active site of the protein ( Table 1). Ivermectin showed the highest binding affinity to the predicted active site of the protein (MolDock score −140.584) and protein-ligand interactions (MolDock score −139.371). Moreover, it formed four H-bonds with Thr307, Glu309, Ile312, and Asn953 amino acid residues present at the predicted active site of the protein. In contrast, favipiravir showed the lowest binding affinity to the predicted active site of the protein (MolDock score −54.595) and proteinligand interactions (MolDock score −68.539). Moreover, it formed one H-bond with the Arg1014 amino acid residue Docking parameters was conducted using Molegro Virtual Docker 6.0.
Frontiers in Microbiology | www.frontiersin.org Frontiers in Microbiology | www.frontiersin.org present at the predicted active site of the protein. Remdesivir showed considerable binding affinity (MolDock score −111.007) and protein-ligand interactions (MolDock score −122.699). It formed one H-bond with Leu303, Glu309, and Asn953. The 3D structural views of the ligand-binding site interactions are provided in Figure 1 and Table 1. By using open protein form, similar binding affinities were detected in all tested drugs, however, some variations in the interacting amino acids were detected (Figure 1 and Table 1).

Binding Interactions of Selected Drugs With SARS-CoV-2 RdRp
The docking results of the selected drugs with SARS-CoV-2 RdRp (PDB ID = 6M71) are provided in    binding affinity (MolDock score −93.634) and protein-ligand interactions (MolDock score −107.868), however, formed two H-bonds only with Ser158. In contrast, favipiravir did not show considerable binding to Mpro (Table 3 and Figure 4).

Binding Interactions of Selected Drugs With SARS-CoV-2 PLpro
The docking results of the selected drugs with SARS-CoV-2 PLpro (PDB ID = 6W9C) are provided in  Figure 5.

Binding Interactions of Selected Drugs With SARS-CoV-2 M Protein
The binding affinities of the selected drugs to the M protein were studied. The docking results are provided in terms of the MolDock score, interaction energy, H-bond energy, and interacting amino acid residues present at the predicted active site Mpro, Main protease (PDB id = 6Y2E); PLpro, Papain like protease (PBD id: 6W9C).  The 3D structural views of the ligand-binding site interactions are provided in Figure 6.

Binding Interactions of Selected Drugs With SARS-CoV-2 N Phosphoprotein
The docking results of the selected drugs with the SARS-CoV-2 N phosphoprotein (PDB ID = 6VYO) are provided in Table 4. Ivermectin showed the highest binding affinity to the predicted active site of the protein (MolDock score  It formed hydrogen bonds with Asn75, Asn77, His145, and Asn154 (Figure 7).

Binding Interactions of Selected Drugs With Human ACE-2 Protein
MLN-4760, an ACE-2 inhibitor, was used as a positive control in the docking process. The docking results (     Figure 2).

RMSD and MM-PBSA Energy
All the studied receptor -ligand complexes showed ligand heavy atoms RMSD within the stable range, except RdRp/ivermectin, spike/remdesivir, NP/remdesivir, and PLpro with either lopinavir or remdesivir (Supplementary Figure 3). The free binding energy of the spike protein (open) was higher in ivermectin (−398.536 kJ/mol) than remdesivir (−232.973 kJ/mol). In the remdesivir/spike structure complex, although the simulation showed high RMSD value of the ligand heavy atoms, this high value is related to the major steric changes the protein underwent. The ligand managed to follow the change, keeping a stable number of two hydrogen bonds in the second half of the trajectory ( Table 7 and Supplementary  Figure 3 and Figures 8A,B). Although there was no significant difference between free binding energy for SARS-CoV-2 RdRp with ivermectin (−179.472 kJ/mol) and remdesivir (−345.437 kJ/mol), the ivermectin/RdRp complex was not stable (Supplementary Figure 3). The free binding energy for ivermectin with ExoN/NSP14 (−418.894 kJ/mol), was higher than that with remdesivir (−362.872 kJ/mol). The free binding energy of lopinavir/Mpro complex was −324.160 kJ/mol while lopinavir/PLpro complex showed low free binding energy (−134.961 kJ/mol). The lopinavir/PLpro complex was not successful since lopinavir left the pocket totally (Figures 8C,D).
Similarly, ivermectin showed a higher free binding energy (−345.675 kJ/mol) with Mpro, and lower binding energy of −221.257 with PLpro. In contrast, remdesivir showed similar free binding energy for Mpro and PLpro; −222.664 kJ/mol and −209.861 kJ/mol, respectively. However, in the remdesivir/PLpro complex, the remdesivir moved from one side of the pocket to the other side (Figures 8E,F). The free binding energy of M protein with ivermectin was the highest among all the complexes (−516.656 kJ/mol) and it also showed higher free binding energy with NP was −334.190 than remdesivir ( Table 7 and Supplementary Figure 2). Remdesivir free binding energy with NP protein (−169.744 kJ/mol) was lower than that detected with ivermectin ( Table 7 and Supplementary  Figure 2). The remdesivir was stabilized in the NP binding pocket in the beginning by 1-4 simultaneous hydrogen bonds, in addition to the hydrophobic attraction (force). Then, during the period t = 21 till t = 30, the pocket was obliterated. However, the ligand was kept near to its pocket, using 1-3 hydrogen bonds, until the pocket was opened. The ligand was then restored to its original position ( Table 7 and Figures 8G-J).
The free binding energy of human ACE-2 with ivermectin and remdesivir were −227.529 and −209.396, respectively. The remdesivir/human ACE-2 complex, just after the initial part of the simulation, showed a tendency to stabilize in place with 4-5 hydrogen bonds (Table 7 and Supplementary  Figure 3). On the other hand, the free binding energy of human TMPRSS2/camostat complex was −266.882. High free binding energy of TMPRSS2/was also recorded with both ivermectin, and remdesivir, −293.245 and −249.480, respectively ( Table 7).

DISCUSSION
Both ivermectin and remdesivir showed high binding affinity to different viral proteins and seem to be potential drugs against SARS-CoV-2; however, laboratory and clinical trials are needed, particularly for ivermectin. Both these drugs possess multidisciplinary actions. Ivermectin showed high binding affinity to the viral S protein as well as the human cell surface receptors ACE-2 and TMPRSS2. In agreement to our findings, ivermectin was found to be docked between the viral spike and the ACE2 receptor (Lehrer and Rheinstein, 2020). The activation of the S protein by TMPRSS2 can make cathepsin L activity and low pH unnecessary for the viral envelope to fuse with the endosomal membrane (Glowacka et al., 2011). The molecular docking of ivermectin with TMPRSS2 suggested an important role of ivermectin in inhibiting the entry of the virus into the host cell, probably by increasing the endosomal pH. Moreover, it efficiently binds to both Mpro and binds also but to lesser extent to PLpro of SARS-CoV-2; therefore, it might also play a role in preventing the post-translational processing of viral polyproteins. The highly efficient binding of ivermectin to the viral N phosphoprotein and M protein is suggestive of its role in inhibiting viral replication and assembly. Accordingly, it may also be involved in inhibiting nuclear transport. In addition to the efficient binding of ivermectin to the N protein, previous studies have revealed that ivermectin inhibits IMPα/β1-mediated nuclear import of the N protein (Rowland et al., 2005;Timani et al., 2005;Wagstaff et al., 2012;Tay et al., 2013;Caly et al., 2020;Yang et al., 2020). A recent in vitro study revealed ∼5000fold reduction of SARS-CoV-2 RNA at 48 h in Vero cells when ivermectin was added to the cells 2 h post-infection (Caly et al., 2020). The highly efficient binding of ivermectin to nsp14 confirms its role in inhibiting viral replication and assembly. It is well known that nsp14 is essential in transcription and replication. It acts as a proofreading exoribonuclease and plays a role in viral RNA capping by its methyl transferase activity (Ma et al., 2015).
Remdesivir showed high affinity to spike but formed unstable complex, however, it showed considerable high affinity to both TMPRSS2 and ACE-2 might denote its possible roles in blocking cellular receptor necessary for viral entry in addition of inhibiting TMPRSS2 induced membranes' fusion required for the SARS-CoV-2 replication (Glowacka et al., 2011;Hoffmann et al., 2020). The affinity of remdesivir to bind human TMPRSS2 was also suggested in a previous study (Wu C. et al., 2020).
Remdesivir is an adenosine analog, it inhibits RNA strand elongation, thereby inhibiting viral replication  (Warren et al., 2016). In addition, it showed high binding affinity to RdRp, as confirmed in the present study. A recent study found that remdesivir able to bind RdRp (Wu C. et al., 2020). In addition, SARS-CoV-2 RdRp in the apo form or in association with the RNA template primer can bind to remdesivir (Yin et al., 2020). The high binding affinity of remdesivir to both Mpro and PLpro of SARS-CoV-2 is suggestive of its potential as a protease inhibitor. Our assumption was confirmed by a recent in silico study (Deshpande et al., 2020), which reported that remdesivir efficiently binds to Mpro. In addition, remdesivir has been found to be more potent than lopinavir/ritonavir both in vitro and in MERS-CoV-infected mice (Sheahan et al., 2020). Remdesivir was found to have high affinity to SARS-CoV-2 main protease (Hall and Ji, 2020).
Remdesivir also showed high affinity to M and nsp14 proteins of SARS-CoV-2, further suggesting its role in inhibiting different steps of viral replication. M and S protein as well as M and N protein interactions are required for viral component assembly in host cells (Siu et al., 2008).
Although lopinavir was reported to show high affinity to SARS-CoV-2 main protease (Kumar et al., 2020), in the current study, it formed unstable complex which also agree with a previous study that found it unable to bind to Mpro and PLpro (Choy et al., 2020;Wu C. et al., 2020).
The binding affinities of chloroquine, hydroxychloroquine, and favipiravir were lower than those of ivermectin and remdesivir for all the tested proteins. However, chloroquine and hydroxychloroquine showed considerably high binding affinities to Mpro, the M protein, ACE-2, and TMPRSS2. Previous studies have revealed that hydroxychloroquine possesses wide-range antiviral activity (Savarino et al., 2006;Yan et al., 2013). In accordance with these findings, it was found that hydroxychloroquine exerts its antiviral effect against SARS-CoV by altering the pH and interfering with the glycosylation of the ACE-2 receptor (Vincent et al., 2005); this is considered to be the major antiviral action of hydroxychloroquine (Savarino et al., 2006). Interestingly, hydroxychloroquine interferes with glycosyltransferases (sugar-modifying enzymes) and inhibits quinone reductase 2 (Kwiek et al., 2004) and the biosynthesis of sialic acid (Vincent et al., 2005). In addition, it interferes with the fusion of the virus and the cellular endosome by altering the endosomal pH (Vincent et al., 2005). Although chloroquine phosphate has been found to successfully inhibit SARS-CoV-2 replication in different clinical trials Yin et al., 2020), this drug did not show binding activity to the target proteins tested in the present study. Similarly, favipiravir is a purine analog that inhibits the elongation phase of RNA synthesis. It has shown promising effects in COVID-19 patients (Cai et al., 2020). Moreover, an in vitro study revealed that compared with chloroquine or remdesivir, only a high concentration of favipiravir (EC 50 ? = ?61.88 µM) could inhibit SARS-CoV-2 replication . However, in the present study, it did not show good binding affinity to any of the tested proteins.
Chloroquine, hydroxychloroquine, and favipiravir could exhibit antiviral effects through different proteins or intermediate protein-protein interactions; however, these were not screened in the present study.
In conclusion, both ivermectin and remdesivir could be considered potential drugs for the treatment of COVID-19. Ivermectin efficiently binds to the viral S protein as well as the human cell surface receptors ACE-2 and TMPRSS2; therefore, it might be involved in inhibiting the entry of the virus into the host cell. It also binds to Mpro and PLpro of SARS-CoV-2; therefore, it might play a role in preventing the post-translational processing of viral polyproteins. The highly efficient binding of ivermectin to the viral N phosphoprotein and nsp14 is suggestive of its role in inhibiting viral replication and assembly. Remdesivir may be involved in inhibiting postentry mechanisms as it shows high binding affinity to N and M proteins, PLpro, Mpro, RdRp, and nsp14. Although the results of clinical trials for remdesivir are promising (Beigel et al., 2020;Wang Y. et al., 2020), similar clinical trials for ivermectin are recommended. Both these drugs exhibit multidisciplinary inhibitory effects at both viral entry and postentry stages.

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 author.

AUTHOR CONTRIBUTIONS
ASA-M contributed conception of the study and critically revised the manuscript. AFE conducted the molecular docking and wrote the first draft of the manuscript. AAA performed the molecular dynamics simulation. ASA-M, AFE, and AAA shared analyzing and discussing the results. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported by the Taif University Researchers Supporting Project Number (TURSP-2020/11), Taif University, Taif, Saudi Arabia. The molecular dynamics simulation was performed in the National Supercomputing Centre (NSCC) Singapore, under projects numbers 32002093, 12000641, and 12000569.