A large-scale computational screen identifies strong potential inhibitors for disrupting SARS-CoV-2 S-protein and human ACE2 interaction.

SARS-CoV-2 has infected millions of individuals across the globe and has killed over 2.7 million people. Even though vaccines against this virus have recently been introduced, the antibody generated in the process has been reported to decline quickly. This can reduce the efficacy of vaccines over time and can result in re-infections. Thus, drugs that are effective against COVID-19 can provide a second line of defence and can prevent occurrence of the severe form of the disease. The interaction between SARS-CoV2 S-protein and human ACE2 (hACE2) is essential for the infection of the virus. Thus, drugs that block this interaction could potentially inhibit SARS-CoV-2 infection into the host cells. To identify such drugs, we first analyzed the recently published crystal structure of S-protein-hACE2 complex and identified essential residues of both S-protein and hACE2 for this interaction. We used this knowledge to virtually dock a drug library containing 4115 drug molecules against S-protein for repurposing drugs that could inhibit binding of S-protein to hACE2. We identified several potential inhibitors based on their docking scores, pharmacological effects and ability to block residues of S protein required for interaction with hACE2. The top inhibitors included drugs used for the treatment of hepatitis C (velpatasvir, pibrentasvir) as well as several vitamin D derivatives. Several molecules obtained from our screen already have good experimental support in published literature. Thus, we believe that our results will facilitate the discovery of an effective drug against COVID-19. Communicated by Ramaswamy H. Sarma.

Introduction

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) belongs to the coronaviridae family—a large family of single-stranded enveloped RNA viruses (Enjuanes et al., 2006; Li, 2016; Perlman & Netland, 2009). This virus is responsible for the current coronavirus 2019 disease (COVID-19) outbreak. SARS-CoV-2 is the seventh reported human-infecting coronavirus out of which SARS-CoV, MERS-CoV and SARS-CoV-2 can cause serious disease. In comparison, the viruses HKU1, NL63, OC43 and 229E of the same family cause only mild disease in humans (Corman et al., 2018). The epidemic, which started in 2019, has caused more than 24 million infections worldwide and has resulted in >2.7 million deaths so far.

An envelope-anchored spike protein (S-protein) mediates coronavirus entry into host cells by first binding to a host receptor and then fusing viral and host membranes. A defined receptor-binding domain (RBD) of S-protein specifically recognizes its host receptor angiotensin-converting enzyme 2 (ACE2). Different lines of research have shown that whether the host would be susceptible to SARS-CoV infection is primarily determined by the affinity between the viral RBD and human ACE2 (hACE2) in the initial viral attachment step (Fang, 2015; Li, 2004; Li et al., 2003, 2005; Mccray, 2007; Moore et al., 2004; Qu, 2005). Thus, inhibition of the interaction between S-protein and hACE2 presents an attractive solution for preventing SARS-CoV-2 infection in humans. Besides, the virus depends on a cycle of infection and replication to multiply in numbers inside the host. This leads to severe disease and death in some individuals. Thus, preventing S-protein from binding to hACE2 could also disrupt this cycle and can reduce the chances of severe disease in patients.

Drugs such as hydroxychloroquine have been shown to inhibit the binding of S-protein and hACE2 in vitro (Liu, Zhou, et al., 2020) and, thus, have been used in the treatment in the fight against COVID-19 (Gautret et al., 2020; Rathi et al., 2020). There are, however, still questions over the effectiveness of this drug against COVID-19 as well as over its adverse side effects (Costedoat-Chalumeau et al., 2007; Cotroneo et al., 2007; Joyce et al., 2013; Nord et al., 2004; Sharma et al., 2016, 2020; Zhao et al., 2018). In addition, a recent clinical trial has focused on using a soluble form of the molecule hACE2 that will sequester the majority of S-proteins to prevent SARS-CoV-2 from binding to the host cell hACE2 (Monteil et al., 2020). Several vaccines against SARS-CoV-2 have undergone human clinical trials (Thanh Le et al., 2020; Liu, Cao, et al., 2020) and have received emergency use authorization.

A recent study measured antibody binding response to S-protein and nucleocapsid protein and the neutralization potency against SARS-CoV-2 (Ibarrondo et al., 2020; Long et al., 2020; Seow et al., 2020). They observed a decline in IgM and IgA binding responses after 20–30 days post onset of syndromes, and 2- to 23-fold decrease in neutralizing antibody during an 18–65 day follow-up period (Seow et al., 2020). In two other studies, similar reduction in IgG titer was observed (Ibarrondo et al., 2020; Long et al., 2020). Although several other questions regarding decline of antibody needs to be answered, the early observations can have important implications in antibody protection against re-infection of SARS-CoV-2 and the durability of vaccine protection.

In this regard, the repurposing of drugs that are already approved for human use presents an attractive solution to look for an effective inhibitor of S-protein-hACE2 binding, since these drugs can readily be used for the treatment of COVID-19 in patients (Elmezayen et al., 2020; Kiplin Guy et al., 2020; Pawar, 2020; Senanayake, 2020). These drugs can be used in conjunction with the vaccines to reduce chances of severe disease in infected individuals. To identify such drugs, we performed a large-scale computational screening to find inhibitors of S-protein-hACE2 binding from a database of drugs already approved for human use that contained ∼4200 molecules. The workflow involved primary screening based on rigid docking of human approved drugs on the recently published crystal structure of S-protein-hACE2 interaction (Wang et al., 2020). Drugs were then ranked by docking score and top 70 drugs were checked for their ability to block S-protein residues required for S-protein-hACE2 interaction. We further confirmed the stability of a few selected top protein-drug complexes using Molecular Dynamics (MD) simulations studies (Scheme 1). We ended up with a small number of strong potential inhibitors from a large number of drug molecules. This included drugs used against hepatitis C virus like velpatasvir, vitamin D derivatives like ergocalciferol, drugs used to prevent asthma symptoms like zafirlukast. These molecules can quickly be tested for their effectiveness against S-protein binding to hACE2 in the laboratory and can lead to a quick and effective treatment against COVID-19.

Scheme 1.

Workflow for identification of potential inhibitors against SARS-CoV-2 S protein and human ACE2 interaction via structure-based virtual screening. The structure of the S-protein used in Autodock and MD simulations was obtained from the crystal structure of S-protein-hACE2 complex (PDB ID: 6LZG)

Methods

The structure of novel coronavirus spike receptor-binding domain complexed with its receptor hACE2 submitted to RCSB PDB (PDB ID: 6LZG) was used (Wang et al., 2020). AutoDock 4.2 was used for virtual screening. AutoDockTools was used for visualizing the binding of a drug to the protein target (S protein) and to identify interacting residues in protein target (Morris et al., 2009). The master database which was used for screening consisted of all drugs present in SuperDrug version 2.0 database (around 4200 drugs) (Goede et al., 2005; Siramshetty et al., 2018) and all approved anti-viral drugs from Merged and Unified data (approximately 900 drugs; https://sites.google.com/view/mud-data/home) that were not present in the Superdrug database. Before docking, both protein target and drug molecules were prepared as described below. Semiflexible docking protocol was followed where protein target was kept rigid, and drug molecules were flexible so that they could attain a degree of freedom torsions bridged by the rotational parameter.

Preparation of drug molecules

Preparing the ligand involved ensuring that its atoms are assigned the correct AutoDock 4 atom types, adding Gasteiger charges if necessary, merging non-polar hydrogens, detecting aromatic carbons if any and setting up the ‘torsion tree’. Since the database contained drug molecules in 3 D structure-data file (SDF) format, they were first converted to PDB format using Open Babel 3.0.0 (O’Boyle et al., 2011), to ensure compatibility with Autodock. Using the script prepare_ligand4.py, available in the MGLTools package, the drug molecules were converted to PDBQT format. This format was similar to PDB format but also includes partial charges (‘Q’) and AutoDock 4 (AD4) atom types (‘T’). The atom types included in our study were ‘A’, ‘C’, ‘HD’, ‘N’, ‘OA’, ‘SA’, ‘F’, ‘Cl’, ‘Br’, ‘I’, ‘S’, ‘P’ and ‘Zn’. The AD4 atom types and their parameters could be found in the ‘AD4_parameters.dat’ file, which was included in the source code of the AutoDock 4 distribution.

Preparation of protein target molecule

The co-crystallized structure of hACE2 and receptor-binding domain (RBD) of S protein was downloaded from RCSB PDB website (PDB ID: 6LZG). InterfaceResidues.py script (developed by Mamoru Suzuki at Institute for Protein Research, Osaka University and available on pymolwiki.org) was used in PyMOL to identify interface residues between chain A (hACE2) and chain B (S protein) (Figure 1) (Delano, 2002). From the PDB file, the chain A was removed and the resulting PDB file contained only the B chain (S protein). The protein target molecule was further processed by removing water molecules, adding all the hydrogen atoms, merging non-polar hydrogen atoms and removing any other ligand/heteroatom present using AutoDock Tools. The charges were assigned using the Gasteiger method, and finally, the protein was saved in PDBQT format.

Figure 1.

Co-crystallised structure of human ACE2 and SARS-CoV-2 S protein. (A) Interacting parts of the respective proteins are shown in stick representation and are colored in purple and orange in hACE2 and RBD of SARS-CoV-2 S protein, respectively. (B) List of Interacting residues of Chain A (hACE2) and Chain B (RBD of S protein). (Visualized using PyMol).

Autogrid4 and Autodock4

The S-protein region that contained residues important for interaction with hACE2 was chosen as the docking site for drug molecules. Autogrid4 parameter was used to attain a rigid grid box of 82 × 126 × 96 with a spacing of 0.375 Å and set around the S-protein site binding to hACE2 protein (). This grid space defined the search space for Autodock to find the binding conformation of a ligand with maximum binding affinity. In the next step, AutoDock4 with the Lamarckian Genetic Algorithm was run to obtain the best docking conformation as well as for calculation of docking score for each ligand. In the process, DLG files were generated that stored all results obtained from docking including docking score that reflected binding affinity of the ligand to the S protein. Lower docking score corresponded to stronger binding interaction of the ligand with the protein molecule. Docking score for all ligand molecules in our analysis was plotted using a custom R script, and top molecules were taken forward for secondary analysis (Figure 2(A)).

Figure 2.

Summary of screening of drugs against SARS-CoV-2 S protein. (A) Docking score distribution of all drugs (in descending order) with a threshold of −5 kcal/mol. Drugs with docking score less than −8.4 Kcal/mol were considered for secondary analysis. (shown by ellipse in dotted red line.) (B) Barplot of the docking score of selected drugs. Values at the top of bins denote the number of residues of S protein to which both human ACE2 receptor and drugs are interacting. (C) Secondary screening results of selected drugs. The second column shows the number of S-protein residues to which both hACE2 and each of the mentioned drugs interact. The highlighted orange colored cell shows the S protein residues (with residue numbers on top) to which both hACE2 and drugs interact. Hydroxychloroquine (HCQ) and azithromycin were chosen as molecules for comparison with the other potential drugs.

Docking validation studies

The docking protocol was validated using two methods: (i) Re-docking of the ligand set in the master database three times, keeping all the settings, including grid-space, the same. The mean docking score and the standard deviation were calculated from these replicate dockings. It was done to ensure that the inhibitor binds exactly to the active site cleft and shows less variation in docking score from one docking to another. (ii) 631 Decoy ligands similar to top20 drugs found in primary screening were obtained from DUD-E online server (http://dude.docking.org/) (Mysinger et al., 2012) and were docked against the binding site of the spike protein along with top20 drugs. As there was no experimentally verified inhibitor against spike protein, top20 drugs from primary screening were assumed as active drugs. The enrichment plot, receiver operating characteristic (ROC) curve and areas under the ROC curves (AUC) were plotted using the ROCR package in R (Bunn, 2013 ; Sing et al., 2005; Stahl, 2000; Sun and Xu, 2014; Wei et al., 2002). The enrichment plot graphs, at any given percentage of the database ranked by the docking calculation, the number of true ligands actually found divided by the number of ligands expected to be found if they were just chosen at random.

It can also be written as ratio of true positive rate and false positive rate in a given subset of the database. A ROC curve is generated by plotting the rate against the false positive rate at each cutoff which is varied from the highest to the lowest predicted scores. The area under ROC curve is a measure of prediction algorithm performance where 0.5 is random prediction and 1.0 is perfect prediction.

Secondary screening of selected drugs

The secondary analysis involved understanding the pharmacological effects of selected drugs (including side effects, if any, and its method of application) and identifying the number of residues of S-protein to which both drug and hACE2 were interacting upon binding to S-protein. For the latter part, Receptor-Ligand option in Discovery Studio Visualizer (Biovia, 2017) was used.

Molecular dynamics (MD) simulation

Selected top-ranked compounds were further subjected to molecular dynamics (MD) simulation using GROMACS 2020.1 software (Hess et al., 2008). Ligand parameterization was performed using CHARMM General Force Field (CGenFF) web-based tool (https://cgenff.umaryland.edu/). MD simulations were performed using the CHARMM36 force field (Huang & Mackerell, 2013). Visual molecular dynamics (VMD) was used to generate GRO files for the complex (Humphrey et al., 1996). The complex was solvated in dodecahedron water boxes containing transferable intermolecular potential with 3 points (TIP3P) water molecules (Jorgensen et al., 1983). The system first performed 50,000 steps of steepest descent with energy minimization, and then, the minimized system was used to perform simulations using an NVT ensemble (same number of particles, volume and temperature) ensemble and NPT ensemble (same number of particles, pressure and temperature). The system was equilibrated with pressure, temperature, the ligand was restrained whenever protein was also restrained, and the simulation time was set to 50 ns, with a time step of 2 fs at 300 K. The trajectories acquired were investigated using rmsd module of GROMACS 2020.1 and molecular dynamics simulation movies with a smoothing window of size 25 were created using visual molecular dynamics (VMD) (Humphrey et al., 1996). The interaction-free energies of each receptor-ligand complex were evaluated using Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) method (Kollman et al., 2000). In this study, whole MD trajectories were used for the calculation of binding free energies, van der Waals energy, electrostatic energy and polar solvation energy by using the g_mmpbsa subroutine implemented in GROMACS by Kumari et al., and made compatible with GROMACS 2020.1 by Salnikov et al. (Baker et al., 2001; Eisenhaber et al., 1995; Kumari et al., 2014; Salnikov et al., 2009; Wagoner & Baker, 2006).

Results

Inhibiting the interaction between S-protein and hACE2 required a detailed understanding of the residues that facilitate this interaction since the drug molecules need to bind to these residues to prevent this interaction. Thus, we first analyzed the crystal structure of S-protein-hACE2 complex (PDB ID: 6LZG (Wang et al., 2020)) (Figure 1(A)) and identified the essential residues for this interaction (Figure 1(B)). Using Discovery Studio Visualizer (Biovia, 2017), we determined these non-bonding interactions to be mainly hydrogen bonding, electrostatic and hydrophobic interactions. We also observed some steric hindrances, indicating further mutations can still facilitate stronger binding. These results were also consistent with the predictions from Wan et al. (2020).

We then derived the structure of S-protein from the complex crystal structure and used Autodock4 to dock 4115 drug molecules for which three-dimensional structures were available. These drug molecules were obtained from SuperDrug2 database (Goede et al., 2005; Siramshetty et al., 2018) and Merged and Unified Data (https://sites.google.com/view/mud-data/home) (see Methods section for more details). We ranked these molecules according to the docking score for binding to the S-protein (Figure 2(A), Supplementary Table 1) and lower value of docking score indicated stronger binding to S-protein (Supplementary Table 1).

We validated our docking methodology by two different methods. First, we redocked the same database three times and checked whether there were significant variations in docking score for a drug across these trials. We observed a positive correlation between percent change in the variation in docking score across these trials (calculated by ((Stdev/abs[mean])×100)) with mean score (Figure 3(A)). This suggested that the drugs with stronger binding to the S-protein and consequently with more negative docking score were binding specifically to the S-protein binding site. In comparison, drugs with weaker binding had less specific interaction with S-protein.

Figure 3.

Validation studies for docking methodology. (A) Plot of percent change in the standard deviation of docking score with respect to mean score vs. mean docking score. (B) The enrichment factor of docking the 631 decoys along with top20 drugs from the primary screening. (C) ROC curve D) ROC curve on semi-log scale.

Second, we performed docking of decoy molecules along with top 20 drugs obtained from primary screening (we ranked the drugs according to BE1 value in Supplementary Table 1), keeping all the settings same (see Methods sections for more details). We analyzed the results using two approaches. The first is an ‘enrichment plot’ (Figure 3(B)). The screening pipeline's effectiveness is assessed by the magnitude of the maximum enrichment factor and the location of the peak. We observed maximum enrichment of 24.27 at 2% of the database (Figure 3(B)). The second plot is ROC curves (Figure 3(C,D)), which show the relationship between the sensibility and the specificity of a test. For our methodology, we found the value of AUC to be 0.861. These results suggested that our docking method was accurate and reproducible.

Next, we performed a secondary screening to identify the potential inhibitors that could be suitable for use. To do so, we considered the top 70 molecules from the primary screening with the lowest docking scores. Drugs in the top 70 could be categorized into four categories: anti-viral drugs, anti-cancer drugs, vitamin D/A derivatives and other drugs (such as anti-diabetic, antibiotic and drugs used for diseases like asthma). Next, we considered the pharmacological effect of each drug in these categories and discarded drugs with known potentially toxic side-effects such as anti-cancer drugs. These drugs can prove to be effective as in cases of tumor cells as well as virus, cells are pushed to enhance the synthesis of nucleic acids, protein and lipid synthesis and boost their energy metabolism (Ciliberto et al., 2020), but they are likely to have very strong side effects, hence we did not consider these drugs further. We also did not consider drugs only approved for topical use, as these drugs may not always be tested for safety in internal use.

We then picked 2–3 drugs from appropriate categories and looked into the interaction of these drug molecules with the S-protein. This facilitated the identification of residues critical for interaction with hACE2, being blocked by each drug molecule (Figures 4 and 5; Supplementary Figures 2 and 3) through a careful analysis of the three-dimensional docking of these drugs with S-protein (). We further considered drugs with docking score less than −8.4 kcal/mol (as drugs with docking score around this value have shown to be effective (Cheng et al., 2018; Mustafa et al., 2017; Onawole et al., 2018; Ramharack & Soliman, 2017; Sheu et al., 2003)) and blocking at least eight S-protein residues required for interaction with hACE2 protein. This resulted in identification of 10 strong potential inhibitor candidates.

Figure 4.

2D Interaction diagram of few top potential drugs that can inhibit human ACE2 and SARS-CoV-2 binding. (A) Pibrentasvir; (B) Velpatasvir; (C) Lonafarnib; (D) Zafirlukast; (E) Pranlukast; (F) Ergocalciferol.

Figure 5.

2D Interaction diagram of few top potential drugs that can inhibit human ACE2 and SARS-CoV-2 binding. (A) Gliquidone; (B) Gitoxin; (C) Ledipasvir; (D) Calcitriol; (E) Hydroxychloroquine; (F) Azithromycin.

Among the top 10 molecules identified through this process, two molecules, Pibrentasvir and Velpatasvir, both showed strongest binding to S-protein with docking score of −9.66 kcal/mol. These drugs also blocked 11 and 10 hACE2 binding sites of S-protein, respectively (Figures 2(B,C) and 4(A,B)). Notably, these drugs are primarily used as inhibitors of the hepatitis C virus (HCV) non-structural protein 5 A (NS5A) replication complex. Although the exact mechanism of action of Velpatasvir has not yet been completely determined, it has been reported to bind to domain I of the NS5A protein, inhibiting the activity of the NS5A protein, which results in the disruption of the viral RNA replication complex, blockage of viral HCV RNA production and inhibition of viral replication (Asselah et al., 2018; Bonaventura et al., 2016; Heo & Deeks, 2018). Zafirlukast and pranlukast, both being leukotriene receptor antagonists used for treating chronic asthma (Adkins & Brogden, 1998; Hur et al., 2018; Keam et al., 2003; Yan et al., 2019), had docking score of −8.71 kcal/mol and −8.61 kcal/mol. These molecules blocked 11 hACE2 binding sites of S-protein (Figures 2(B,C) and 4(D,E)).

Multiple recent epidemiological studies have identified a negative correlation between vitamin D level in individual populations and SARS-CoV-2 infection and mortality (Daneshkhah et al., 2020; Ilie et al., 2020; Mitchell, 2020), although no exact mechanism has been proposed yet. Our screening identified a role for several vitamin D molecules in inhibiting S-protein-hACE2 interaction, thus providing a hint toward how vitamin D can reduce SARS-CoV-2 infection. For example, ergocalciferol, a vitamin D2, showed strong binding (−8.52 kcal/mol, blocking 8 residues) with S protein in our study. The same was observed with another vitamin D molecule, calcitriol. Further, our top 20 drugs go in line with findings from many other independent studies involving virtual screening of ligands against spike protein of SARS-CoV-2 (Kalhor et al., 2020; Vital De Oliveira et al., 2020).

We further compared our results with hydroxychloroquine (a synthetic quinoline derivative) and azithromycin. Docked structures of these two drugs are shown in Figure 5(E,F). The docking score of these molecules is −4.83 kcal/mol and −4.32 kcal/mol, respectively, much lower than the drug candidates that we report here. Further, each of these molecules showed interactions with 7 and 8 residues, respectively, of S-protein critical for its binding to hACE2, which was again lower than the numbers observed for top drugs (Figure 2(B,C)). These observations in regards to hydroxychloroquine are in accordance with an in-vitro experiment performed by Gordon et al. (2020) that demonstrated low efficacy of hydroxychloroquine against SARS-CoV-2. The same study found progesterone to have a better anti-viral activity but only to a small extent. In our study, Levonorgestrel, a synthetic progesterone, showed strong binding with a docking score of −7.88 kcal/mol, but it was able to block just 1 residue of S-protein necessary for interaction with hACE2 (Supplementary Figures 2(I) and Li, 2016). This could perhaps explain the relatively lower efficacy of progesterone in inhibiting SARS-CoV-2 infection (Gordon et al., 2020).

Next, we analyzed different types of non-covalent interactions between a drug molecule and the S-protein for the top 19 drugs in our study, along with hydroxychloroquine and azithromycin. As shown in Figures 4–6, Supplementary Figures 2 and 3, we found 15 types of non-covalent interactions, including some unfavorable interactions. Almost all the drugs had more than one hydrogen bond and even more than five in a few cases. Our top-performing drugs, Velpatasvir and Pibrentasvir, showed a total of 16 and 26 non-bonding interactions, respectively. Velpatasvir formed 5 conventional and carbon hydrogen bonds with residues Ser494, Gly496, Thr500 and Asn501, π-σ bond with residues Tyr505, and one π-anion interaction with residues Glu484 and 9 Van der Waals interactions were found as well. Pibrentasvir, on the other hand, formed 9 conventional and carbon hydrogen bonds with residues Arg403, Glu484, Gln493, Ser494, THR500 and Gly502, 3 halogen bonds with residues Ser494, Gly496 and Asn501, 2 π-π stacked bond with residues Tyr505, one π-anion interaction with residues Glu406 and 7 Van der Waals interactions. Both zafirlukast and pranlukast, both being leukotriene receptor antagonists, formed a total of 18 non-covalent bonds. Zafirlukast formed 5 hydrogen bonds with residues Tyr449, Leu492, Gln493 and Ser494, 5 alkyl/π-alkyl bonds with residues Leu452, Tyr453, Tyr489 and Tyr495 one π-π T-shaped bond with residues Tyr505 and 7 Van der Waals interactions were found. Pranlukast formed the same number of hydrogen bonds as zafirlukast with residues Tyr449, Ser494, Gly496 and Tyr505, one π-σ bond with residues Phe490, one π-π T-shaped bond with residues Tyr449, one π-alkyl bond with residues Arg403 and many Van der Waals interactions were found as well. On the other hand, hydroxychloroquine had just 2 conventional hydrogen bonds with residues Gly485 and Gln493, 4 weaker carbon hydrogen bonds with residues Glu484 and Cys488, 3 alkyl/π-alkyl bonds with residues Leu452, Leu455 and Tyr489 and 6 Van der Waals interactions. A similar trend is observed in azithromycin, as it formed 2 conventional hydrogen bonds with residues Tyr489 and Gln493, 1 carbon hydrogen bond with residues Glu484, 5 π-alkyl bonds with residues Phe486, Tyr489 and Phe456 and 4 Van der Waals interactions. Azithromycin also had an unfavorable acceptor-acceptor interaction with Gly485. The total number of non-covalent interactions for these two drugs were found to be 15 and 13, which were less than all of our top 10 drugs.

Heatmap showing types and number of non-covalent interactions between SARS-CoV-2 S-protein residue and drugs. Since a residue can have more than one interaction, and residues critical for hACE2 protein are only shown in Figure 2(C), therefore, the number of interactions per drug can be more than the number of residues shown in Figure 2(C).

To evaluate the conformational behavior, stability and flexibility of the top drug molecules to bind S-protein, we subjected the structures of S protein-drug complexes to MD simulation. Root mean square deviation (RMSD) during the 50 ns simulation was calculated considering the proteins backbone with respect to the initial conformations (Figure 7). Five out of six drugs achieved a steady state after 5 ns while lonafarnib reached a steady-state only after 20 ns. The binding has also been visualized through movies (Supplementary movies S1 to S7). One of the top drugs, Tirilazad, showed unstable binding (). Tirilazad binding was unstable till 36 ns, although stable afterward, but it can again change conformation if MD is extended beyond 50 ns. This instability in the case of Tirilazad can be possibly due to the presence of two strong unfavorable acceptor-acceptor interactions with LEU-492 and TYR-505, causing it to change the docking site (as shown in Movie S7) although the initial docking model showed good score. The free binding energy of all SARS-CoV-2-RBD/drug complexes was calculated for the whole trajectories, and the results of energy components of the complexes are given in Figure 8. As expected from the RMSD plot, zafirlukast showed the highest affinity to SARS-CoV-2 S protein with a binding energy of −20.680 kcal/mol, followed by Velpatasvir and Pibrentasvir with the binding energy of −19.924 kcal/mol and −18.822 kcal/mol respectively. A slightly lower value of lonafarnib (binding energy = −13.576 kcal/mol) can be due to a partial shift in position (as shown in Movie S3), causing a change in RMSD value at t = 12 s. This could have decreased the overall binding energy even though the drug showed good docking score in the primary screening.

Figure 7.

RMSD plots of the SARS-CoV-2 S protein bound to pibrentasvir, velpatasvir, lonafarnib, zafirlukast, pranlukast, and ergocalciferol during 50 ns of simulations. The x-axis indicates time in ns and the y-axis indicates the RMSD value in nm.

Figure 8.

Binding energy components (kcal/mol) from MM/PBSA for pibrentasvir, velpatasvir, lonafarnib, zafirlukast, pranlukast, and ergocalciferol complexed with S-protein of SARS-CoV-2. The values show mean ± 1 standard deviation calculated through the python script MmPbSaStat.py provided in the g_mmpbsa package.

Discussion

We utilized the crystal structure of virus-receptor complex and a database of drugs already approved for human use to perform a computational screening for potential and effective drugs that could be repurposed to treat COVID-19. Our work showed that at least 10 better drug molecules are potentially more effective than hydroxychloroquine (HCQ) and azithromycin for inhibiting S-protein-hACE2 interaction. As the interaction between S-protein and hACE2 is crucial for SARS-CoV-2 transmission, we are confident that targeting this interaction can effectively reduce COVID-19 transmission and can bring down the disease severity. Potential drugs include velpatasvir, pibrentasvir, zafirlukast, pranlukast, lonafarnib and ergocalciferol. These drugs were not only able to efficiently bind to RBD of S-protein but were also able to competitively block the residues vital for hACE2 and S protein interaction. Thus, these molecules present strong, attractive candidates for the treatment of COVID-19.

A study by Zhang et al. showed that there were various types of SARS-CoV-2, but their phylogenetic relationships revealed two major genotypes, namely, Type I and II; Type I strains can be further divided into Type IA and IB (Zhang et al., 2020). The genomes of the two types mainly differ at three sites, which are positions 8750, 28112 and 29063, using the genome coordinates of sample MN938384.1 as a reference (Chan et al., 2020). However, none of these sites lie in the coding region of S protein (between nucleotides 21531 and 25352); hence a drug inhibiting the binding of S-protein and hACE2 receptor will be potent for all types of SARS-CoV-2.

One caveat of our work is that our results are entirely based on computational screening, and theoretical binding energy sometimes does not correlate well with IC50 (Cournia et al., 2017); hence we need to verify experimentally in the laboratory before human use. However, we believe this computational screening will help us find a quick and effective antidote against COVID-19. Besides, these results will accelerate the design of a tailored drug molecule against COVID-19 in the long run.