Applying polypharmacology approach for drug repurposing for SARS-CoV2.

Exploring the new therapeutic indications of known drugs for treating COVID-19, popularly known as drug repurposing, is emerging as a pragmatic approach especially owing to the mounting pressure to control the pandemic. Targeting multiple targets with a single drug by employing drug repurposing known as the polypharmacology approach may be an optimised strategy for the development of effective therapeutics. In this study, virtual screening has been carried out on seven popular SARS-CoV-2 targets (3CLpro, PLpro, RdRp (NSP12), NSP13, NSP14, NSP15, and NSP16). A total of 4015 approved drugs were screened against these targets. Four drugs namely venetoclax, tirilazad, acetyldigitoxin, and ledipasvir have been selected based on the docking score, ability to interact with four or more targets and having a reasonably good number of interactions with key residues in the targets. The MD simulations and MM-PBSA studies showed reasonable stability of protein-drug complexes and sustainability of key interactions between the drugs with their respective targets throughout the course of MD simulations. The identified four drug molecules were also compared with the known drugs namely elbasvir and nafamostat. While the study has provided a detailed account of the chosen protein-drug complexes, it has explored the nature of seven important targets of SARS-CoV-2 by evaluating the protein-drug complexation process in great detail. Graphical abstract: Drug repurposing strategy against SARS-CoV2 drug targets. Computational analysis was performed to identify repurposable approved drug candidates against SARS-CoV2 using approaches such as virtual screening, molecular dynamics simulation and MM-PBSA calculations. Four drugs namely venetoclax, tirilazad, acetyldigitoxin, and ledipasvir have been selected as potential candidates. Supplementary Information: The online version contains supplementary material available at 10.1007/s12039-022-02046-0. © Indian Academy of Sciences 2022.

Introduction

The severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2) resulted in a pandemic[1-3] and continues to cause turmoil in our everyday lives.[4] Tackling the challenges caused by viral infections often involves the development of therapeutics which not only control the pristine form of the virus but also its variants.[5] The SARS-CoV-2 pandemic triggered a collapse of the healthcare system and has instigated the search for therapeutics. Though there are few therapies available for the treatment of SARS-CoV-2 infection, the effective cure is still a matter of question.[6-9] The major factor of this enhanced transmission may be traced to the emergence of new mutant strains.[10,11] Around 30 SARS-CoV-2 proteins have been reported so far[12] and the manifestation of mutations in terms of infectivity, transmission, replication, and transcription is an intriguing fact.[13-15] The genome sequence of SARS-CoV-2 is 96.2% identical to the bat CoV RaTG13 sequence, although it shares 79.5% identity to SAR-CoV.[16]

SARS-CoV-2 encodes four structural proteins, and sixteen non-structural proteins (NSPs). The four structural proteins make up a mature SARS-CoV-2 viz, envelope (E), membrane (M), nucleocapsid (N) and spike (S).[17] The Spike (S) homotrimeric glycoprotein located on the virion surface allows entry of the virus into the human cells and employs the cellular transmembrane serine protease 2 (TMPRSS2) for S protein priming which helps in internalization of the virus.[18,19]

Figure 1 displays the viral targets selected in the current study based on their role in viral replication and RNA synthesis. The 3 cysteine-like protease (3CLpro)[20] and papain-like protease (PLpro)[21] are involved in the cleaving of polyprotein to release independent functional proteins, where the former, due to its high sequence similarity among the coronavirus, may be studied for developing a general-purpose antiviral drug. Further RNA-dependent RNA polymerase (RdRp (NSP12)),[22] Helicase (NSP13),[23] Exoribonuclease (NSP14)[24] and Endoribonuclease (NSP15)[25] play a pivotal role in SARS-CoV-2 replication and RNA synthesis, while 2'-O-methyltransferase (NSP16) stabilizes the RNA, and protect the viral RNA from host immune response.[26] Thus, the current study has a robust representation of proteins involved in RNA synthesis, cleaving of polyprotein to release functional proteins, and replication. Spike proteins, which are essentially responsible for the viral entry was not included as they have been extensively studied, besides the current study focuses on selecting the proteins involved in viral transcription and translation. Thus, the chosen proteins, as depicted in Figure 1 are capable of providing common features for developing multi-targeting drugs against COVID-19.

Figure 1

Cartoon representation of the selected protein targets (A) 3CLpro PDB ID:6M03 (B) PLpro PDB ID: 6W9C (C) NSP12: Modelled (D) NSP13: Modelled, (E) NSP14: Modelled, (F) NSP15 PDB ID: 6W01, (G) NSP16 PDB ID: 6W4H, along with labelled domains of each protein. The selected binding site residues are encircled and highlighted which depict the amino acid residues of the binding site residues.

The 3CLpro protein consists of three domains (domain I, domain II, domain III) where domain I and domain II have a five-stranded antiparallel β-barrel structure with a chymotrypsin-like folding scaffold, whereas the domain III has a five α-helices cluster.[27,28] PLpro protein has two domains which are ubiquitin-like (UBL) domain in N-terminal and it is necessary for pathway antagonism and ubiquitin-specific protease (USP) domain in C-terminal which contains the protease active sites.[29,30] The structure of NSP12 protein contains RdRp domain, RdRp associate nucleotidyltransferase (NiRAN) domain and an interface domain connecting the two domains and these domains, consist of the active sites and thus, play a key functional role of NSP12 protein.[31] There are five domains for NSP13 which consist of zinc-binding domain (ZBD) which is involved in coordinating three structural zinc ions, a helical stalk domain, a beta-barrel 1B domain and two RecA like helicase subdomains which are 1A and 2A. This RecA like” helicase subdomains contain residues that are involved in binding of nucleotides and also hydrolysis.[32] The NSP14 protein structure contains two domains which are exoribonuclease (ExoN), in N-terminal and N7‐methyltransferase (N7‐MTase) domain in C-terminal. The replication proofreading activity of Exon domain makes it a potential target for developing inhibitors against SARS-CoV-2.[33] NSP15 consists of three domains which are N-terminal domain (ND), middle domain (MD), and C-terminal catalytic NendoU domain (NendoU). The EndoU domains are involved in biological functions which are associated with RNA processing and also in innate immune response.[34] Finally, the NSP16 consist of a methyltransferase (MTase) domain and the alpha-D helix in MTase domain is important for SAM-binding (S-adenosyl-L-methionine-dependent methyltransferase) and RNA cap-binding for methyltransferase activity.[35]

Vaccine and drug development for the prevention and treatment of COVID-19 has become an area of interest for the scientific community worldwide. Several vaccines have been developed globally, but fail to eradicate the virus completely due to the emergence of new viral strains. Thus, employing drug repurposing[36] with a polypharmacology[37] approach may decipher the process to tackle the therapeutic emergency. As a response to the COVID-19 pandemic, several studies have utilized various molecular modeling approaches for repurposing pre-existing drugs.[38-41] There are a few repurposed drugs available in the market for COVID-19 treatment, however, several limitations have been observed with the existing drugs.[9,42-44] To examine the relevance of the selected seven targets, a docking study has been performed on the 15 clinically repurposed drugs that were extensively employed for the COVID-19 treatment (Figure 2 and Figure S1, SI). Figure 2 provides an Autodock Vina analysis of the 105 complexes (all complexes of 15 clinically used drugs with 7 targets) and 83 of these combinations have shown fair binding affinity.

Figure 2

The distribution of docking scores of known popular drugs with seven SARS-CoV-2 drug targets. The drugs are shown in X-axis and the proteins are shown on Y-axis. The yellow-blue color spectrum shows increasing binding affinity. The drug molecules which have a docking score of less than −6 kcal/mol are represented with white space in corresponding proteins.

This supports the selection of drug targets having a good ability to identify the promiscuous candidates in the docking study. Extensive docking analysis was carried out on the known drug candidates, as provided in the drug bank and drug central using Autodock Vina[45] to see the interactions and binding affinity of 4015 approved drugs against well-known seven NSPs of SARS-CoV-2, namely 3CLpro, PLpro, RdRp, NSP13, NSP14, NSP15, and NSP16. Analysis of virtual screening results followed by 100 ns MD simulation and MM-PBSA calculations indicate that venetoclax, tirilazad, acetyldigitoxin and ledipasvir showed reasonable stability and retain consistent interaction with their respective multiple targets. The schematic workflow adopted in this study has been displayed in Figure 3.

Figure 3

Schematic representation of workflow adopted in this study to identify repurposable drug candidates against multiple SARS-CoV-2 targets.

Selection of target proteins and their structural details

A total of seven SARS-CoV-2 protein targets have been selected in this study based on their role in COVID-19 pathophysiology, especially viral RNA replication. The therapeutic targets were selected based on (i) the significance of particular proteins in viral replication and transcription, (ii) the presence of a binding pocket that presents hotspots for the drug development, and (iii) availability of viral protein with a good resolution from RCSB Protein Data Bank (PDB) database (https://www.rcsb.org).[46] After analyzing the information on bound inhibitor, missing residues and resolution of the structure, four crystal structures of SARS-CoV-2 were taken from PDB. Table S1 (SI) provides the PDB ID, active site residues of the chosen targets namely: (i) 3CLpro (ii) PLpro (iii) NSP15, (iv) NSP16 and three homology modeled proteins namely: (v) NSP12 (vi) NSP13 (Helicase), (vii) NSP14 were obtained from SWISSMODEL server[47] (https://swissmodel.expasy.org) due to unavailability of resolved crystal structures during this study.

The 3CLpro is a key enzyme of the virus, as it cleaves viral replicase polyprotein at 11 sites and generates NSP 4 to 16 that are essential for viral replication. The CLpro active site contains the catalytic dyad Cis145 and His41. 3CLpro of SARS-CoV-2 shares 96% similarity with that of SARS-CoV.[48] The PLpro (multi-domain protein) cleaves viral replicase polyprotein at three sites to release NSP 1 to NSP3 which antagonize the host innate immunity. The PLpro sequence identity between SARS-CoV-2 and SARS-CoVis 83%.[49]

The RdRp enzyme (NSP12) is a core catalytic subunit essential for viral RNA synthesis replication. RdRp proteins have been targeted for the development of potential antiviral molecules against different viruses including hepatitis C virus (HCV),[50] ZIKA virus (ZIKV),[51] and coronaviruses (CoVs).[52] The nucleoside triphosphate (NTP) entry channel for RdRp of SARS-CoV-2 is formed by a stretch of amino acid Ser759 to Asp761 and other amino acid residues such as Lys545, Arg553 and Arg555. RdRp of SARS-CoV-2 is highly similar to other coronaviruses such as SARS-CoV with which it has 96% sequence identity.[53]

NSP13 belongs to the 1B helicase superfamily with NTPase and it is involved in duplex RNA/DNA unwinding and 5′-RNA capping activities.[22] The site-directed mutagenesis studies on SARS-CoV NSP13 have revealed that six residues, Lys288, Ser289, Asp374, Glu375, Gln404, and Arg567 are essential for NTPase activity.[54] Due to their essential role in viral replication and conservation across all coronavirus species, NSP13 has emerged as a potential target for antiviral drug development.[22,54]

It has been observed that NSP14 is key for viral mRNA stability as it evades the host immune response.[24,55] The residues Ile332, Phe367, Val389, Asp331, Gly333, Trp292, Phe401, Tyr420, Phe426, Thr428, and Phe506 are important for N7-MTase activity. In the present study, the SAM binding domain of NSP14 has been considered as a binding pocket for drug binding.[55] Processes such as capping of viral mRNA are catalyzed by NSP16, which shares above 90% similarity among the novel and older SARS-coronavirus.

Uridylate-specific endoribonuclease (NSP15) cleaves 3′ of uridylates via a ribonuclease A (RNase-A) and it restrains exposure to the host immune system.[25] It has very high sequence similarity with NSP15 protein of SARS and MERS and the catalytic C-terminal domain (active site) of the enzyme-containing the catalytic triad His235, His250 and Lys290.[34,56] Encouragingly due to the SAM binding site's conserved nature, inhibitors targeting this pocket may be developed as pan-antiviral inhibitors.[57] The binding pocket of the selected seven SARS-CoV-2 proteins active sites are shown in Figure 1 and Table S1, SI.

Computational

Homology modelling of NSP12, NSP13 and NSP14 protein structures

The modelled structure of NSP12 (YP_009724389.1), NSP13 (YP_009724389.1) and NSP14 (BAE93399.1) protein was downloaded from SWISS-Model repository (https://swissmodel.expasy.org/repository/).[58] The selected modelled structures were validated using Ramachandran (RC) plot[59] and ERRAT[60] validation tools available in the SAVES server (https://saves.mbi.ucla.edu/).

Preparation of protein structures

All the seven selected SARS-CoV-2 proteins were prepared in Autodock Tools (ADT version 1.5.7rc1) by removing solvent molecules and non-interacting ions, adding polar hydrogens and partial charge assignment. The protein structures were saved in Autodock PDBQT. Based on the information of active site residues and their adjacent residues, the grid box size was generated. Table S2 (SI) provides the grid box parameters for all the seven SARS-CoV-2 protein targets.

Preparation of drug structures

A total of 4015 known drug candidates retrieved from DrugBank[61] and DrugCentral[62] were considered for the present study. Redundant molecules and ion-containing molecules were removed manually during the SMILES notations inspection. Further, the geometry optimization was performed using the raccoon.py script of MGL tools and converted to 3D structures. The selected non-redundant molecules were further converted to PDBQT format.

Virtual screening

The virtual screening was performed using AutodockVina1.1.2. The grid box was set up around the active site of the protein (Table S2, SI). Total nine conformers were generated and the best conformer was identified based on docking score and the binding pose. The top-ranked 20 candidates for each target were identified from virtual screening based on the high docking score. The top 25 drugs were prioritized according to the calculated physicochemical properties (Table S3, SI) and pharmacokinetic properties (Table S4, SI). This analysis leads to the identification of 25 drug candidates, which bind to at least two or more of the selected seven targets, by imposing a threshold value of -6.0 kcal/mol.

Further, the drugs having good docking scores and interactions with four or more targets were selected and subjected to post-docking analysis. Discovery studio visualizer[63] was used to analyze the protein-drug complex interactions. The schematic representation of the workflow adopted in this study has been shown in Figure 3.

Molecular dynamics simulations and MM-PBSA calculations

As described in the foregoing sections, 18 protein-drug complexes have been identified and subjected to MD simulations. MD simulations were performed with Gromacs 5.0.4 package[64] using CHARMM force field[65] and SPC water model.[66] The topological parameters of the drug molecules were generated using the CGenFF server.[67] MD simulations were performed under periodic boundary conditions (PBCs) with a cubical box maintaining a distance of 1.0 nm and 1.5 nm and the boundary edges between the protein-drug complexes depending on the size of the proteins. The protein-drug complexes were solvated with SPC water molecules and the system was neutralized by adding the counterions into the solvated box depending on the charge of the system. Then, the systems were subjected to 1000 step steepest descent minimization followed by 2000 conjugate-gradient minimization for initial energy minimization with the frozen drug to avoid a further structural clash in the solvated system. Then, the whole system was subjected to 5000 steepest descent minimizations followed by 6000 conjugate-gradient minimizations with a maximum step size of 0.01 fs for final energy minimization. After subjecting the systems to position restraint equilibrations, heat under canonical ensemble from 0 to 303 K for 500 ps was generated using a modified Berendsen thermostat.[68] Later, the systems were equilibrated for 1 ns under isothermal-isobaric conditions (with a constant pressure of 1.0 bar). Finally, the production run of 100 ns was performed without any restraints, followed by an integration time step of 2 fs. The coordinates were saved every 2 ps under constant conditions of 300 K temperature and 1 atm pressure.[69] The LINCS algorithm[70] was used to restrain the bond lengths and the long-range electrostatics were calculated using the particle mesh Ewald (PME),[71] while the SETTLE algorithm[72] was employed to constrain the geometry of water molecules. The trajectory of each system was analyzed using Origin Pro[73] and VMD software.[74] Molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method using g_mmpbsa package,[75] was employed to calculate the binding free energies of protein-drug complexes using the following equation.Where, Ecomplex, Eprotein and Edrug are the total MMPBSA energy of protein-drug complex, total solution free energies of the isolated protein and drug molecules, respectively. The total free energy of each individual can be expressed using equation 2.Where, z is either protein, drug or protein-drug complex. is the ensemble average value of the molecular mechanic’s potential energies in a vacuum. T and S represent absolute temperature and entropy, respectively and together TS contributes to the entropic contribution to the total free energy in a vacuum. is the sum of both bonded and nonbonded interactions of the molecules and is the average free energy of solvation, consisting of polar (Epolar) and nonpolar (Enonpolar) components.

The Poisson-Boltzmann equation is solved using a linear model to estimate the polar solvation energy (Epolar)whereas the nonpolar part of the solvation energy (Enonpolar) is estimated using solvent accessible surface area (SASA) method and can be expressed using equation 3.Where, γ is a parameter representing surface tension of the solvent, A is the calculated SASA of the molecule and b is a fitting parameter. The last 20 ns trajectory (i.e., 80 ns to 100 ns) was used for the MMPBSA calculation. MmPbSaStat.py program was used to calculate the binding energies and MmPbSaDecomp.py was used to extract the residue-specific contributions towards binding.[75] In this current investigation, we retrieved 400 frames from 80–100 ns trajectories from all the eighteen complexes, to calculate the binding free energy. To analyse the docking score distribution of all the four drugs with their respective targets, different conformers were generated from the 100 ns trajectories by taking a pose at interval of 10 ns. The drug molecules were removed and protein was prepared without changing their conformation. The putative 180 (18 * 10) points generated along the trajectory of MD simulations were subjected to docking.

Results and Discussion

In the current study, a virtual screening analysis of the known drug molecules obtained from Drug bank[61] and Drug central[62] was carried out against the seven SARS-CoV2 viral targets (3CLpro, PLpro, NSP12, NSP13, NSP14, NSP15 & NSP16). The top 20 compounds which amount to the list of top 0.4% of the molecules have been selected for further studies. It was observed that none of the known drug molecules which are popularly used in the treatment of COVID-19 were listed out in the top 0.4% (Figure 2) though most of them show a fair binding with multiple targets. Around 25 of the known drugs appear to have significant binding affinity for at least 2 or more targets as shown in Table 1. Among these, 4 drugs namely venetoclax, acetyldigitoxin, tirilazad and ledipasvir were considered for further studies. The two drugs ergotamine and dihydroergotamine were omitted from the study as they appear to be less promising candidates for drug repurposing. Thus, the 18 selected complexes were subjected to MD simulations and MM-PBSA calculations. In addition, an analysis of the docking poses, RMSF, RMSD, H-bond have been carried out along the 100 ns MD trajectory of these 18 complexes, which provided interesting insights into the nature of binding.

Table 1

List of 25 drug molecules and their therapeutic indications.

Sl. No.	Drugs	Known indication	Targets
1.	Dihydroergotamine	Migraine	3CLpro, PLpro, NSP13, NSP14, NSP15, NSP16
2.	Ergotamine	Migraine	3CLpro, PLpro, Nsp12, NSP13, NSP14, NSP15, NSP16
3.	Venetoclax	Anticancer	3CLpro, PLpro, NSP12, NSP13, NSP14, NSP16
4.	Acetyldigitoxin	Cardiovascular	3CLpro, NSP13, NSP14, NSP16
5.	Dihydroergocristine	Migraine	NSP12, NSP13, NSP15, NSP16
6.	Ledipasvir	Antiviral	PLpro, NSP12, NSP14, NSP16
7.	Tirilazad	Cardiovascular	3CLpro, PLpro, NSP12, NSP16
8.	Antrafenine	Antiinflammatory	PLpro, NSP14, NSP15
9.	Conivaptan	Hyponatremia	NSP13, NSP14, NSP15
10.	Entrectinib	Anticancer	3CLpro, PLpro, NSP16
11.	Paritaprevir	Antiviral	3CLpro, NSP13, NSP16
12.	Ubrogepant	Migraine	3CLpro, PLpro, NSP15
13.	Bisdequalinium	Antibacterial	NSP12, NSP13
14.	Digitoxin	Cardiovascular	NSP13, NSP14
15.	Eltrombopag	Anemia	NSP14, NSP16
16.	Glycyrrhizic acid	Antidiabetic	NSP13, NSP15
17.	Quinupristin	Antibacterial	NSP12, NSP16
18.	Teicoplaninaglycone	Antibacterial	3CLpro, NSP16
19.	Temoporfin	Anticancer	3CLpro, PLpro
20.	Zinostatin	Anticancer	NSP14, NSP15
21.	Guamecycline	Anticancer	3CLpro, NSP12
22.	Anazolene	Antiallergic	NSP15, NSP14
23.	Deslanoside	Cardiovascular	NSP12, NSP13
24.	Dexamethasone	Anticancer	NSP15, 3CLpro
25.	Lumacaftor	Cystic fibrosis	NSP15, PLpro

The last column lists the SARS-CoV-2 targets which show noticeable binding affinity with the drug as indicated by the docking score.

Homology modeling

From the models available at the SWISS-Model repository, model template PDB ID: 6M71 was chosen for NSP12 which showed 99% sequence identity with ERRAT value of 92%. Ramachandran plot of the modelled structure reveals that there is nearly 91% favored region, 7.7% allowed region and 0.7% outlier region. For NSP13, model template PDB ID: 5RLC was selected which showed 96% sequence identity and also ERRAT value of 96%. From the RC plot, it can be seen that the selected model showed 89.6% favored region, 9.8% allowed region and 0.2% outlier region. Similarly, for NSP14, model template PDB ID: 5NFY was selected which had 90% sequence identity with ERRAT value of 90% and the RC plot validates the selected model having 86.0% favored region, 11.6% allowed region and 1.7% outlier region. The template models were selected based on their sequence identity and query coverage which was followed by calculating the RC plot with ERRAT value that validated the model structure for further study. The RC plot and ERRAT values are shown in Table S5, SI.

Virtual screening analysis

Docking studies were performed for over 4015 approved drugs on the binding pocket of selected SARS-CoV-2 target proteins. The distribution of binding energies of 4015 approved drugs for the seven viral targets is shown in Figure 4. In this study, the top 20 drug molecules (Table 1) from each of the seven targets were identified based on the highest docking score obtained from AutodockVina. Venetoclax was found to have significant interaction with 6 targets while tirilazad, acetyldigitoxin and ledipasvir interacted significantly with 4 targets. Thus, the identified 18 complexes were subjected to MD simulations.

Figure 4

The distribution of docking scores of different drug molecules with seven selected SARS-CoV-2 drug targets. The drugs are shown in X-axis and the proteins are shown in Y-axis. The yellow-green color spectrum shows increasing binding affinity. The drug molecules which has less than − 6 kcal/mol or no docking score are represented with white space in corresponding proteins. The four selected drugs are highlighted in boldface.

Molecular dynamics simulations

Molecular dynamics (MD) simulations for the top 18 protein-drug complexes were carried out to study the conformational stability and sustainability of the protein-ligand complex. MD simulations can effectively validate the results obtained from molecular docking approaches, by examining the variation of the binding interactions in the complexes at various points along the trajectory. Similarly, the computed RMSD values will indicate the stability and rigidity of the complex (Figure 5). Figure 5A shows the RMSD for the seven targets in presence of venetoclax. NSP14 showed an increase in RMSD after 80 ns suggesting conformational flexibility. It was observed that NSP14 showed larger conformational flexibility when it is bound to venetoclax while in the case of other drugs such as acetyldigitoxin, and ledipasvir, the protein RMSD shows a convergence suggesting conformational stability. Significant changes in conformational nature were observed for PLpro bound to tirilazad, as for the remaining cases, the structural changes in proteins were observed to be less significant. In Figure S2 (SI), the apo protein backbone RMSD has been illustrated. It can be observed that the protein structure has been largely intact even after binding of the drug molecules, considered in the study. Figure 5B illustrate the variations in the RMSD values corresponding to the drugs when they are bound to different protein targets. Here the convergence reveal that the drugs are in the binding pocket of the target protein during the simulation time period. However, substantial fluctuations in the RMSD are an indication pointing to the possibility that the drug obtained an alternative conformation in the binding site. Depending upon the relative free energy differences, different conformations are sampled either less frequently or often. Even though the docking protocol gives relative energies of different docking poses, the discussion is usually based on the most stable binding pose. So, computing the free energies using configurations obtained by molecular dynamics simulations is necessary to obtain an accurate estimate for the binding free energies. It is worth noting that tirilazad in PLpro target displayed fluctuations in RMSD suggesting its conformational flexibility. However, ledipasivir exhibited less significant conformational flexibility in its targets. Since the binding site microenvironments are specific for each target, different conformations may be required for optimal binding with different targets. RMSFs are the measures of conformational flexibility of various residues and larger values for RMSFs suggest larger conformational changes while near-zero values suggest their rigid nature. RMSFs for various target proteins bound to different drugs are shown in Figure S2, SI. The protein bound to tirilazad showed the highest RMSF value suggesting that this induced a larger conformational change in the target protein. In NSP14, ledipasivir and venetoclax induce larger structural changes in protein in comparison to drug acetyldigitoxin. In the remaining targets, different drugs were observed to exhibit similar effects of the RMSFs (Figure S3, SI). The compactness in the overall protein can be compared from the radius of gyration which supports the RMSD results (Figure 5).

Figure 5

RMSD plot generated from 100 ns trajectory (A) Protein backbone plots are shown in the left panel (B) Ligand RMSD plots of four selected drugs with its respective proteins are shown in the right panel.

In order to validate our findings, we also performed MD simulation for the reported drug molecules elbasvir and nafamostat (interacting with all the targets of venetoclax). Comparison with the known reference molecules, elbasvir and nafamostat, reveals a high similarity to the shortlisted drug molecules in the study.

Hydrogen bond (H-bond) interaction is one of the major players in governing drug stability at the protein's active site. All the 18 protein-ligand complexes reveal the existence of multiple hydrogen bonds, with the number of such bonds fluctuating between two to five (Figure S4, SI). Analysing the frames along the molecular dynamics trajectory indicate that the number of hydrogen bonds were fluctuating between three to four, among the six complexes of the venetoclax. Similarly, the number of hydrogen bonds in the complexes involving tirilazad, acetyldigitoxin and ledipasvir ranges from two to three, three to four and two to five, respectively.

Analysis of non-covalent interactions of the protein-drug complex along the MD trajectory

The interaction pattern of the initial docked complex before MD simulation and after MD simulation of venetoclax, tirilazad, acetyldigitoxin and ledipasvir in complex with the respective protein targets were analysed. The interactions at 100 ns are depicted in Table S6 (SI). In the initial venetoclax-3CLpro complex, there was one H-bond, one π-σ bond, one C-H….O bond, two π-alkyl, one alkyl interactions as shown in Figure 6A and Table S6, SI. However, H-bond with residue Thr24, one π-σ bond with Thr25, one π-alkyl with Met49, one alkyl with Met165 were retained after 100ns MD simulations. In addition, one new H-bond with Thr26, two π-alkyl bonds with His41 and one alkyl with Cys145 were exhibited in this complex. Though, C-H bond with Ser46 was replaced by Asn142 during MD simulation.

Figure 6

2D representation depicting docking interactions of venetoclax with SARS-CoV-2 protein targets (A) 3CLpro (−9.1 kcal/mol) (B) PLpro (−8.2 kcal/mol) (C) NSP12 (−10.2 kcal/mol) (D) NSP13 (−10.8 kcal/mol) (E) NSP14 (−12.1 kcal/mol) and (F) NSP16 (−11.8 kcal/mol).

Venetoclax-PLpro complex displayed five H-bonds, one C-H….O bond, one π-π T-shaped bond, one π-alkyl and three alkyl interactions were found in the initial docked complex as shown in Figure 6B. On the other hand, none of the interactions was retained except one C-H….O interaction with Asp164 during MD simulation. However, some new interactions were established after MD simulations such as two H-bonds with Arg166, one alkyl interaction with Pro298, and one C-H interaction with Lys157 in venetoclax-PLpro complex (Table S6, SI). Venetoclax was not able to retain H-bonds, CH….N bond (His810), C-H….O (Lys798, Asp618) bond, alkyl, π-π T-shaped, π-anion and halogen interactions with NSP12 after MD simulations. However, venetoclax retained only two π-alkyl bonds with His810 and Lys798 as depicted in Figure 6C and Table S6 (SI). Venetoclax complex maintained only two π-alkyl interactions with Tyr180 and Pro408 residues of NSP13 during MD simulation and for venetoclax-NSP14 complex, none of the interactions was reserved which were existed in the initial docked complex. Furthermore, there were some new interactions were observed with NSP14 such as one H-bond with Gln262, one π-π stacked interaction with Lys288 and one π-alkyl interaction with Val290 as shown in Figure 6D and Table S6 (SI).

After MD simulation of venetoclax-NSP16 complex, only one H-bond with Asp6873 was retained whereas some new interactions formed such as π-π T-shaped bonds with residues Leu6898, Pro6932 and Val6937, and one π-anion interaction with Ser6872 residue of NSP16 (Table S6, SI).

From the initial docked complex of tirilazad-3CLpro, it can be observed that there is one H-bond, one cation-π, two π-alkyl and seven alkyl interactions as shown in Figure 7A and Table S7 (SI). However, after 100 ns of simulation, only one H bond with residue Arg188 was retained. Interaction of tirilazad with PLpro showed one H-bond, two C-H…O bond, one π-σ bond, two π-anion, one π-alkyl, and three alkyl interactions were observed as shown in Table S7, SI and Figure 7B. Furthermore, one π-anion with Asp164, one π-alkyl with Tyr264, two alkyl interactions with Leu162 and Pro248 were intact after 100 ns MD simulations (Table S7, SI). In addition to one C-H…O bond with Asp164, one new C-H….O bond with Pro248 was observed after MD simulation. However, in tirilazad-PLpro complex, tirilazad missed H-bond with Trp800, π-sulfur interaction with Cys813 and one π-alkyl bond with Phe812 of NSP12 after MD simulations (Figure 7C and Table S7, SI). Whereas, tirilazad retained two alkyl bonds out of six with Lys593 and Val588, one π-alkyl bond with Phe594 and C-H….O bond with Thr591 as depicted in Table S7, SI. A new alkyl bond was also formed with Leu758 of NSP12. In the case of initial tirilazad-NSP16 complex, there were two H-bonds, two C-H bond, three alkyl, one π-alkyl and one π-anion interactions present (Table S7, SI and Figure 7D). However, one H-bond with Lys6968, one C-H….O bond with Try6930, one alkyl interaction with Met6840, one π–alkyl interaction with Tyr6930 and one π–anion interaction with Asp6879 were retained after MD simulations. Additionally, some new interactions were observed such as one H-bond with Ser7000, two C-H bonds with Ser6999 and Asp6897, and two π–alkyl interactions with Met6929 and Leu6898.

Figure 7

2D representation depicting docking interactions of tirilazad with SARS-CoV-2 protein targets (A) 3CLpro (−8.7 kcal/mol) (B) PLpro (−7.9 kcal/mol) (C) NSP12 (−9.7 kcal/mol) and (D) NSP16 (−11.4 kcal/mol).

In the acetyldigitoxin-3CLpro complex, initially, there were seven H-bonds and five alkyl interactions as shown in Figure 8A and Table S8, SI. However, only four H-bonds with residue Asn142, Gly143, His163 and Glu166, and four alkyl interactions with Met49, Cys145, Leu27, Met165 were retained after 100ns MD simulations. In addition, one new H-bond with Thr26 was found in this complex. The acetyldigitoxin with NSP13 docked complex exhibited five H-bonds, two C-H bond and five alkyl interactions as shown in Figure 8B and Table S8, SI. On the other hand, none of the H-bond was retained but three new H-bonds with Thr380, Asp383 and Ser539 were found in post-MD interactions. One C-H….O bond with Ser535 and two alkyl bonds with Ala312 and Cys309 were reserved and two new alkyl interactions with Ala407 and Arg339 were also observed in post-MD complex of Acetyldigitoxin-Nsp13 as mentioned in Table S8, SI). In theacetyldigitoxin-NSP14 complex, (Figure 8C and Table S8, SI) one new H-bond with Leu144 and two new alkyl interactions with His378 and His401 were observed other than existing interactions initial complex. Only one alkyl interaction with Phe504 was retained whereas one new alkyl interactions with His345 of NSP14 was established after MD simulations (Table S8, SI). In acetyldigitoxin-NSP16 complex, initially, there were four H-bonds, one C-H….O bond, four alkyl, and two π-alkyl interactions were present (Table S8, SI and Figure 8D). Furthermore, one new H-bond with Thr6970, one C-H….O bond with Gly6946, and two π–alkyl interactions with Pro6932 and Leu6898 were established other than existed interactions in acetyldigitoxin-NSP16 complex after MD simulations (Table S8, SI).

Figure 8

2D representation depicting docking interactions of acetyldigitoxin with SARS-CoV-2 protein targets (A) 3CLpro (− 9.6 kcal/mol) (B) NSP13 (− 10.5 kcal/mol) (C) NSP14 (− 12.0 kcal/mol) and (D) NSP16 (− 11.3 kcal/mol). The docking score obtained from Autodock Vina for each protein are given in squared bracket.

Initially, there were two H-bonds, one C-H….O bond, one π-alkyl, two alkyl, one π-anion and one halogen interaction in the ledipasvir-PLpro complex, as shown in Table S9 (SI) and Figure 9A. On the other hand, after MD simulation, some new interactions were observed (one H-bond with residue Glu203, two alkyl interactions with Tyr264 and Leu162, and two C-H bonds with Thr301 and Asp164). However, one π-alkyl interaction with Tyr264, one π-anion and one halogen interaction with Glu167 was retained during MD simulation of the ledipasvir-PLpro complex.

Figure 9

2D representation depicting docking interactions of ledipasvir with SARS-CoV-2 protein targets (A) PLpro (− 8.2 kcal/mol) (B) NSP12 (− 10.2 kcal/mol) (C) NSP14 (− 11.8 kcal/mol) and (D) NSP16 (− 11.4 kcal/mol). The docking score obtained from Autodock Vina for each protein are given in squared bracket.

Ledipasvir-NSP12 complex exhibited three C-H….O bond, one π-σ bond, one π-anion, one π-donor H-bond, and four alkyl interactions (Table S9, SI and Figure 9B). However, one π-σ interaction with Ala688, one π-anion with Asp761, and two alkyl interactions with Lys500 and Ala685 were retained after 100 ns MD simulations (Table S9, SI). Additionally, two new C-H bonds with Arg569 and Tyr689 were observed in the ledipasvir-NSP12 complex after MD simulations. But ledipasvir-NSP12 complex missed π-donor H-bond after molecular dynamics.

In the ledipasvir-NSP14, initially, four H-bonds, one C-H….O bond, two alkyl, eight π-alkyl, one π-π T-shaped bond and one halogen interactions were present (Table S9 (SI) and Figure 9C). However, five π-alkyl interactions with Pro429, Val290, Pro335 and Trp292, one C-H bond with Try6930, one π-π T-shaped bond with Phe367 were retained after MD simulations. One new alkyl interaction with Trp385 and two new C-H bonds with Val389 and Phe426 were observed in the post-MD ledipasvir-NSP14 complex.

In the ledipasvir-NSP16 complex, initially, there are seven H-bonds, one C-H….O bond, two π-σ bonds, two π-alkyl, one π-donor H-bond and one halogen interactions were present (Table S9 (SI) and Figure 9D). However, only three H-bond with Asp6873, Asp6912, and Phe6947, one π-σ bond with Leu6898, and one halogen interaction with Tyr6930 were engaged after MD simulations. But three H-bonds, one π-σ bond, π-alkyl and π-anion bonds were not retained in the post-MD complex. Additionally, few new interactions were examined such as four H-bond with Lys6844, Asn6841, Lys6898 and Glu6945, one π-σ bond with Met6929, and one C-H….O bonds with Glu7001 during the MD simulation of the ledipasvir-NSP16 complex. To compare with the known drugs, a study has been undertaken with elbasvir and nafamostat, which reveal that interaction with the active site residues are virtually the same as that of the 4 other selected candidates in the aforementioned section (the corresponding 2D interaction diagram of elbasvir and nafamostat has been illustrated in Figures 10, 11).

Figure 10

2D representation depicting docking interactions obtained from Autodock Vina of elbasvir with SARS-CoV-2 protein targets (A) 3CLpro (− 8.1 kcal/mol) (B) PLpro (− 8.3 kcal/mol) (C) NSP12 (− 9.5 kcal/mol) (D) NSP13 (− 9.6 kcal/mol) (E) NSP14 (− 11.2 kcal/mol) and (F) NSP16 (− 9.3 kcal/mol).

Figure 11

2D representation depicting docking interactions obtained from Autodock Vina of nafamostat with SARS-CoV-2 protein targets (A) 3CLpro (− 8.3 kcal/mol) (B) PLpro (− 7.7 kcal/mol) (C) NSP12 (− 8.1 kcal/mol) (D) NSP13 (− 7.6 kcal/mol) (E) NSP14 (− 10.1 kcal/mol) and (F) NSP16 (− 9.2 kcal/mol).

To examine the nature of interactions evolving during molecular dynamics simulations, we ventured on collecting frames at the interval of 10 ns, for each of the 18 complexes. The 180 frames, thus retried were subjected to docking with two objectives in mind; a) the stability of the complex and b) the variation of binding energy during the MD simulation. Figure 12 illustrates the variations of the binding energies of all the complexes, while the detailed analysis of the exact count of the non-covalent interactions[76,77] at each of the frames was documented in the supporting information (Table S6-S9, SI). The electrostatic potential map of all the drug molecules along with the interacting amino acids are depicted individually in the six figures given in the supporting information (Figures S5-S10, SI).

Figure 12

The docking score distributions of different conformers generated from regular time intervals of MD trajectories of (A) venetoclax, (B) tirilazad, (C) acetyldigitoxin and (D) ledipasvir complexes respectively. About 10 poses were extracted from the 100 ns trajectories by taking a pose at interval of 10 ns. The drug molecules were removed and protein was prepared without changing their conformation. The respective drug molecules were docked in the active site of 18 * 10 conformer.

From the docking study carried out at an interval of 10 ns from 100 ns trajectories, it was observed that among the six protein targets interacting with venetoclax, NSP16, and NSP14 have shown substantial variations in the docking scores along the trajectory, which may be traced to the fluctuating number of hydrogen bonds and other non-covalent interactions (Figure 12). The current exercise was carried out to see whether the protein-ligand complexes are intact or not during the simulations. Clearly, it was observed that in most cases the protein-ligand complexes are intact, albeit mostly losing some of the binding interactions as a function of time. While most of the protein-drug complexes retain the bulk of the binding affinity, even after the 100 ns simulation, the venetoclax complex with NSP16 drifting of the docking score may be traced to the loss of as many as 3 hydrogen bonds during the time evolution of the trajectory. The variations in the docking scores can be explained by comparing with the loss (mostly) or gain of the number of hydrogen bonds during the various frames of the MD trajectory (see Table S6-S9, SI).

MM-PBSA binding free energy calculation

The binding affinity of the drug obtained from Autodock Vina is based on the single protein-drug configuration corresponding to the least energy binding pose of the drug within the binding site of a target protein. However, the drug can have multiple effective conformations within the binding site and the reliable binding free energy estimate should be obtained as the sum over these. The binding free energies computed using MM-PBSA[78] the multi-targeting potential of this target approach corresponds to average binding free energies obtained over various protein-ligand configurations from MD trajectories. Even though the calculation of binding free energies from MM-PBSA approach is computationally very demanding, this approach has been successful in ranking protein-ligand complexes.[78-80] However, it is worth noting that the MM-PBSA approaches were shown to yield binding free energies which are quantitatively much larger when compared to docking energies.[79,80] MM-PBSA approach has also been used to identify various lead compounds for COVID-19. However, the current study differs in the aspect that we aim at identifying drugs that are having multi-targeting potential.

Below we discuss the binding free energies computed for the four selected drugs with their high-affinity targets. A recent study on remdesivir interaction with RNA dependent RNA polymerase using MM-PBSA approach reported binding free energy of −78.5 kcal/mol.[81,82] Since remdesivir is one of the recommended drugs for the treatment of COVID-19, we use this value as a reference to comment on the relative binding free energies for the drugs studied in this work. The binding free energies computed using MM-PBSA approach for venetoclax further confirms the multi-targeting potential of this target. Table 2 shows the computed binding free energies for this drug with 6 different targets and binding free energies are in the range -56 kcal/mol to -103 kcal/mol. If we analyze the individual contributions to total binding free energies it can be clearly established that the protein-drug interactions are majorly driven by hydrophobic interactions. In general, the sum of electrostatic and polar solvation energy is positive and so electrostatic interactions are acting against the complex formation. The non-polar solvation free energies are negative suggesting that they contribute to the stabilization of protein-drug complex formation. Based on the reference value, it is observed that the drug, venetoclax has superior binding affinity for the 4 targets such as 3CLpro, NSP14, NSP12, NSP13 (as the binding free energies were lower than the reference value, −78.5 kcal/mol) while it has moderate affinities for the remaining two targets. Similarly, the tirilazad showed significantly larger binding affinities for the two targets PLpro and NSP16 when compared to the remaining targets, NSP12 and 3CLpro. Ledipasvir showed larger binding affinities for the targets, NSP12, NSP14 and NSP16 when compared to the target, PLpro. However, the drug acetyldigitoxin showed only moderate binding affinities for the four targets 3CLpro, NSP13, NSP14 and NSP16. Even though this drug showed superior binding affinities for all targets based on the scoring function as implemented in Autodock Vina, with the rescoring using MM-PBSA binding free energies, we only observe moderate binding. The binding free energies of these molecules are essentially in the same range, compared to that of the known drug molecules elbasvir and nafamostat. The validation study by taking known drug molecules as reference ones increase the confidence levels in the application of these models. Thus, the first level screening using Autodock Vina and rescoring using more reliable scoring functions based on MM-PBSA binding free energies can be a potential protocol to remove false positives and to improve the reliability of the results lead compounds. The study suggests that the multilevel scoring approach using two or more scoring functions (with varying accuracy) sequentially can be effectively used to identify lead compounds for therapeutics. Overall, the current study using Autodock Vina and MM-PBSA based scoring functions suggest that the drugs, venetoclax, tirilazad and ledipasvir have multi-targeting potential and can be used to combat the current COVID-19 challenge which is posing a major threat to world healthcare systems due to its ability to mutate at a faster rate. Table 2 displayed the average MM-PBSA energy (kcal/mol) terms for the binding of four multi-targeting candidates.

Table 2

MM-PBSA analyses of the four selected drug candidates against the multiple targets of SARS-CoV-2.

Drug name	Protein	vdW. energy	Elec. energy	Polar solvation energy	SASA energy	Binding energy
Venetoclax	3CLpro	− 187.08	− 24.72	129.59	− 21.10	− 103.32
	PLpro	− 133.73	− 28.05	120.60	− 15.06	− 56.24
	NSP12	− 131.44	− 25.94	81.83	− 13.03	− 88.58
	NSP13	− 278.11	− 52.01	275.11	− 33.15	− 88.17
	NSP14	− 245.48	− 46.32	174.52	− 24.35	− 141.63
	NSP16	− 188.05	− 39.82	181.72	− 21.29	− 67.447
Tirilazad	3CLpro	− 121.73	− 27.08	110.30	− 13.06	− 52.59
	PLpro	− 151.48	− 6.45	72.31	− 18.39	− 104.01
	NSP12	− 102.36	− 9.23	56.26	− 10.32	− 30.68
	NSP16	− 135.24	− 18.43	77.64	− 15.40	− 91.44
Acetyldigitoxin	3CLpro	− 94.62	− 32.69	99.41	− 12.83	− 40.73
	NSP13	− 100.81	− 39.69	117.42	− 14.40	− 37.28
	NSP14	− 59.15	− 14.75	59.30	− 8.70	− 23.30
	NSP16	− 38.70	− 12.68	64.74	− 10.45	− 55.57
Ledipasvir	PLpro	− 114.89	− 28.08	87.50	− 14.59	− 70.03
	NSP12	− 241.98	− 79.58	207.79	− 29.28	− 116.06
	NSP14	− 273.26	− 48.60	223.06	− 33.70	− 132.50
	NSP16	− 259.25	− 100.71	266.84	− 29.25	− 122.38
Elbasvir	3CLpro	− 167.08	− 21.72	119.59	− 19.20	− 99.32
	PLpro	− 158.14	− 27.30	118.92	− 17.14	− 83.65
	NSP12	− 194.16	− 40.30	178.91	− 25.30	− 80.85
	NSP13	− 194.25	− 54.67	185.75	− 22.94	− 86.11
	NSP14	− 294.48	− 96.95	306.26	− 35.04	− 120.21
	NSP16	− 198.05	− 49.83	191.72	− 29.30	− 77.54
Nafamostat	3CLpro	− 38.44	− 5.63	20.90	− 5.90	− 28.45
	PLpro	− 0.002	21.87	13.50	− 0.19	35.18
	NSP12	− 123.44	− 23.64	71.32	− 11.43	− 78.48
	NSP13	− 123.10	314.19	394.97	− 19.55	566.50
	NSP14	− 92.99	− 82.88	112.56	− 12.90	− 76.22
	NSP16	− 168.05	− 19.82	161.72	− 11.29	− 47.447

Each of the 18 complexes was subjected to 100ns MD simulation. The energy terms (in kcal/mol) were calculated from the data obtained from the last 20 ns trajectory.

vdW.: Van der Waals; Elec.: Electrostatic; SASA: solvent accessible surface area.

Venetoclax interaction with the SARS-CoV-2 targets

After considering a number of factors discussed in the foregoing sections, one may come to the conclusion that venetoclax has been a highly potential drug candidate from polypharmocology point of view, and essentially viral infections it can be a very promising option. A differential analysis has been carried out in order to validate and examine the robustness of the study carried out. As a first step, we have calculated the distance between the interacting catalytic residues with the venetoclax. The average distance was in between 0.2 nm–3 nm. In case of Plpro, the distance has been fluctuated between 2 nm–4 nm. This analysis clearly shows that the venetcolax has been bound well within the binding pocket (Figure S11, SI). Moreover, we have performed the PCA – FEL analysis for the venetoclax MD complexes. The first 10 eigenvectors projections of proteins were extracted and analyzed for their cosine content and FEL was defined based on PC1 and PC2 having a cosine content less than or equal to 0.2. The FEL of 6 venetoclax complexes was mapped from two PCs based on cosine content analysis in order to get minimum energy and energetically favored conformations. The energetically favored conformed were retrieved from the MD trajectories and interactions were analyzed and depicted in Figure S12, SI.

Conclusions

In the present study, various molecular modeling strategies were applied to identify promising multi-targeting repurposable potential drug candidates against SARS-CoV-2 infection. An effective polypharmacology approach has been employed in this study, by considering seven of the essential proteins involved in SARS-CoV-2 infection, 3CLpro, PLpro, RdRp (NSP12), NSP13, NSP14, NSP15, and NSP16, to develop efficient therapeutics capable of confronting the resistant strains. Some of the polypharmacological hits reported in this study have never been considered for COVID19 treatment. After identifying the known drugs which displayed high promiscuity in binding with several of the druggable targets, we tried to filter the top 20 hits from each target. Four drugs namely, venetoclax, tirilazad, acetyldigitoxin and ledipasvir were selected based on high docking score, visual examination of interactions, good binding pose and drugs showing interaction with four or more targets. The extensive MD simulation and MM-PBSA analysis showed the stability of protein-drug complexes and the sustainability of the key interactions with the respective targets. The van der Waals and electrostatic energies were found to have the prime driving forces to the overall negative binding free energy and stabilities. Our proposed repurposing candidates are exhibited multi-targeting potential which may assist the experimental study with higher chances to identify new antiviral drugs against SARS-CoV-2. The current study reveals that targets of key importance in COVID-19, display substantial binding with several of the known drugs, which indicates that polypharmacology approaches may be optimal choices. While vaccine development is indispensable, there is tremendous pressure to develop efficient anti-viral molecules, especially as one witnesses the emergence of a number of immune escape variants.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 2526 kb)