In-silico drug repurposing for targeting SARS-CoV-2 main protease (Mpro).

COVID-19, caused by novel coronavirus or SARS-CoV-2, is a viral disease which has infected millions worldwide. Considering the urgent need of the drug for fighting against this infectious disease, we have performed in-silico drug repurposing followed by molecular dynamics (MD) simulation and MM-GBSA calculation. The main protease (Mpro) is one of the best-characterized drug targets among coronaviruses, therefore, this was screened for already known FDA approved drugs and some natural compounds. Comparison of docking and MD simulation results of complexes of drugs with that of inhibitor N3 (experimentally obtained) suggests EGCG, withaferin, dolutegravir, artesunate as potential inhibitors of the main protease (Mpro). Further, in silico docking and MD simulation suggest that EGCG analogues ZINC21992196 and ZINC 169337541 may act as a better inhibitor.Communicated by Ramaswamy H. Sarma.

Introduction

A new coronavirus, named SARS-CoV-2, has started from Wuhan, China and spread over the world (Shereen et al., 2020; Wu et al., 2020; Zhou et al., 2020). It is a severe acute respiratory syndrome-related coronavirus (Gorbalenya et al., 2020). The coronavirus disease, COVID-19 has affected 213 countries and territories around the world, among them USA, Brazil, Russia, India, UK, Spain, Italy and Peru are the worst-hit countries having over 4.89 million infections combined. As of June 14, 2020, around 7.89 million people have been infected, and 0.432 million have died (https://www.worldometers.info/coronavirus/). The number of corona cases is increasing exponentially and not a single vaccine or drug is available for the treatment. Initially, antimalarial drugs chloroquine and hydroxychloroquine (Principi & Esposito, 2020; Singh et al., 2020) was proposed as treatment, however, recent studies have not shown any advantages and is associated with some side effects which can be life-threatening. The development of a vaccine will still one year away and there is a strong need for repurposing an already known drug which can block the activity of virus.

Coronaviruses are a large family of RNA viruses that are encapsulated inside a membrane envelope. Envelop have proteins appearing like spikes sticking out from their surface (Walls et al., 2020). Main protease (Mpro) is an enzyme involved in processing of polyprotein which is translated from viral RNA. Inhibition of this protease will block the replication of the virus and thus is an important drug target. Michael acceptors can act as cysteine protease inhibitors, one such inhibitor is N3. It can reversibly bind to Mpro which then blocks its activity by forming a covalent bond with catalytic cysteine after the nucleophilic attack (Jin et al., 2020; Santos & Moreira, 2007). The crystal structure of the Mpro of SARS-CoV-2 with N3 inhibitor has been recently reported, and this can be taken as a control for discovering new inhibitors. Alpha-ketoamide based inhibitors have also been reported; one of the inhibitors has shown low micromolar EC-50 against SARS-CoV-2 (Zhang et al., 2020).

Mpro is a homodimer, each monomer consisting of three domains. Domain 1 consists of residues 10–99, domain 2 consists of residues 100–182 and residues 198–303 form domain 3 of Mpro. Domain 3 is responsible for regulating the dimerization of Mpro. Catalytic dyad of protein is formed by residues His 41 and Cys145. Dimerization of enzyme is necessary for catalytic activity since it helps in making of S1 pocket of the substrate binding site (Anand et al., 2003; Jin et al., 2020). Thus, any drug-like candidate which strongly binds to S1 site or inhibits the dimerization process can potentially inhibit the replication of virus. Docking and molecular dynamics (MD) simulation have emerged as a tool to predict the putative drugs if the target protein structure is available.

Recently, several research groups have reported some inhibitors of SARS-CoV-2 main protease (Adem et al., 2020; Khan et al., 2020; Ton et al., 2020; Verma et al., 2020; Zhang et al., 2020). Antiviral drugs and natural compounds are some of the widely explored class of compounds due to their good binding interaction with the main protease. In our study, we have carried out in-silico screening of a range of FDA-approved antiviral, antimalarial drugs, some Ayurveda medicines and some natural polyphenols targeting Mpro. We have arranged the obtained hits according to their interaction energies and top compounds have been subjected to MD simulations in much more realistic environment. MM-GBSA methodology was used to obtain more accurate interaction energies and they are compared with the interaction energy of N3 with Mpro.

Materials and methods

Preparation of protein structure and ligand database for docking

The crystal structure of the free enzyme of the SARS-CoV-2 main protease (PDB: 6Y2E) was recently published by Zhang et al. (2020). Its substrate-binding pocket was chosen for the screening of compounds as obtained from Mpro-N3 complex (Figure S1). The ligand database was prepared by retrieving the structures of FDA approved antiviral and antimalarial drugs from Drugbank (Wishart et al., 2006), structures of Ayurveda compounds from PubChem (Kim et al., 2019) and structures of polyphenols were drawn using Marvin Sketch (MarvinSketch 17.23.0, ChemAxon). Structures of EGCG scaffold were obtained from the Zinc database (Sterling & Irwin, 2015).

Figure 1.

(a) Cα root mean square deviation plot for last 500 ns of MPro-Drug complex, (b) Average RMSD values for last 300 ns of MPro-Drug complex.

Docking analysis

AutoDock Vina (Trott & Olson, 2010) was used for the screening of ligand database. AUTODOCK 4 (Morris et al., 2009) was used to prepare the ligand structure for docking. Docking generates several poses of ligand inside the binding pocket, the poses showing minimum binding energy were chosen for carrying out further studies. OpenBabel (O’Boyle et al., 2011; Open Babel-4.2.1) was used for obtaining PDB structures of ligands from PDBQT and for adding the protons. Various interactions like hydrophobic and hydrogen bonding interactions between ligand and protein were visualized using LigPlot (Laskowski & Swindells, 2011).

MD simulations

MD simulations were carried out on Amber18 package (Case et al., 2018). Ligand parameters were obtained from general Amber force field (GAFF2) (Träg & Zahn, 2019) and AM1 BCC method was used to derive the charges. Bonded and non-bonded parameters for protein were obtained from the Amber 14SB force field. Docked protein–ligand complexes were solvated in a 10 nm truncated octahedron with TIP3P water molecules (Jorgensen et al., 1983) and counter ions were added to neutralize the system. The study utilized periodic boundary conditions and PME summation (Petersen, 1995) for electrostatic calculations.

The shake methodology was applied to restrict covalently bonded hydrogen atoms. Constant pressure condition was maintained using Berendsen thermostat (Berendsen et al., 1984). Time step of 2 femtosecond with 9 Å cut off for non-bonded interactions were applied. Protein–ligand complexes were energy minimized in two steps, first using 250 steps of steepest descent followed by 750 steps of conjugate gradient method. Equilibration was also carried out in two steps, heating the systems at 300 K followed by simulation of complexes by decreasing the force slowly up to 0.1 N. Next, fully unrestricted equilibration was performed for 5 ns. The convergence of various system properties was monitored. The production run was performed for 1500 ns under NPT ensemble. Simulation of protein without any ligand was also carried out. For Mpro-N3 complex, available crystal structure (PDB: 6LU7 (Jin et al., 2020)) was used and simulation was carried out as mentioned above. 500 ns simulations of protein were also performed with compounds having scaffold of EGCG.

Analysis of MD simulation trajectory

For analysis, Amber trajectory was converted to Gromacs (Abraham et al., 2015) trajectory using python script. For checking the stability of protein due to ligand binding, root mean square deviation (RMSD) was monitored with respect to the reference structure (first frame) for Cα atoms of protein. In order to see the residue level fluctuations, root mean square fluctuation (RMSF) was observed for last 300 ns for Cα atoms of protein. For obtaining distance distribution, distances between center of mass of active site residues and center of mass of ligands were calculated for last 300 ns and distribution was obtained from the histogram analysis tool of MS-Excel. Clustering analysis was done by taking the last 300 ns trajectory, using the single linkage method. In this, a cluster was formed by subsequent addition of structures whose distance is less than 0.1 nm to any element of the cluster, if the distance is more than 0.1 nm a new cluster was formed.

MM-GBSA approach

MM-GBSA method was employed for calculating the free energy difference between bound and unbound form of receptor and for this MMPBSA.py tool of Amber18 was utilized. MMPBSA.py is a post-processing tool which extracts snapshots of complex, protein and ligand from the same trajectory of solvated complex. Free energy was calculated according to the following thermodynamic cycle (Scheme 1):

Scheme 1.

Thermodynamic cycle utilized for MM-GBSA calculation.

Solvation free energies have contributions from electrostatic and non-electrostatic interactions; electrostatic contribution was obtained by solving the Generalized Born (GB) equation. Here, a modified GB model (Onufriev & Case, 2019) was used for electrostatic part calculations. Non-electrostatic part was considered to be proportional to the total solvent accessible area of the molecule and was obtained from the linear combination of pairwise overlaps (LCPO) algorithm (Weiser et al., 1999). Free energy change in vacuum was obtained by taking average interaction energy of protein and ligand from molecular mechanics and considering entropy change upon binding. For the comparison of the binding energies of various ligands to the same receptor, entropy contributions were not taken into account.

Result and discussion

In-silico screening of antiviral, antimalarial drugs, polyphenols and ayurvedic compounds through molecular docking with Mpro

In order to find the potential binder to the substrate site of Mpro (PDB: 6Y2E (Zhang et al., 2020)), a total of 125 FDA approved drugs and compounds were screened utilizing molecular docking studies. Out of the total, 31 compounds have binding energy below −6 kcal/mol, 52 compounds have binding energy between −6 and −7 kcal/mol, 35 compounds have binding energy between −7 and −8 kcal/mol and only 7 compounds have binding energy above and equal to −8 kcal/mol (Table S1). Compounds having binding energies above or equal to −8 kcal/mol were selected for carrying out MD simulation to understand their molecular level interactions with the receptor, their binding energy and interacting residues are listed in Table 1 and their structures are presented in Figure S2. The selected seven compounds were Epigallocatechin gallate (EGCG), Delavirdine, Dolutegravir, Indinavir, Artesunate, Tinosporin B and Withaferin A. Their interactions with the receptor are shown in Figure S4. Epigallocatechin gallate (EGCG) is a polyphenolic compound found in green tea and has been shown to have antiviral activity. The calculated binding energy between EGCG and Mpro from docking is −8.3 kcal/mol and EGCG.

Table 1.

Binding parameters and interacting residues for Mpro-ligand complexes obtained from docking simulations.

S. No.	Ligand	Binding energy (kcal/mol)	Interacting residues
1	Tinosporin B	−8.0	25, 41, 44–46, 49, 142, 145, 164–166, 189
2	EGCG	−8.3	25, 41, 46, 49, 140, 141, 144, 142, 145, 165, 166, 189
3	Withaferin	−8.9	24, 25, 41, 44–46, 49, 140–142, 164, 166
4	Delavirdine	−8.1	25–27, 41, 46, 49, 140–143, 145, 163, 166
5	Dolutegravir	−8.6	24, 25, 41, 45, 46, 142, 143, 145, 166
6	Artesunate	−8.0	24–27, 41, 46, 49, 140–145, 163, 165, 166
7	Indinavir	−8.0	25, 26, 41, 49, 140–142, 145, 163, 164, 166, 189

Figure 2.

Cα RMSF plot for (a) N3, (b) tinosporin B, (c) EGCG, (d) withaferin A, (e) delavirdine, (f) dolutegravir, (g) indinavir and (h) artesunate complexes.

Figure 4.

EGCG analogs obtained from screening for Mpro.

Delavirdine, Dolutegravir and Indinavir are FDA approved antiviral drugs used for the treatment of Human Immunodeficiency virus (HIV). Binding energies obtained for delavirdine, dolutegravir and indinavir are −8.1 kcal/mol, −8.6 kcal/mol and −8.0 kcal/mol, respectively.

Artesunate is one of the most rapid-acting antimalarial drugs and it is a semi-synthetic derivative of artemisinin. The calculated binding energy for protein–artesunate complex is −8 kcal/mol.

Two Ayurveda compounds also showed significant binding to the target, one is tinosporin B which is a diterpenoid furanolactones extracted from plant tinospora cardyfolia. Tinosporin has anti-inflammatory, anti-microbial, antihypertensive and antiviral properties. Another Ayurveda compound which qualified our binding energy criterion is Withaferin A, it is extracted from the leaves of the Indian plant Withania somnifera (ashwagandha) and exhibits chemo preventive, anti-cancer and immunomodulatory actions. Calculated binding energies for tinosporin A and withaferin A were found to be equal to −8.0 kcal/mol and −8.9 kcal/mol, respectively.

MD simulations were performed to ascertain the stability of docked complexes. RMSD values of Cα were plotted with respect to the first frame for all the complexes (Figure 1(a)). For a better comparison, the average values of RMSD for last 300 ns were plotted along with standard deviation (Figure 1(b)). For protein with N3 inhibitor, the RMSD value was observed to be equal to 0.34 nm (±0.04). The RMSD value for the Artesunate complex increased from 0.25 nm (±0.01) for protein without ligand to 0.59 nm (±0.02), which indicated increase in structural fluctuations due to ligand binding. Least RMSD changes were observed for the EGCG complex with a value of 0.22 nm (±0.03) showing the stable docked complex formation. For others, RMSD values remain in the range of 0.24–0.38 nm. Tinosporin B (0.27 nm), EGCG, Withaferin A (0.28 nm), Delavirdine (0.28 nm) and Dolutegravir (0.24 nm) complexes showed smaller Cα RMSD values as compared to N3 complex, whereas, Indinavir (0.38 nm) and Artesunate exhibited larger Cα RMSD values. After 1100 ns, Cα RMSD values for all the complexes were not changing much, signifying a stable complexation between protein and ligand.

Next, residue level fluctuations were seen through root mean square fluctuation (RMSF) analysis (Figure 2). In N3-complex, for all the residues, except for terminal residues, RMSF values were not changing much and values observed below 0.3 nm (Figure 2(a)). Increased fluctuations were observed for the residues Gln 189, Thr 190, Ala 191, Gln 192, Ala 193 and 194 as compared to N3 bound protein. This trend was observed for all the complexes. For EGCG, withaferin A, dolutegravir and artesunate complexes (Figure 2(c–e,h)), RMSF values were found to be similar to N3 bound protein. For tinosporin B, delavirdine and indinavir, larger fluctuations were also observed for the residues 44–60, which includes some ligand interacting residues Cys 44, Thr 45, Ser 45 and Met 49 (Figure 2(b,f,g)). Our RMSF analysis suggests that upon binding of tinosporin B, delavirdine or indinavir, fluctuation is more in a small region of protein and EGCG, withaferin A, dolutegravir and artesunate complexes remain stable.

Furthermore, to ensure the stable binding of ligand inside the binding pocket during simulation, distances between center of mass of ligand and center of mass of active site residues His 41,Cys 145,Glu 166 and Gln 189 were calculated for all the complexes for last 300 ns and frequency distributions were plotted (Figure 3). Distribution was centered within 1 nm for N3, EGCG, Withaferin A, Dolutegravir and Artesunate for all the four residues, whereas, broad distance distributions were observed for Tinosporin B, Delavirdine and Indinavir complexes above 2 nm. A broad distribution at higher distance values can be due to weak protein–ligand complexation at the active site. Therefore, only EGCG, withaferin A, dolutegravir and artesunate show tight binding to the active site based on the distance distribution analysis.

Figure 3.

Distance distribution between ligands and (a) residues 41, (b) residue 145, (c) residue 166 and (d) residue 189.

Subsequently, clustering analysis was done on the last 300 ns trajectory and a representative structure was obtained which then visualized through LigPlot to see the change in molecular level interaction during simulations (Figure S3).

The obtained representative structures from clustering analysis were compared with the docked structures before the start of simulation, for the active site interactions (Figure S3). For Delavirdine and indinavir, no active site interactions were seen and binding site for these ligands were changed completely. Tinosporin B with one intact contact also showed weak binding interactions. For protein–EGCG complex, number of active site interactions increased to six s, indicating a good binding of ligand to the receptor. Withaferin A and dolutegravir showed four and two active site interactions, respectively, and also seen to be well within their desired binding pockets. Artesunate exhibited five active site interactions. Therefore, tinosporin B, indinavir and delavirdine have been ruled out to be tested as a drug candidate since these were not stable inside the binding pocket. EGCG, withaferin A, dolutegravir, artesunate and their analogs can be considered as good antidotes for COVID-19 and further studies can be carried out on them.

Screening of EGCG analogs and MD studies

The binding of EGCG to the active site of Mpro is found to be efficient; therefore, all the compounds having scaffold of EGCG were also screened through molecular docking and MD simulations. A total of 132 compounds were docked to Mpro and three compounds (Figure 4) showing highest binding energies were subjected to 500 ns MD simulation. Their binding energies are given in Table 2. RMSD and clustering analyses were performed (Figures S4 and S5). Stable Cα atoms RMSD and conserved binding pocket during simulations indicate their ability to serve as inhibitors of main protease.

Table 2.

Binding energies obtained from MM-GBSA method.

S. No.	Complex	Binding Free energy (kCal/mol)
1	N3	−28.79 ± 3.73
2	Tinosporin B	−13.04 ± 7.90
3	EGCG	−25.53 ± 4.18
4	Withaferin A	−29.89 ± 4.67
5	Delavirdine	−11.13 ± 8.53
6	Dolutegravir	−18.70 ± 6.43
7	Indinavir	−19.00 ± 12.38
8	Artesunate	−24.72 ± 5.01
9	ZINC3870415	−18.52 ± 5.53
10	ZINC21992196	−35.32 ± 6.75
11	ZINC169337541	−28.72 ± 3.67

MM-GBSA approach for calculating binding free energy of complex

MM-GBSA approach utilizes molecular mechanics, the generalized Born (GB) model and solvent accessibility method for free energy calculations. Binding free energies for all the complexes including Mpro-EGCG analogs were calculated using MM-GBSA approach (Table 2). For N3, tinosporin B, EGCG, withaferin A, delavirdine, dolutegravir, indinavir and artesunate complexes, last 300 ns trajectory from 1500 ns MD trajectory was taken and for complexes of EGCG analogs, last 100 ns trajectory from 500 ns was selected for this calculation. For N3-Mpro complex, binding energy was found to be equal to −28.79 (±3.73) kcal/mol. Among all the complexes, withaferin A, EGCG, dolutegravir, artesunate, ZINC21992196 and ZINC169337541 showed significant binding energy when compared with Mpro-N3 complex. Tinosporin B, delavirdine, indinavir and ZINC3870415 showed low binding energies. The similarity ensemble approach (SEA) predicted the probable targets for ZINC21992196 and ZINC 169337541, which indicates that these compounds may act as protease inhibitors (Keiser et al., 2007). Our MM-GBSA results support our previous findings and predict that EGCG, withaferin A, dolutegravir and artesunate may act as potential inhibitor of SARS-CoV-2 Mpro. Favorable interaction of Mpro with EGCG analogs also propose that compounds or drugs having same scaffold can also be tested for their inhibitory actions.

Conclusion

Our docking and MD studies suggest EGCG, withaferin A, doutegravir, artesunate and their analogs may act as an inhibitor for the main protease of SARS-CoV-2. EGCG which is a major component of green tea is capable of showing antiviral activity against many DNA and RNA viruses, like adenovirus, influenza virus, human immunodeficiency virus (HIV) and coronavirus, by inhibiting various stages of viral infection (Benelli et al., 2002; Ho et al., 2009; Nance et al., 2009; Song et al., 2005; Weber et al., 2003). Similarly, withaferin A which is a component of ayurvedic plant Withania somnifera (ashwagandha) also exhibits antiviral properties against influenza virus (Cai et al., 2015). Artesunate, a semisynthetic derivative of artemisinin, is a well-known malaria drug and has the potential to fight against many viruses (Efferth et al., 2008; Lisewski et al., 2014). Dolutegravir is an antiretroviral drug used along with other medications to treat HIV/AIDS (Min et al., 2010). Supported by previous literature and our results we propose EGCG, withaferin A, dolutegravir and artesunate as potential drugs for COVID-19 and the compounds having similar scaffolds can also be tested for their potency. Furthermore, experimental work needs to be done to verify their use as a drug.

Click here for additional data file.