Identification of potential molecules against COVID-19 main protease through structure-guided virtual screening approach.

The pandemic caused by novel coronavirus disease 2019 (COVID-19) infecting millions of populations worldwide and counting, has demanded quick and potential therapeutic strategies. Current approved drugs or molecules under clinical trials can be a good pool for repurposing through in-silico techniques to quickly identify promising drug candidates. The structural information of recently released crystal structures of main protease (Mpro) in APO and complex with inhibitors, N3, and 13b molecules was utilized to explore the binding site architecture through Molecular dynamics (MD) simulations. The stable state of Mpro was used to conduct extensive virtual screening of the aforementioned drug pool. Considering the recent success of HIV protease molecules, we also used anti-protease molecules for drug repurposing purposes. The identified top hits were further evaluated through MD simulations followed by the binding free energy calculations using MM-GBSA. Interestingly, in our screening, several promising drugs stand out as potential inhibitors of Mpro. However, based on control (N3 and 13b), we have identified six potential molecules, Leupeptin Hemisulphate, Pepstatin A, Nelfinavir, Birinapant, Lypression and Octreotide which have shown the reasonably significant MM-GBSA score. Further insight shows that the molecules form stable interactions with hot-spot residues, that are mainly conserved and can be targeted for structure- and pharmacophore-based designing. The pharmacokinetic annotations and therapeutic importance have suggested that these molecules possess drug-like properties and pave their way for in-vitro studies.Communicated by Ramaswamy H. Sarma.

Introduction

The coronavirus (COVID-19) is a newly emerged human-infectious coronavirus (CoV), pandemic and a global health emergency. Unfortunately, at present there is no well-defined treatment or therapeutics against COVID-19 is available but the preventive measures are being recommended worldwide. However, the clinical trials for already marketed drugs such as lopinavir, ritonavir, hydroxychloroquine, azithromycin, (Tirumalaraju, 2020c) chloroquine (ClinicalTrials.gov, n.d.), Remdesivir (Tirumalaraju, 2020b) etc. along with antibiotics are being evaluated to treat the secondary infections (www.clinicaltrials.gov). All of the drug options come from experience treating SARS, MERS or some other new influenza virus previously (Lu, 2020). These drugs would be helpful but the efficacy needs to be further confirmed. Few COVID-19 vaccines are also under clinical trials such as Moderna’s mRNA-1273, first US clinical vaccine funded by NIH’s NIAID (National Institute of Allergy and Infectious Diseases) (Tirumalaraju, 2020a). Thus, there is an unmet requirement for the specific anti-COVID-19 therapeutics to limit the severity of the deadly disease. Various clinicians and researchers are engaged in investigating and developing antivirals using different strategies combining experimental and in-silico approaches (Elfiky, 2020; Enayatkhani et al., 2020; Enmozhi et al., 2020; Islam et al., 2020; Jin et al., 2020; Khan, Jha, et al., 2020; Khan, Zia, et al., 2020; Qamar et al., 2020; Sinha et al., 2020) with the goal of identifying novel, selective and potent therapeutic agents.

An attractive drug target among coronaviruses is the main protease (Mpro, 3CLpro), due to its essential role in processing the polyproteins that are translated from the viral RNA (Boopathi et al., 2020). The present study focused on the main proteases in CoVs as potential target proteins for COVID-19 treatment. Mpro/3CLpro is active in its dimer state but till now there is no crystal structure available for the dimer form. Its monomer inhabits the 306 amino acids including 3 domains, folded into helices and β-strands. The electron density map for the monomer protein is clearly visible (Figure 1(A)). The domain I (residues 1–101) and II (residues 102–184) includes an antiparallel β-barrel structure; and domain III (residues 201–303) includes five α-helices arranged into a largely antiparallel globular cluster, and is connected to domain II by means of a long loop region (residues 185–200) (Jin et al., 2020). The catalytic dyad (H41 and C145) is responsible for the catalytic activity of SARS-COV-2 and is placed at the junction of domain I and domain II (Figure 1(A)).

Figure 1.

Assessment of APO and COM structures of SARS-COV-2. (A) Overview of the APO structure (PDB-ID: 6M03), (B) Superimposition of APO (Yellow) and COM1 (Cyan) structure with compound N3 represented in VdW, (C) Binding site overlay highlights the conformational differences in the residues. The residues are shown in licorice representation along with inhibitors shown in purple. The HBs are shown via red color dotted lines.

Recently, crystal structures for monomeric Mpro in both APO (PDB-ID: 6M03) and HOLO (PDB-ID: 6LU7, bound with N3 inhibitor) forms were crystallised. The Mpro of 2019-nCov shares 96% similarity with the Mpro of the SARS-CoV (Qamar et al., 2020). It is reported that 12 residues vary in both SARS-CoV-1 and SARS-CoV-2 but the residue S46 in SARS-CoV-2 (COVID-19) (corresponding residue A46 in SARS-CoV-1) is part of the binding pocket of the N3 molecule or active site (Qamar et al., 2020). Another co-crystal (PDB-ID 6Y2F/6Y2G), α-ketoamide inhibitor (13b) is also reported recently, providing the structural and residue-based architecture of catalytic sites (Zhang et al., 2020). These co-crystals are paving the route for the application of virtual screening (VS) to get more efficacious molecules (Al-Khafaji et al., 2020; Kumar et al., 2020; Lobo-Galo et al., 2020).

Knowledge gained from the previous outbreaks and existing antivirals gain attraction as the fastest route to fight the current coronavirus epidemic, henceforth this emergency put drug repurposing on fast track. Drug repurposing approach is being widely applied to quickly identify therapeutic solutions due to availability of their pharmacokinetic, toxicological and manufacturing data. It includes drugs that are either FDA approved, investigational, withdrawn or shelved compounds. Although there are studies of the repurposing and marketed drugs which proposed several candidates for SARS-CoV-2 treatment (Aanouz et al., 2020; Elmezayen et al., 2020; Pant et al., 2020). With this aim, we have first used molecular dynamics simulations to standardize our computational model specially focused on stable architecture of the binding site, which was used for VS to eventually facilitate the rapid identification of potent molecules. The findings from this study may provide an opportunity to explore these compounds for anti-COVID-19 therapeutics.

Materials and methods

Protein structure preparation

The crystal structures APO (6M03) and COM1 (6LU7) were optimized and then minimized using the Protein Preparation Wizard module of Maestro (Anang et al., 2018; Maestro, 2017) in which OPLS3 (Optimized Potentials for Liquid Simulations) force field was used (Jorgensen et al., 1996).

SiteMap analysis

The SiteMap (SiteMap, 2017) program of Schrodinger Suite was also used for calculating binding sites on crystal 6LU7. The method was implemented as an unbiased approach to undermine the presence of any secondary or allosteric binding site. SiteMap which is a ligand independent method, will also help in calculating the druggability of the identified site (Halgren, 2009; Mattapally et al., 2018; Srivastava et al., 2018; Thakur et al., 2020). The OPLS-2005 force field (Jorgensen et al., 1996) was employed, and a standard grid was used with 15 site points per reported site and cropped at 4.0 Å from the nearest site point.

Ligand selection and preparation

The ligand structures were taken from SELLEKCHEM database (http://www.selleckchem.com/), the DrugBank database (https://www.drugbank.ca/) and the Repurposing hub (https://clue.io/repurposing) (see Supplementary Figure 1). The ligands were prepared using Schrödinger’s (version 2017-1), LIGPREP (LigPrep, 2017), The optimization was done using the OPLS3 force field (Jorgensen et al., 1996). The known pharmacological activity of the hits were curated from the above-mentioned databases along with the PubChem database (https://pubchem.ncbi.nlm.nih.gov/).

Molecular dynamics simulation of Mpro

The APO of Mpro (PDB-ID: 6M03) along with COM1 and COM2 (docked 13b pose in Mpro monomer) were subjected to molecular dynamics simulation for 200 ns using Desmond v3.6 module from Schrodinger suite (Bowers et al., 2006). In addition, we also performed short simulation of hits molecules to observe their stability and its impact on protein. The systems were built via Systems builder using OPLS3 force field and solvated with TIP3P water solvent model. All the complexes were placed in the orthorhombic periodic boundary conditions with a size of repeating buffered units at 10 Å. Counter ions were also added to neutralize the systems. An energy minimisation step was done for each system using a steepest descent integrator for 2000 steps. The NPT ensemble was employed for the simulations with the Nose-Hover chain thermostat and the martyna-tobias-klein barostat. RESPA integrator was used with a time step of 0.002ps. For short range coulombic interactions, a 9.0 Å cut-off was considered. Bonds to hydrogen were constrained using the MSHAKE algorithm of Desmond. The coordinates were saved at intervals of 20 ps that are referred to as ‘frames’ in this study.

Virtual screening of virtual libraries on Mpro

The site of bound peptide ‘N3’ was chosen as primary site for ligand docking which is also confirmed by SiteMap as the most druggable site. The grid was generated using the centroid of N3 by using the Receptor Grid Generation panel in Glide. Docking studies were carried out using the VS Workflow (VSW) of Glide Schrodinger Suite (Friesner et al., 2006; Mittal et al., 2020).The ligands chosen from the database were passed through three stages of the screening workflow starting from high-throughput screening (50% filtered), followed by standard precision (30%) and finally extra precision (10%) stages. The final poses were processed using the Prime MM-GBSA panel at the end (Schrodinger suite, LLC, New York, NY, 2016-3) (Supplementary Figure 1). The control molecules N3 and 13b were also docked on the same grid. The results of the docking were then quantified on the basis of docking scores and MM-GBSA (ΔGbind). The pharmacokinetic properties of the molecules were evaluated by SwissADME server (Daina et al., 2017; Joshi et al., 2020). Also, the top filtered compounds were further cross-checked their stability through MD simulations as mentioned earlier.

Free energy analysis through MM-GBSA (molecular mechanics/generalized born surface area)

The average binding energy was calculated for equilibrated MD trajectory. The binding energy was calculated by where the difference in the minimized energies between ligand and protein complexes is denoted by ΔEMM. ΔGsolv is the difference in the GBSA solvation energy of the complexes and sum of solvation energies for the protein and ligand, whereas the differences in surface area energy of the complex and sum of that in protein and ligand (Kellici et al., 2019; Mittal et al., 2020).

Figures

All the images were generated using VMD and Schrodinger Suite (Asthana et al., 2014, 2015; Humphrey et al., 1996) and graphs were plotted using XMGRACE (Mittal et al., 2019; Turner, 2005).

Results and discussions

Comparative structural analysis

The availability of co-crystal 6LU7 in which a peptide (N3: a covalently-bonded inhibitor) is bound and knowledge of interacting residues with molecule 13b (α-ketoamide inhibitor) made the understanding of protein interesting. Since, the residues reported for compound 13b involve the same residues as for N3 molecules (Zhang et al., 2020). Therefore, monomer crystals were used for further analysis. The architecture of binding site with any possible conformational changes was concurred after superimposition of the crystal structures COM1 and APO. It shows that overall structure of protein is well aligned (RMSD APO-vs.-COM1 is 0.4 Å) except C-terminal region (Figure 1(B)). We also observed some noticeable differences in the binding site architecture of COM1 and APO (Figure 1(C)). It was found that all the residues of the binding site are well aligned with that of APO except the residues T25, M49, M165, R188, Q189 and T190 that shows their side chain conformational changes, while P168 shows the backbone movement also (Figure 1(C)). The inhibitor forms the hydrogen bond (HB) interactions with residues G143, H163, E166, Q189, and T190.

Exploring the druggability of binding sites, docking and benchmark setup for compound screening

To bring in more robustness in confirming the final binding site before performing the VS, we performed ligand independent binding site search on MPro. The top score of the SiteMap program also confirms the co-crystal site as the primary binding site (marked as Site 1) with the highest Dscore of 1.09 (best druggability score) (Figure 2 and Table 1). The volume of the pocket is 287.09 Å3. This binding site of Mpro is encompassed by domain I and II and the loop region connecting domains II and III. Sitemap result shows that this pocket is relatively smaller in size with a size score of 120 (reference value: 130), and more exposed to solvent with an exposure value of 0.614 (reference value: 0.52) (Halgren, 2009). Furthermore, we performed the focused docking on the Site 1 with the N3 and 13b molecules to set up the benchmark for VS execution based on two main criteria i.e. docking score and MM-GBSA values. The docking score and MM-GBSA values for N3 (COM1) is −10.6 kcal/mol and −64.32 kcal/mol, respectively, while for 13b (COM2) it is −6.66 kcal/mol and −63.26 kcal/mol, respectively (Table 2).

Figure 2.

Binding site identification: The possible binding sites and poses found by SiteMap. The yellow, red, and blue regions indicating the hydrophobic, ligand acceptor and ligand donor sites, respectively. (A) The identified sites are shown in dotted boxes. (B) the zoom-in view of most appropriate site.

Table 1.

SiteMap analysis on Mpro monomer.

Title	SiteScore	Size	Dscore	Volume (Å3)	Exposure	Enclosure	Contact	Phobic	Philic	Balance	don/acc
site_1	1.02	120	1.09	287.09	0.61	0.65	0.87	1.21	0.71	1.70	0.86
site_2	0.64	41	0.59	116.62	0.76	0.55	0.78	0.27	1.03	0.27	0.60
site_3	0.65	30	0.45	106.33	0.63	0.68	1.00	0.17	1.44	0.12	0.56
site_4	0.61	25	0.56	73.75	0.75	0.59	0.72	0.73	0.79	0.91	4.17

Table 2.

Hits selected from Proteases library along with their pharmacokinetic parameters and structures.

Name	Docking score	MM-GBSA ΔG Bind	aMol Weight (g/mol)	Target	Developmental phase	bHBA	cHBD	dRotatable bonds	ePSA (Å2)	flogP	glogS	Structure
N3	−10.6	−64.32	680.79	Mpro	Experimental	9	5	22	197.83	3.73	−4.89
13b	−6.656	−63.26	595.69	Mpro	Experimental	7	5	17	167.89	3.10	−3.90
Leupeptin Hemisulfate	−9.257	−80.784	426.554	Serine Protease	Experimental	5	6	18	166.27	0.58	−3.72
Pepstatin A	−9.919	−69.603	685.892	HIV Protease	Experimental	9	8	22	223	2.46	−4.19
Nelfinavir	−8.822	−68.943	567.8	HIV Protease	Approved	5	4	10	101.9	4.61	−5.5
Bortezomib (Velcade)	−8.291	−65.989	384.237	Proteasome	Approved	6	4	9	124.44	0.89	−3.9
Ixazomib (MLN2238)	−6.658	−59.781	361.029	Proteasome	Approved, Investigational	4	4	7	98.66	2.57	−4.5
MG-101 (ALLN) (calpain)	−7.126	−58.882	383.525	Cysteine Protease	Investigational	4	3	13	104	2.62	−3.07
Carfilzomib (PR-171)	−6.795	−58.155	719.91	Proteasome	Approved, Investigational	8	4	20	158.47	4.2	−5.2
L-685,458	−7.56	−57.905	672.853	Gamma-secretase	Experimental	6	5	22	159.85	4.76	−6.57
Calpeptin	−6.74	−57.113	362.463	Cysteine Protease	Investigational	4	2	12	84.5	3.23	−3.77
Z-FA-FMK	−6.829	−57.077	386.417	Cysteine Protease	Experimental	5	2	12	84.5	2.46	−3.3
Atazanavir	−6.853	−56.572	704.855	HIV Protease	Approved, Investigational	7	5	18	171.22	4.54	−5.3
ITF2357 (Givinostat)	−8.991	−56.502	475.965	HDAC	Investigational	5	3	9	90.9	3.51	−4.9
Indinavir	−7.946	−53.636	613.8	HIV Protease	Approved	7	4	12	118.03	2.81	−4.1
LAQ824 (Dacinostat)	−7.965	−53.51	379.452	HDAC	Phase I	4	4	10	88.59	2.52	−3.42
Anagliptin	−6.295	−53.377	383.447	DPP-4	Investigational	6	2	6	115.42	−0.54	−3.3
Aloin (Barbaloin)	−6.933	−52.231	418.394	Tyrosinase	Withdrawn	9	7	3	167.91	−0.14	−2.46
LY411575	−6.477	−52.07	479.475	Gamma-secretase	Experimental	6	3	5	98.7	2.99	−4.63
RG2833 (RGFP109)	−6.61	−51.173	339.431	HDAC	Experimental	3	3	8	84.2	3.04	−3.38
Ritonavir	−6.949	−50.393	720.944	HIV Protease	Approved	6	4	18	145.78	5.22	−5.8
E-64	−7.085	−46.011	357.405	Cysteine Protease	Experimental	6	5	11	172	−0.43	−0.8
GI254023X	−6.492	−43.484	391.504	Immunology and Inflammation Related	NF*	4	3	10	98.7	2.65	−3.41
Isorhamnetin	−6.901	−42.236	316.262	Tyrosinase inhibitor	Experimental	7	4	2	116	1.65	−3.36
TAPI-1	−7.768	−41.196	499.602	TNF-alpha protease inhibitor I	Experimental	6	6	13	163	1.75	−2.87
Trelagliptin	−6.936	−39.477	357.382	DPP-4	Investigational	5	1	3	93.67	1.3	−3.2
Amprenavir (Agenerase)	−6.61	−39.138	505.627	HIV Protease	Approved	6	3	11	131.19	2.43	−4
Lopinavir (ABT-378)	−7.607	−38.997	628.801	HIV Protease	Approved	5	4	15	120	4.69	−5.5
Delanzomib (CEP-18770)	−6.352	−37.975	413.275	Proteasome	Investigational	6	5	9	131.78	2.57	−3.8
Telaprivir	−6.78	−37.58	679.85	HCV proteases	Approved, Withdrawn	8	4	14	179.56	2.58	−4.3
Lycorine	−6.502	−35.836	287.31	HCV proteases	Experimental	5	3	0	62.2	0.8	−2.54
Darunavir	−6.9	−32.78	547.66	HIV Protease	Approved	7	3	11	140.42	2.82	−3.9

The general recommended ranges are as follows:

Molecular weight, <500.

Number of hydrogen bond acceptors, <10.

Number of hydrogen bond donors, <5.

Number of rotatable bonds, <10.

Polar surface area, <140 Å2.

Predicted octanol/water partition coefficient, −0.4 to +5.6.

Predicted aqueous solubility, < −5.0.

NF, not found.

Molecular dynamics simulation of Mpro in ligand bound and unbound systems

We performed MD simulations for both APO and COM proteins for stability of the systems with respect to time evolution. The RMSD (root mean square deviation) plot for APO-vs.-COM1 systems showed that APO protein is stable after 20 ns simulation while COM1 attains stability after 50 ns of simulation (Figure 3(A)). We observed a sudden deviation in the COM1 system that started around ∼35ns and lasted till ∼55ns. This increase in deviation is mainly attributed to extended loop conformation (residue 185 to 200), which is flexible in nature. The APO-vs-COM2 system showed that both the systems are stable throughout the trajectory with considerably minor deviations (Figure 3(B)). The radius of gyration (Rg) of the simulated systems also showed that all the systems have attained compactness (Figure 3(C)). To explore the fluctuations of the systems residue wise, the RMSF (root mean square fluctuation) analysis was carried out. From the RMSF plots we observed that binding of N3 made Mpro protein more flexible especially in the domain III as compared to 13b (Figure 3(D–E)). Since domain III is involved to form a homodimer so this might be the reason for fluctuation in this region in the presence of ligands. The ligand RMSD plots showed that they are stable during the dynamics. The atoms ranging from 1 to 8 and 40 to 49 of N3 molecules are observed to be highly fluctuated (RMSF > 1.0 Å) while all the atoms of 13b are observed with RMSF > 1.0 Å especially atoms 39–43 (Figure 4). This suggests that the flexible atoms of the ligands can be replaced by another atom (or group of atoms) which may restrict its movement in the pocket and ability to increase the potency that inhibit the Mpro with optimum efficiency. In case of COM2, the highest fluctuation is observed for E290 that is responsible for dimer formation of the Mpro by salt-bridge interactions between E290 of one monomer and R4 of the other (Zhang et al., 2020). The binding site RMSD of systems APO, COM1 and COM2 are very stable which signifies that residues are undergoing least conformational changes (Figure 3(F)). Also, in comparison to APO, inhibitor N3 (COM1) has induced more stability in the binding site as compared to inhibitor 13b (COM2). However, it can be noted that there is slight deviation in RMSD of COM1 and COM2 by the convergence of simulation (Figure 3(F)). From the MD analysis it is very well evident that the binding site is well compact and the most stable state was extracted for the VS purpose.

Figure 3.

Changes in the Mpro structures (APO-vs.-COM) and its dynamics with respect to time. (A) The RMSD was calculated throughout the MD trajectory simulation time of 200 ns using backbone atoms of APO (black) with (A) COM1 (red), (B) COM2 (blue), respectively, (C) Radius of gyration of Mpro simulated systems. The RMSF values for the APO with (D) COM1 and (E) COM2 were plotted using C-alpha atoms. The domains and the respective differences between COMs-vs.-APO has been shown in terms of difference graph in brown (panel D) and in green (panel E). (E) RMSD and (F) Solvent accessible surface area (SASA) of binding sites of N3 (red) and 13b (blue) compared with APO (black) over 200 ns.

Figure 4.

The RMSD and RMSF plots for the control ligands N3 and 13b are shown with respect to their MD trajectory.

The analysis of stable protein-ligand interactions and the contribution of a particular residue towards the ligand binding in the pocket is an important aspect to identify the hot-spots. From the MD simulation of COM1, the residues T26, G143, C145, Q189 and T190 are observed to be involved in HB formation for more than 60% (occupancy) of the simulation time (Figure 5(A)). The residue H41 is observed to form water mediated interaction with the ligand for 56% of the simulation time. The residues N119, S144, Q189, T190, and Q192 also form transient water mediated interactions. In the COM2, the residues E166 and T190 are highly involved in the HB stability for more than 60% of their occupancy (Figure 5(B)). Residue E166 is also involved in transient water mediated and salt bridge interactions. The residues H41, C145, H164, and Q189 also establish water mediated interactions for certain time during the simulation time. The hydrophobic contributions by residues H41, M49 and P168 were observed, however its occupancy is less than 40% in the simulation time. Our MD simulations shows that the pocket is hydrophilic in nature but there is no stable water molecule involved during the dynamics. The binding site SASA (solvent accessible surface area) shows that APO protein is quite more solvent exposed than both the complexes (Figure 3(G)). The residues C145, Q189, and T190 are the common key residues that are observed from the protein-ligand interactions occupancy plots.

Figure 5.

A schematic representation of ligand–atom interactions with the protein residues that occur more than 30% of the simulation time in the trajectory of 200 ns is shown for both (A) N3 and (B) 13b molecules. Protein–ligand interactions monitored throughout the simulation are represented as the bar plot. The interactions are categorized as HBs, hydrophobic, ionic and water bridges. The stacked bars are normalized over the course of the trajectory: for example, a value of 0.6 suggests that the specific interaction is maintained for 60% of the simulation. Values over 1.0 or 100% are possible as some protein residue may make multiple contacts of same subtype with the ligand.

From the overall analysis of MD simulation, we conclude two main aspects, 1) we obtained the stable computational model (COM1) as it is showing the intact binding site with minor conformational changes and, 2) the insights of the binding pocket, where hot-spot residues and their respective interactions were elucidated qualitatively and quantitatively. Both aspects were applied for VS of the library.

Virtual screening of two different libraries

Protease library on the Site1 of Mpro protein

It has been a proven and successful strategy to inhibit the viral proteases for the treatment of viral infections such as in the cases of human immunodeficiency virus (HIV) and hepatitis C virus (HCV) and hence screening of protease inhibitors could be a useful approach against COVID-19 Mpro (Ghosh et al., 2016; Yang et al., 2006). With this logic, 227 protease molecules from Selleckchem were curated and prepared for screening on Site1. Based on our filtering criteria of optimum MM-GBSA and docking scores, the top 30 initial hits were filtered out (Figure 6(A)). The docking scores, MM-GBSA, pharmacokinetic descriptors, their known targets, along with the structures of the hits are shown in Table 2. Other than the quantitative parameters (docking score and MM-GBSA energy), the additional qualitative parameters were chosen to screen the docked conformations which show similar interaction patterns and interactions with respect to controls. However, some of these molecules (including our control molecules: N3 and 13b) violate the few pharmacokinetic properties but they can be used as immediate starting point for the experimental validation and optimization as per their potency and efficacy that can facilitate the identification of more potent drug-like molecules. For all selected molecules (top hits), the residue mapping based on interaction pattern was executed (Figure 7). The top three molecules found in our study are Leupeptin Hemisulfate, Pepstatin A and Nelfinavir on the basis of MM-GBSA and docking scores and their ligand interaction diagrams are shown in Figure 8(A–C) respectively. However, few of the molecules are observed to be common in our screening and some in-silico studies published targeting Mpro and they are Nelfinavir, Lopinavir, Indinavir, Ritonavir, Darunavir (Das et al., 2020; Khan, Zia, et al., 2020; Muralidharan et al., 2020; Nutho et al., 2020) highlighted in cyan colour (Table 2).

Figure 6.

Scatter plot to calculate the docking scores (x-axis) and predicted MM-GBSA binding free energies (kcal/mol) (y-axis) for hits obtained from VS of (A) protease library (B) FDA/repurposed molecules library. The points are coloured by MM-GBSA ΔG values. In plot (A) the total conformers obtained after docking from 227 molecules protease library: 399 and in plot, (B) total conformers reported after virtual screening from 13947 FDA/repurposed library: 531.

Figure 7.

Interaction fingerprint: Protein-ligand interaction fingerprinting for hits obtained from protease library at site 1 residues involved in the interaction site are within 4.0 Å along with control molecules. The presence and absence of residues are marked as 1 and 0 respectively. The residues are coloured on the basis of polar (blue), hydrophobic (green), acidic (red) and basic (purple) properties. The residues which are showing > 90% occupancy are considered as key residues. The occupancy was also calculated and mentioned at the bottom line.

Figure 8.

Interaction map of top hits: (A–C) The docked poses of Leupeptin Hemisulfate, Pepstatin A, Nelfinavir respectively within the binding site. Residues lining the pocket (Cyan) under 4.0 Å and its respective inhibitors (purple) are shown in licorice representation. Red and green dotted line shows the HB and pi-pi interactions and red solid line means salt bridge.

FDA/repurposing library on the Site1 of Mpro protein

Another library of FDA/repurposing molecule curated from different sources such as repurposing hub, Selleckchem and DrugBank were used to prepare for screening purpose. Due to high failure rates, considerable costs and slow pace of new drug discovery and development, re-use of ‘old’ drugs to treat diseases becomes the quickest route to find the active molecules, as it requires the use of known bioactive moleculeswith potentially lower production costs and shorter timelines (Pushpakom et al., 2019). From the multi-step route of VS, we filtered out the best 41 hits (Figure 6(B)). The docking scores, MM-GBSA, pharmacokinetic descriptors, their known targets, along with the structures of the hits are shown in Table 3 and their residue wise interaction mapping is also carried out (Figure 9). From this library, the best three molecules are Birinapant, Lypressin and Octreotide and their ligand interaction diagrams are shown in Figure 10(A–C) respectively. Nevertheless, these molecules violate the pharmacokinetic properties but could be optimised to drug-like properties based on the experimental results. The finding of some common top molecules in this study are also reported in other screening studies on the same target, claiming the robustness of our protocols. These common molecules are Lypressin Octreotide, Mitoxantrone, Hesperidin, Echanoside, Pralmorelin, Epicatechin, Diosmin, Flavitan (flavin adenine dinucleotide), Curcumin, Saquinavir, Montelukast, Baicalin, Thymopentin (Das et al., 2020) highlighted in cyan colour (Table 3).

Table 3.

Hits selected from FDA/Repurposed library along with their pharmacokinetic parameters and structures.

Title	Docking score	MM-GBSA ΔG Bind	Phase	Source	aMW	bRB	cPSA (Å2)	dHBD	eHBA	flogP	glogS	Pharmacological activity	Structures
Birinapant	−8.141	−105.15	Approved	Drugbank	806.9	15	194.2	8	10	2.1	−4.6	Peptidomimetic activator of SMAC and inhibitor of IAP; potential antineoplastic activity.
Lypressin	−7.859	−102.499	Approved	Drugbank	1056.2	19	452.9	16	24	2.2	−4.3	Antidiuretic hormone
Octreotide	−7.202	−94.415	Approved	Drugbank	1019.2	17	332.3	13	12	0.4	−4.9	Potent inhibitor of growth hormone, glucagon, and insulin
PD-173212	−7.168	−93.65	Pre-clinical	Repurp hub	599.85	15	65.72	1	7	7.7	−9.08	Calcium channel blocker
Mitoxantrone	−8.232	−93.159	Approved	Drugbank	444.5	12	163.1	8	10	0.9	−2.8	Antineoplastic activity
Dipyridamole (Persantine)	−7.184	−91.983	Approved	Drugbank	504.6	12	145.4	4	12	1.5	−2.7	Nucleoside transport inhibitor and a PDE3 inhibitor; inhibits blood clot formation
Amikacin	−8.187	−88.745	Approved	Drugbank	585.607	10	331.9	13	17	−8.6	−1.1	Aminoglycoside antibacterial agent
Hesperidin	−12.344	−88.5	Approved	Drugbank	610.5	7	234.2	8	15	−0.2	−2.4	Neurological conditions, antioxidant and anti-inflammatory effects
Echinacoside	−11.473	−87.149	Investigational	Drugbank	786.7	14	324.4	12	19	−0.9	−2.3	Treatment of Neurological and other Disorders
Ginsenoside Rg1	−10.349	−86.87	Approved	Repurposing hub	444.7	4	40.4	2	2	6.3	−6.4	Antitumor, anti-inflammatory, antioxidation, and inhibition of cell apoptosis.
Pralmorelin (GHRP-2)	−9.344	−86.865	Phase II	Repurposing hub	818	21	256	9	8	2.4	−5	Acts to endogenously increase growth hormone release from the pituitary.
PSI-7976	−8.19	−86.73	Pre-clinical	Repurp hub	529.45	10	170.845	3	14	1.07	−3.71	HCV inhibitor
Boceprevir	−7.028	−86.398	Approved, withdrawn	Drugbank	519.6	10	150.7	4	5	1.9	−4.4	HCV protease inhibitor (genotype 1).
Acteoside	−10.631	−85.38	Investigational	Drugbank	624.5	11	245.2	9	14	1	−2.8	A neuroprotective agent, an antileishmanial agent, an anti-inflammatory agent, a plant metabolite and an antibacterial agent.
TMC353121	−11.266	−83.08	Approved	Repurposing hub	558.727	13	107.7	4	6	4.09	−5.45	Potent respiratory syncytial virus (RSV) fusion inhibitor
TXA127	−9.5	−82.253	Phase II	Repurp hub	899.02	28	421.153	10	20	−2.93	−2.12	Angiotensin receptor agonist
Evista (Raloxifene Hydrochloride)	−7.086	−81.815	Approved, investigational	Drugbank	473.5	7	70	2	5	5.4	−6	Treatment of osteoporosis in postmenopausal women and those on glucocorticoids
Sunitinib	−7.381	−80.705	Approved	Drugbank	398.4	7	77.2	3	3	3.2	−4.1	Used to treat certain types of advanced or progressive tumors of the stomach, intestines, esophagus, pancreas, or kidneys.
Rutin	−13.104	−76.635	Approved	Drugbank	610.5	6	267	11	17	1.7	−3.1	Radioprotective and antiplatelet activity. Rutin reduces hepatic and blood cholesterol levels.
MM-102	−10.424	−76.45	Pre-clinical	Repurposing hub	669.8	21	178.3	7	7	2.85	−5.49	Peptidomimetic MLL1 inhibitor
Felypressin	−8.486	−76.45	Pre-clinical	Drugbank	669.8	19	264.7	12	13	−1.1	−4.4	Non-catecholamine vasoconstrictor
Epicatechin	−11.036	−74.878	Investigational	Drugbank	290.2	1	110.2	5	6	1	−2.6	Prevent the onset of type II diabetes and many cardiovascular diseases
Quercetin 3-Rutinoside	−11.362	−71.34	Approved	Drugbank	610.5	6	265.6	10	16	0.1	−2.2	Antiallergic, anti-inflammatory, antiproliferative, and anticarcinogenic properties
Lanreotide	−8.097	−69.96	Approved	Drugbank	1096.3	17	355	13	12	1.87	−5.3	Management of acromegaly and symptoms caused by neuroendocrine tumors,
Diosmin	−10.032	−67.23	Approved, investigational	Drugbank	608.5	7	234.2	8	15	−0.4	−2.6	Venous disease
KD025 (Slx-2119)	−7.023	−66.999	Phase II	Repurposing hub	608.5	8	104.82	3	5	2.95	−5.7	Use idiopathic pulmonary fibrosis (IPF), selective ROCK2 inhibitor.
Cefoperazone	−7.074	−66.95	Approved, investigational	Drugbank	645.6	9	220.2	4	11	−0.1	−3.4	Effective against Pseudomonas infection and various bacterial infections.
Pilaralisib	−7.691	−66.86	Investigational	Drugbank	541	7	148.3	4	8	3.98	−5	Potential antineoplastic activity
Troxerutin	−10.993	−66.24	Investigational	Drugbank	742.6	15	293.2	10	19	−0.5	−2.7	Chronic venous insufficiency; vasoprotective
Vitexin	−7.858	−65.995	Investigational	Drugbank	432.4	3	177	7	10	1.38	−2.8	Platelet aggregation inhibitor, an antineoplastic agent and a plant metabolite.
Tenofovir Alafenamide	−7.379	−65.61	Approved	Drugbank	476.4	12	143.8	2	7	1.49	−3.3	Treat chronic hepatitis B, and prevent HIV-1 infections
Argatroban Monohydrate	−7.932	−64.8	Approved, investigational	Drugbank	508.6	8	177.7	6	9	−0.9	−3.5	Synthetic direct thrombin inhibitor derived from L-arginine
Flavitan	−10.237	−62.71	Approved	Drugbank	785.5	13	356	9	20	−2.7	−3.2	Vitamin B2 deficiency or metabolic disorder including stomatitis, eczema, etc
Curcumin	−8.941	−61.71	Approved, Experimental, Investigational	Drugbank	368.3	8	93	2	6	3	−4.3	Chemopreventive and anticancer agents.
Saquinavir	−7.461	−56.66	Approved, Investigational	Drugbank	670.8	13	166.7	5	7	4	−5.4	Anti-HIV protease inhibitor
Montelukast	−7.418	−53.21	Approved	Drugbank	586.1	12	70.4	2	4	7.2	−7.8	Leukotriene antagonist, an anti-asthmatic drug and an anti-arrhythmia drug
Dapagliflozin	−8.205	−50.56	Approved	Drugbank	408.8	6	99.3	4	6	2.5	−3.4	odium-glucose cotransporter 2 inhibitor
Baicalein	−8.169	−47.81	Pre-clinical	Repurposing hub	270.2	1	87	3	5	1.7	−3.3	Anti-inflammatory, antioxidant, antiviral, and antitumor effects, lipoxygenase inhibitor
Pitavastatin	−8.435	−47.71	Approved	Drugbank	421	8	90.6	3	5	3.75	−5	Novel statin that induces plaque regression and elevates HDL-cholesterol levels
Thymopentin	−8.026	−40.28	Investigational	Drugbank	679.9	22	327.6	11	14	−3.3	−4	Treatment of rheumatoid arthritis, AIDS, and other primary immunodeficiencies.

The general recommended ranges are as follows:

Molecular weight, <500.

Rotatable bonds <10.

Polar surface area, <140 Å2.

Number of hydrogen bond donors, <5.

Number of hydrogen bond acceptors, <10.

Predicted octanol/water partition coefficient, −0.4 to +5.6.

Predicted aqueous solubility, < −5.0.

NF, not found.

Blue highlighted rows are the molecules found common in screening results of other papers.

Figure 9.

Interaction fingerprint: Protein-ligand interaction fingerprinting for hits obtained from FDA and Repurposing library along with our controls at site 1 residue involved in the interaction site are within 4.0 Å. The presence and absence of residue are marked as 1 and 0 respectively. The residues are coloured on the basis of polar (blue), hydrophobic (green), acidic (red) and basic (purple) properties. The residues which are showing > 90% occupancy are considered as key residues. The occupancy of each residues is calculated at the bottom line.

Figure 10.

Interaction map of top hits: (A–C) The docked poses of Birinapant, Octreotide and Lypressin, respectively within the binding site. Residues lining the pocket (cyan) under 4.0 Å and its respective inhibitors (purple) are shown in licorice representation. Red dotted line and red solid line shows the HBs and salt bridge interactions.

Characterization of binding site

The residues that have shown the highest occupancy (>90%) are re-analysed. The residues H41, M49, N142, C145, H164, M165, E166, D187, R188 and Q189 for all the hits obtained from protease library (Figure 7) with respect to the control compounds N3 and 13b. The residues G143, E166, and Q189 are observed to be majorly involved in HBs interactions for most of the hits while H41 is observed to form pi-pi stacking and pi-cation interaction in some cases (Figure 7). Similarly, for the hits obtained from FDA/repurpose library, the residue interaction mapping showed that the residues H41, M49, N142, C145, M165, E166, R188, and Q189 have higher occupancy (>90%) in the binding site (Figure 9). The HBs analysis found that residues T26, G143, E166, and Q189 are highly involved in the HB interactions (Figure 9). From the residue interaction mapping analysis of hits, we observed that the binding site of Mpro is composed of 44.7% of polar residues, 42.1% of non-polar residues, 10.5% of acidic residues and 2.6% of basic residues which encompasses a diverse class of molecules. This indicates that the binding site (site1) of Mpro is hydrophilic in nature and solvent exposed which is in concordance with our site map results. However, the flexibility and adaptability of the pocket towards the ligand explored by molecular dynamics simulation matches nicely with the previous findings. Also, we find H41, M49, N142, C145, H164, M165, E166 and Q189 are the most conserved residues and their high occupancy in all our top hits. As reported by Wang (2020) that MERS, SARS and COVID-19 have four (H41, H163, M165 and Q189) common residues, while HCV NS3/4A and COVID-19 have only one common hotspot residue (H41), confirming our findings that the hits having interaction with these residues might have broad spectrum value.

Analysis of top ranked molecules through MD simulations

We identified six molecules with a lower free energy of binding combined with a higher theoretical drug discovery value compared to the co-crystallized ligand N3 and 13b (Figure 11). The MD analysis of top molecules highlight that Nelfinavir and Birinapant have granted maximum stability to protein throughout simulation as compared to other four selected hits (Supplementary Figure 2). However, the other four hits were also found to be stable which is well reflected in the result (Supplementary Figure 2). The MM-GBSA calculations from the equilibrated MD trajectories of the hits and control molecules showed that all the hits (except for Leupeptin Hemisulfate) has shown considerably better binding energy than the control molecules (Supplementary Table 1). Some molecules have shown a minor decrease in their binding energies, also some of them have presented a substantial improvement in comparison to the co-crystallized ligand. We observed that their docking scores ranging from −7.0 to −9.0 kcal/mol, however, some of our molecules have shown substantially good docking energy as well as their free energy is also considerably higher side which is consistent and reflected well through MD simulations. Our identified hits interacted with the protease with at least two HBs with an average of over four (Figures 8 and 10), appearing as the hot-spot residues whose knowledge is critical for structure-based hits identification.

Figure 11.

A plot of the MM-GBSA values (primary y-axis) and docking score (secondary y-axis) for the control molecules and 6 finalized potential hits along with their 2D structures.

Furthermore, the identified molecules have shown relevant pharmacokinetic descriptors with logP values, MW, PSA and HB donors/acceptors in the range described by Lipinski. Though there are many approved and successful drugs which are documented of not following the standard Lipinski rule, especially for protein-protein interactions disruptors. It is also encouraging that the identified hits are available at commercial suppliers such as MCULE and Selleckchem. The identified top six molecules presented comparatively better docking and post-processing MD simulation analysis of free energy i.e. MM-GBSA scores with respect to the control molecules and therefore represent excellent candidates for further investigation in-vitro (Figure 11). From the overall analysis, the identified set of molecules can be considered as early hit molecules, albeit we accept that no experimentally supported hit-optimization was conducted, though these molecules looks promising.

Conclusion

As the cases increase day by day there is an extremely urgent need of the designing small compound or peptide drugs to cure the 2019-nCoV. Based on crystal structures and MD simulations a computational model of Mpro was chosen for VS. It is an important step before performing VS, as it not only highlights the dynamical changes in the pocket which should be taken in account for reducing false-positives results, but also highlight the overall impact of binding of the compounds on the protein. Thereafter, interaction fingerprinting of binding sitesof static (from crystals)-vs.-dynamic (from MD simulations) changes, a structure-based VS of known and existing molecule libraries were performed to identify the potential molecules that possibly repurposed against Mpro. A high-resolution crystal structure of COVID-19 Mpro in complex with N3 and knowledge of key residues from 13b molecules is available on time, allowing us to conduct drug repurposing. Through MD simulations we have observed and quantified the dynamical changes and nature of the binding site at residue level. The investigation of solvation sites or the role of water molecules in the binding site could be carried out to explore the druggability of this pocket that might help in rational designing of molecules. The comparative analysis of identified hits with N3 and 13b further supports our findings. The knowledge of interactions in terms of HBs, hydrophobic contacts, salt-bridges, and pi-pi stacking and their conformations was considered exclusively for filtering of the molecules, other than quantitative values such as Dscore, docking and ΔG of the complexes.

Currently, antineoplastic, immunomodulators, nucleotide inhibitors, antimalarial, ribonucleoside inhibitors, steroid hormones, protease inhibitors, antiretrovirals are being evaluated in clinical trials against COVID-19. We identified similar category of molecules in our VS approach. The hits found in our study belongs to different chemical classes and have the potential to accommodate inside the pocket that can be further quickly explored by antiviral experimental assays for binding affinity and inhibition activity. There will be no prior need for synthesis of these molecules as they are easily available and these predicted hits or their scaffolds can be a good starting point for developing novel hits if confirmed in in-vitro studies. The potential hits listed in this study are promising candidates and can facilitate the hunt of anti-COVID-19 Mpro drug discovery.