Drug repurposing for identification of potential inhibitors against SARS-CoV-2 spike receptor-binding domain: An in silico approach

SARS-CoV-2 belongs to one of the largest RNA virus genomes, ranging in size from 27 to 32 kb. The viral genome of SARS-CoV-2 contains non-structural protein genes with open-reading frame (ORF) 1a, ORF1b and four key structural proteins that are encoded by spike (S), envelope (E), membrane (M) and nucleocapsid (N) genes1. Among all the structural proteins, CoV spike (S) glycoprotein has been reported to play the most crucial role in viral attachment, fusion and entry and also used as a key target for the production of antibodies, vaccines and entry inhibitors. S glycoprotein is considered as the most important potential therapeutic target which is recognized by the cellular receptor and primed by the host cellular protease2. S1 and S2 are the two subunits of S protein; the S1 facilitates the entry of virus into the host cells by its binding to the host receptor through its receptor-binding domain (RBD), whereas S2 subunit facilitates the fusion of the viral and host membranes. The RBD in SARS-CoV-2 S protein was reported to be the gateway of infection that bound robustly to human angiotensin-converting enzyme 2 (ACE2) receptor3.

The SARS-CoV-2 genome-specific vaccines and their therapeutic antibodies are currently being tested. Alternatively, drug repurposing approaches of the existing therapeutic agents were carried out for clinical studies of COVID-19 therapeutics. Drugs such as remdesivir, favipiravir, hydroxychloroquine, ivermectin and lopinavir/ritonavir were repurposed for the treatment of COVID-19, based on their previous clinical history for potential therapeutics of other virus infections and pathologies4. These clinical therapies can be divided into two categories depending on their target: first, acting directly against SARS-CoV-2 either by inhibiting pivotal viral enzymes/proteins responsible for the replication of genomes or by preventing the entry of viruses into human cells; second, boosting the innate response by modulating the human immune system, this plays a key role against viruses, or by inhibiting the inflammatory processes that cause lung injury.

Based on the drug targets, clinical trials have been done on several classes of drugs including favipiravir and remdesivir (RNA polymerase inhibitors), lopinavir/ritonavir (protease inhibitors), chloroquine along with its hydroxyl derivative (aminoquinolines), xiyanping injection and corticosteroids (anti-inflammatory agents) and ACE2 inhibitors5. In addition to the clinical benefits of aminoquinolines, protease inhibitors and RNA polymerase inhibitors over COVID-19 therapeutics, there have been a few controversial findings reported from various research findings67.

Remdesivir was reported to inhibit SARS-CoV-2 which could be a prospective treatment of SARS-CoV-2, responsible for COVID-198. It was permitted to enter the clinical trials immediately under COVID-19 emergency conditions, on the basis of its safety and antiviral activities9. In an in vitro investigation, a hepatitis C virus protease inhibitor, namely simeprevir, was found to be a favourable repurposable drug for the treatment of COVID-19, which was also shown to synergize remdesivir in suppressing the replication of SARS-CoV-210.

Understanding the importance of RBD for COVID-19 therapeutics, limitations of vaccine development at a short interval, research controversies and prioritizing the urgent need of COVID-19 therapeutics, we performed an integrative computational approach to inhibit the RBD of S protein from SARS-CoV-2 using 56 drugs which were commercially available and used for therapeutics of various diseases. In this approach, an attempt was made to repurpose the drugs with potential antiviral activity against COVID-19 which could inhibit the RBD of S protein and ultimately prevent its entry into the human cells through ACE2 receptor.

Material & Methods

Sequence, structure and domain architecture analysis: The sequence, structure and functional information of RBD of S protein of SARS-CoV-2 were retrieved from National Center for Biotechnology Information (NCBI), GenBank with accession number MN908947. The spike (S) protein consisted of 1273 amino acids (aa), of which 229 were aa codes for RBD protein of SARS-CoV-2 that directly interacted with human ACE2. The RBD lies within the S1 region of a coronavirus S protein that triggers host-cell receptor binding activity11. The experimental structure of RBD with Protein Data Bank (PDB) ID: 6M0J was retrieved from Research Collaboratory for Structural Bioinformatics (RCSB), PDB with a resolution of 2.45 Å. The co-crystallized human ACE2 was removed using BIOVIA Discovery Studio 4.5 Visualizer (BIOVIA, San Diego, CA, USA).

Retrieval of drugs: The information of 56 commercially available drugs that are mostly prescribed for viral diseases, was retrieved from various sources6121314. The three-dimensional (3D) structure of all these 56 drugs was retrieved from the NCBI PubChem database in Structure Data Format, which was later converted to PDB format using online web server namely Simplified Molecular Input Line Entry System as per the requirement of AutoDock 4.2 docking tool (http://autodock.scripps.edu/).

Prediction of binding site: The key interacting residues were analyzed from the experimental structure of RBD-ACE2 complex (PDB ID: 6M0J) and available literature, which were identified as active site residues that take part in the binding site formation.

Molecular docking using AutoDock 4.2 tool: Molecular docking studies were carried out using experimental structure of RBD S protein against all the 56 drugs using AutoDock 4.2 tool. To allocate Kollman charges for the protein and Gasteiger partial charges for all the inhibitors, ADT v.1.5 was used. The grid with dimension space and parameters were as follows: x-centring: -37.872, y-centring: 28.878 and z-centring: 2.979, were chosen to allow full-extended conformation of the ligand. Based on the binding energy values, ligand efficiency, intermolecular hydrogen (H)-bonds and other hydrophobic and electrostatic interactions, the best resulting docked complexes were identified and processed for further computational analysis. The existence of intermolecular interactions between RBD–drug complexes was depicted through LigPlot+ (https://www.ebi.ac.uk/thornton-srv/software/LigPlus/).

Quantum chemical calculation using density functional theory (DFT): A density functional theory (DFT)-based study using the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) energy was carried out to investigate the reactivity and efficiency of potential drugs with antiviral activity against RBD by concerning the Becke, 3-parameter, LeeYang-Parr (B3LYP) correlation function of DFT15. The energy calculations were made using ORCA Program version 4.016. The electronic energy, frontier HOMOs, LUMOs, gap energy and dipole moment were calculated for the potential drugs. The DFT was calculated using the following equation:

E = min_n {∫ V_nuclei (r̄)n (r̄)d³r̄ + F[n (r̄]}

Where, V_nuclei+n r }n ≡ trial density and F ≡ universal functional.

Molecular dynamics (MD) simulations: GROMOS96 43A1 force field in the GROMACS version 5.1.4 molecular dynamics (MD) simulation package was used to analyze the Apo (protein only) and Holo (protein-ligand) states to understand the dynamic behaviour, binding mode and specificity of RBD inhibitors. RBD-simeprevir complex, which represented the highest binding affinity from our docking analysis, was further processed for MD simulations to explore its inhibitor specificity, dynamic behaviour and mode of binding activity. The 'pdb2gmx' programme of GROMACS package was used for the generation of topology file. The steepest descent method with a 1000 kJ/mol tolerance was used for energy of minimization as required for releasing the conflicting interactions. In the primary phase, constant number of particles, volume and temperature (NVT) ensemble was carried for temperature equilibration by restraining the positions of backbone atoms for 1000 picoseconds (psec) followed by the secondary phase where pressure equilibration was used in constant number of particles, pressure and temperature (NPT) ensemble for 1000 psec. A run of 30 nanoseconds (nsec) MD time period was set for both Apo and Holo states using periodic boundary conditions with constant temperature. Visual MD (VMD 1.9.1) was used to analyze the resultant trajectories which are inbuilt in GROMACS. The root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), total energy and solvent accessible surface area (SASA) were analyzed using gmx_rmsd, gmx_rmsf, gmx_gyrate, gmx_tenergy and gmx_sasa. All 2D plots were graphed using GRaphing, Advanced, Computation and Exploration 5.1.23 version (https://www.its.hku.hk/services/research/hpc/software/grace) for data analysis of MD simulations.

Principal component analysis (PCA): To recognize the coordinated movements in Apo and Holo states of RBD, principal component analysis (PCA) was performed through essential dynamics (ED) process using gmx_covera and gmx_aneig tools as per the protocol within the software package of GROMACS. After diagonalizing and calculating the covariance matrix representing the molecules' concerted motion, a set of eigenvectors and eigenvalues were obtained.

Results

Annotation of binding site of RBD of spike protein: The literature survey and crystal structure analysis of RBD–ACE2 complex denoted the residues Lys417, Gly446, Tyr449, Tyr453, Leu455, Ile472, Phe486, Asn487, Tyr489, Gln493, Gly 496, Thr500, Asn501, Gly502 and Tyr505 of RBD, which interact with human ACE2 during cell entry, thus annotated as active site residues that take part in binding site formation17.

Molecular docking: The binding free energies of all the 56 (RDB-drug complexes) interactions are summarized in Table I. The docking results reflected different binding free energies for different drug-target complexes, ranging from −0.23 to −8.52 kcal/mol. RBD-chloroquine phosphate docking complex represented the least binding energy, whereas RBD-simeprevir complex represented the highest binding affinity among the 56 drugs under investigation. The drugs remdesivir and favipiravir represented the binding energy −4.68 and −4.32 kcal/mol which were lower than simeprevir. To identify the potential drugs with better antiviral activity for SARS-CoV-2 therapeutics, a cut-off binding energy with a range from −5.98 to −8.52 kcal/mol was considered for further computational analysis. The selected range provided eight drugs consisting of chalcone, grazoprevir, enzaplatovir, dolutegravir, daclatasvir, tideglusib, presatovir and simeprevir. The drugs bortezomib, entecavir, ribavirin and trifluridine were considered due to the highest number of conventional H-bond interactions with RBD. Remdesivir was also considered under priority research area for COVID-19 therapeutics. Similarly, the drug favipiravir was considered for further analysis based on its applications for COVID-19 therapeutics.

Hydrogen bond analysis and inter-atomic distance calculation: The study showed that all docked complexes exhibited variable numbers of intermolecular H-bonding patterns. The docking analysis depicted eight H-bonds (average of ~2.179 and ~2.699 Å) in RBD complex with entecavir and ribavirin, which were the highest number among all the complexes (Fig. 1A and B). Remdesivir and favipiravir represented four H-bonds (average of ~2.907 and ~2.659 Å) (Fig. 1C and D). H-bond network of RBD-simeprevir complex resulted in four number of H-bonds in spite of highest binding affinity among all drugs. The H-bond network was between Lys417, Gln493 and Ser494 of RBD and simeprevir with an average distance of ~2.810 Å (Fig. 2A). The residues Tyr453, Leu455 and Lys417 were involved in hydrophobic interactions. The post-MD docking studies of RBD-simeprevir complex represented a binding affinity of −8.74 kcal/mol. The H-bond network was between Gly496 and Ser494 of RBD and simeprevir with an average distance of ~3.12 Å.

Close

Fig. 1

Intermolecular hydrogen bonding, electrostatic and hydrophobic interactions formed between (A) RBD-entecavir complex, (B) RBD-ribavirin complex, (C) RBD-remdesivir complex, (D) RBD-favipiravir complex. The images are drawn by LigPlot+ tool. RBD, receptor-binding domain.

Export to PPT

Close

Fig. 2

(A) Intermolecular hydrogen bonding, electrostatic and hydrophobic contacts formed between RBD-simeprevir complex drawn by and LigPlot+ tool. (B) LUMO and HOMO plots of simeprevir which exhibited higher reactivity against RBD. The positive electron density is indicated by red colour while blue colour indicates negative electron density. HOMO, highest occupied molecular orbital; LUMO, lowest unoccupied molecular orbital.

Export to PPT

Quantum chemical calculation: Owing to the importance of quantum computation, quantum chemistry was employed to study the frontier molecular descriptors such as HOMO and LUMO, gap energy and dipole moment for all the 14 shortlisted drugs that were predicted with better potential activity against RBD of S protein in SARS-CoV-2 (Table II). The effective reactivity for all the 14 drugs, which showed band energy gap (ΔE), i.e. the difference between E_LUMO and E_HOMO, ranged from 9.254 to 13.126 kcal/mol. Simeprevir displayed greatest reactivity against RBD among all the screened drugs based on its lowest band energy gap, which was calculated to be 9.254 kcal/mol (Fig. 2B). Taking together the results of molecular docking and DFT analysis, the drug simeprevir was further postulated for MD simulation along with RBD.

Trajectory analysis of MD simulations: The MD simulation of the Apo and Holo (RBD-simeprevir complex) states of RBD was carried out to evaluate the dynamics and stability of RBD protein, RMSD, Cα-RMSF, Rg, total energy and SASA from the trajectories resulted from MD simulations using GROMACS tools. The RMSD depiction for backbone residues was developed and plotted against a time scale of 30 nsec to access the dynamic stability of RBD. The backbone RMSD (Holo) was observed with a stable deviation after 20 nsec of simulation (Fig. 3A) when compared to its Apo state. The Apo state depicted a significant deviation from 0 to 30 nsec (0.1-0.25 nm) compared to the Holo state, with a stable RMSD with a value ranging from ~0.23 to ~0.22 nm from 25 to 30 nsec. This means that the drug simeprevir can help stabilize the protein by changing the conformation. The result of RMSD was further validated through fluctuation of residues using RMSF. The mobility of different residues was observed in RBD through RMSF plots (Fig. 3B). In Apo state, it was observed that the amino acid residues between 370-379 and 430-435 exhibited greater deviations in their Cα atoms in comparison to other regions. Around 10 residues (475-485) displayed greater deviations in Holo state of RBD as compared to its Apo state. This plot signifies that binding of simeprevir decreases the mobility of residues in Holo state than in Apo state. Rg was calculated to analyze the overall compactness for both the states. Gyration radius versus time graphs were plotted to check the compactness (Fig. 3C). The Apo state's Rg ranged from ~1.84 to ~1.85 nm, whereas the Holo Rg ranged from ~1.82 to ~1.79 nm, which represented higher Rg value in Apo state than in Holo state. The energy plot depicted the decreased mobility of residues in Holo state compared to Apo state, which was confirmed form RMSF plot (Fig. 3D). The hydrophobic interactions mediate the exposure of amino acids to certain solvent. The frequency of these interactions with the solvent and the core protein residues is directly proportional to the exposed surface area. The sketch of SASA (Fig. 3E) showed a reduction in the accessible solvent surface in the Holo state of RBD relative to its Apo state. SASA's findings showed the alteration of hydrophilic and hydrophobic interaction areas resulted by the binding of simeprevir to RBD, which could potentially prevent the host-viral interactions and ultimately making binding surface unavailable for the virus with human counterparts. Throughout the simulation time, the SASA graphs of the Holo state represented SASA with ~91 to ~110 nm², which was lower than the Apo state with a value of ~115 to ~ 120 nm² (Fig. 3E).

Close

Fig. 3

Conformational stability of RBD (Apo and Holo states) from SARS-CoV-2 spike protein throughout 30 nanoseconds (nsec) time period of MD simulations. (A) Backbone-RMSD of RBD. (B) Cα-RMSF profile of RBD. (C) Radius of gyration (Rg) profile of RBD. (D) Total energy of RBD and RBD-simeprevir complex (Apo and Holo state) during 30 nsec MD simulations. (E) Solvent accessible surface (SASA) analysis of RBD-simeprevir complex during 30 nsec MD simulations. The Apo and Holo are displayed by black and red lines, respectively. RMSD, root mean square deviation; RMSF, root mean square fluctuation; MD, molecular dynamics.

Export to PPT

H-bond analysis: The intermolecular hydrogen bonds of the RBD-simeprevir complex were tracked using the gmx_hbond tool of GROMACS (Fig. 4A). The simulation of Holo state represented an inconsistent number of intermolecular hydrogen bonds throughout the simulation time period. It represented four H-bonds (with an average atomic distance of ~1.67 nm). The number of H-bonds was directly proportional to the stability of the drug-target complex over the entire simulation time period. During simulations, a few crucial H-bonds such as Gln493 and Lys417 were broken, but later novel H-bonds, van der Waals and hydrophobic contacts were compensated (Fig. 4B). In spite of certain novel interactions with residues Leu452, Phe 456, Tyr489, Leu492 and Tyr505, it did not compensate with residue Ser494 (H-bond). This reflects its potentiality against the targeted protein as Ser494 which is one of the crucial residues in boosting the ACE2 binding18.

Close

Fig. 4

(A) Deviation of H-bonds contributed in interaction during 30 nsec simulation in RBD-simeprevir complex. (B) Post-MD simulations intermolecular hydrogen bonding, electrostatic and hydrophobic contacts formed between RBD-simeprevir complex drawn by Diglot+ tool. (C) The cloud represents the projection of trajectories eigenvectors (EV1 and EV2) (Black: Apo; Red: Holo). (D) Projection of the motion of the Apo and Holo states of RBD in phase space along the first two principal eigenvectors (EV1 and EV2).

Export to PPT

Principal component analysis (PCA): The movement of the RBD Apo and Holo states in phase space was captured by the trajectory projections from PC1 and PC2 (Fig. 4C), which were well aligned with RMSF (Fig. 3B). The trace values shown for Apo and Holo were 13.65 and 12.31 nm² for RBD. The lower Holo state trace values confirmed the overall decrease in RBD flexibility relative to its Apo state. The motion direction was shown by the vectorial representation of the solitary components. The majority of internal movements are shown by the initial vectors, while EV1 and EV2 represent a large number of overall movements. Following the plotting of eigenvalues against eigenvectors, steep curves of eigenvalues (Fig. 4D) were obtained.

Discussion

Attempts for the development of vaccines and direct-acting potential antiviral drugs are being carried out for effective treatment of COVID-19 therapeutics. Existing reports from various research findings have suggested that drugs such as chloroquine, hydroxychloroquine, arbidol, remdesivir, favipiravir, azithromycin and nelfinavir have shown efficacy and safety for COVID-19 treatment121314. Pharmaceutical and medicinal experts, however, raised the queries on their efficacy because both these drugs were originally designed to treat other diseases. This depicted their implementation through drug-repurposing. Considering the pros and cons of the COVID-19 investigational drugs, the current investigation was based on in silico screening of 56 commercially available drugs against RBD, through a computational approach to find out potential candidates for COVID-19 therapy. This investigation revealed chloroquine phosphate and simeprevir to have the least and highest antiviral activity based on their binding energy. RBD-simeprevir complex, which was approached for MD simulation predicted higher fluctuation pattern in a few residues of RBD, which might be due to their direct interaction with human ACE2 as reported by Wu et al1. The integrative computational approach of docking, quantum chemical calculation and MDS was used, which was able to detect the linkage between residual movements in RBD and its way of interactions with these drugs. The capability of simeprevir for COVID-19 therapeutics was eventually hypothesized by the Apo and Holo state trajectories. The overall findings obtained from various computational tools affirmed the pivotal role of the 14 drugs. This investigation was an in silico approach based on the information of drugs and experimentally derived crystal structure of RBD. The findings of this in silico investigation could be a supporting evidence for in vivo and in vitro studies needed to be carried out to confirm the efficacy and antiviral drug potency of simeprevir against RBD SARS-CoV-2 S protein.