Comparative molecular investigation of the potential inhibitors against SARS-CoV-2 main protease: a molecular docking study.

Recent outbreak of novel coronavirus and its rapid pandemic escalation in all over the world has drawn the attention to urgent need for effective drug development. However, due to prolonged vaccine and drug development procedure against a newly emerged devastating SARS-CoV-2 virus pathogen, repurposing of existing potential pertinent drug molecules would be preferable strategy to reduce mortality immediately and further development of new drugs to combat overall global Covid-19 crisis in all over the world. Herein, we have filtered 23 prospective drug candidates through literature review. Assessing evidences from molecular docking studies, it was clearly seen that, Epirubicin, Vapreotida, and Saquinavir exhibited better binding affinity against SARS-CoV-2 Main Protease than other drug molecules among the 23 potential inhibitors. However, 50 ns molecular dynamics simulation indicated the less mobile nature of the docked complex maintaining structural integrity. Our overall prediction findings indicate that Epirubicin, Vapreotida, and Saquinavir may inhibit COVID-19 by synergistic interactions in the active cavity and those results can pave the way in drug discovery although it has to be further validated by in-vitro and in-vivo investigations.Communicated by Ramaswamy H. Sarma.

Introduction

Following detection of SARS-CoV-2 in late December 2019 in Wuhan, capital of Central China's Hubei Province, it has quickly spread across the globe (Tang et al., 2020) and this highly spreading and devastating infection, COVID-19 has killed 3,70,120 people worldwide (https://www.worldometers.info/coronavirus/#countries, n.d.). Named after the related structured severe acute respiratory syndrome-related coronaviruses, SARS-CoV-2 has a genome sized∼30kb with the specific pattern, ORF1ab encoding the large replicase polyprotein pp1ab and forming the non-structural proteins (nsp1-16) as well as the structural proteins (S, E, M, and N) (Khailany et al., 2020). Moreover, the SARS-CoV-2 produced 6 accessory proteins, encoded by ORF3a, ORF6, ORF7a, ORF7b, ORF8, and ORF10 genes (Khailany et al., 2020). Two cysteine proteases, namely the main protease (Mpro) or the 3C-like protease (3CLpro) and the secondary papain-like protease 2 (PL2pro) mainly process the replicase complex (RNA-dependent RNA polymerase and the helicase) that ultimately facilitate the viral replication, where Mpro and PL2pro cleaves at 11 and three sites of the polyprotein, respectively. Wu et al. (2020) reported that among the non-structural proteins of SARS-CoV-2, 3CLpro plays a key role in the replication and maturation step.

Based on these facts, the Mpro is a more promising drug target against SARS-CoV-2 than the PL2pro, as well as for potential drug targets; scientists generally choose proteins and enzymes of CoVs vital for the replication process. In addition, previous studies have also shown that the main proteases of different coronaviruses are highly conserved in terms of both sequences and 3D structures (Liu & Wang, 2020). These features, together with its functional importance, have made 3CLpro an attractive target for the design of anti-corona viral drugs (Liu & Wang, 2020) than other corona viral proteins.

Until now, no clinically approved drug is available for coronavirus (SARS-CoV-2) treatment (Battegay et al., 2020). Scientists are trying to solve the puzzle through drug repurposing techniques, which is a promising option considering the global pandemic caused by the virus (Pushpakom et al., 2019). Drug repurposing, an alternative option for approved or investigational drugs outside their defined indication, has a lower risk of failure and drastically reduces the time frame for development, facilitating prompt clinical decisions at lower costs development (Cherian et al., 2020).

In the recent reports, we have come across some of the scientific endeavors, most of which are related to drug repurposing, to find effective drugs, nucleotide analogs (Chien et al., 2020; Kirchdoerfer, 2020), and peptide binder (Zhang et al., 2020). For example, Ritonavir and Lopinavir, which were under phase II/III clinical trials for MERS-CoV (Chan et al., 2015; Kim, Won, et al., 2016), have activity against HCoV-229E, HCoV-NL63 and animal CoVs (Kim, Won, et al., 2016), and can be an excellent repurposed drug for SARS-CoV-2 because the drug binding sites and 96% protein sequences were similar in SARS-CoV and SARS-CoV-2 main protease (Liu & Wang, 2020). Also, researchers have currently reported using the drug repurposing approach based on the molecular docking and dynamics study where the key target proteins are 3CL protease, RNA dependent RNA polymerase (RdRp), and spike proteins (Elfiky, 2020b; Muralidharan et al., 2020; Smith & Smith, 2020; Tahir Ul Qamar et al., 2020; Yu et al., 2020).

In this study, we explore potential drug candidate(s) against SARS-CoV-2 by performing the virtual screening based on drug repurposing concept, using available commercial medicine (23 potential drugs against other viruses) that can eventually, be predicted to act against the main protease (Mpro), or 3-C like protease (3CLpro). Due to the lack of confidence in giving the accurate binding energy, binding mode, solvation effect, constant motion effect, and entropy energy during molecular docking (Sethi et al., 2019), a molecular dynamics study was preferred.

Materials and method

Molecular docking preparation

In this study, the crystal structure of the main protease of SARS-CoV-2 was extracted (PDB ID: 6LU7) which has similar binding sites of SARS-CoV (Rose et al., 2016). The protein structure was cleaned and prepared by Discovery Studio (San Diego: Accelrys Software Inc, 2012) and Pymol (DeLano, 2002) software package. The cleaned protein structure was prepared and minimized with GROMOS 43B1 force field with the aid of SWISS PDB Viewer (Kaplan & Littlejohn, 2001). The ligand molecules were retrieved from PubChem database (Kim, Thiessen, et al., 2016) after a deep literature review process. The ligand structure was optimized and cleaned with the help of mmff94 force field (Halgren, 1996) with the steepest gradient algorithm.

Molecular docking study

Molecular docking was performed to understand the binding affinity with the ligand and drug compounds. The ligand were converted in PDBQT format to make them acceptable format for employing docking in AutoDock Vina (Trott & Olson, 2010). The center of the grid box was X:-26.289, Y: 12.606, Z: 58.9496 and dimensions were (X: 50.8171, Y: 67.5389, Z: 59.6721) Ǻ. Finally the docked pose were analyzed with the aid of Discovery studio and Pymol.

Molecular dynamics simulation

The molecular dynamics simulation was run in YASARA application (Krieger et al., 2004). The AMBER14 force field (Case et al., 2007) was employed for the minimization of the drug-protein complex and transferable intermolecular interaction potential 3 points (TIP3P) were used for the addition of Na and Cl ions. The box size of the cell was set as (96.9654 × 96.9654 × 96.9654) Ǻ with a periodic boundary condition. A cut off the radius of 8 angstroms was used for the calculation of short-range Van der Waals and Coulomb interaction. The Particle Mesh Ewalds was used to calculate long-range electrostatic interaction with a physiological condition at 300 K, pH 7.4, 0.9% NaCl (Krieger et al., 2006; Krieger & Vriend, 2015). The simulation was run for 50 ns and root mean square deviation (RMSD), root mean square fluctuation (RMSF), the radius of gyration (Rg), solvent-accessible surface area (SASA). Moreover, all the trajectory further used to calculate MM-Posson-Boltzman surface area or MM-PBSA binding energy for binding free energy calculation where positive energy denotes favorable binding (Dash et al., 2019; Krieger & Vriend, 2015).

Result

Among 23 drug molecules only three showed better binding energy while interacting with Main Protease (Table 1).

Table 1.

List of potential drug compounds and their molecular docking score energy.

Compound ID	Compound Name	Docking Energy (kcal·mol−1)
CID-4748	Perphenazine	−6.1
CID-37542	Ribavirin	−6.3
CID-41867	Epirubicin	−9.5
CID-64143	Nelfinavir	−7.1
CID-65016	Amprenavir	−7.1
CID-71306	Vapreotida	−9.1
CID-92727	Lopinavir	−9.4
CID-131536	Fosamprenavir	−6.6
CID-148192	Atazanavir	−7
CID-164522	Bepotastine	−7.1
CID-213029	Darunavir	−7.1
CID-392622	Ritonavir	−7.1
CID-441243	Saquinavir	−9.5
CID-454216	Valrubicin	−7.8
CID-492405	Favipiravir	−5.1
CID-2826718	Caspofungin	−7.1
CID-5311054	Collistin	−6.6
CID-5362440	Indinavir	−7.8
CID-6918173	Icatibant	−7.1
CID-10445549	Galidesivir	−5.9
CID-54682461	Tipranavir	−7.6
CID-121304016	Remdesivir	−7.2
CID-135413536	Aprepitant	−7.3

Epirubicin, formed five hydrogen bond at Phe140, Glu166, Gln189, Asn142, Pro168 and one hydrophobic bond at Pro168 position. On the other hand, Vapreotida and main protease complex stabilized by several hydrogen bond at Gln110, Asp245, Arg105, Asp245, Ile249 and one pi-pi-stacked bond Phe294, one pi-pi-T shaped at His246, one alkyl Pro252 and two pi-alkyl bond at Val202, Ile249. Besides, main protease and Saquinavir complex stabilized by six hydrogen bond at Ile249, Phe294, Gln110, Thr292, Pro293, Phe294 and six hydrophobic bond at Asn245, Phe294, Ile249, Val202, Pro293, Ile249 (Table 2 and Figure 1).

Table 2.

Docking simulation results with docking score energy and interaction with amino acids.

Compound Name	Docking energy (kcal·mol−1)	Hydrogen Bond	Hydrophobic Bond	MM-PBSA (kj/mol)
		Ligand atom-Amino Acid Distance(Å)	Interaction-Amino Acid Distance(Å)
Epirubicin	−9.5	H-Phe140:O(1.88),H-Glu166:OE2(2.68),HE21-Gln189:O(2.12),HA-Asn142:O(2.65),HA-Pro168:O(2.68),	Pi-Sigma-Pro168(4.70)	93.405 ± 1.1
Vapreotida	−9.1	H-Gln110:OE(2.08),H-Asp245:OD1(2.98),H-Arg105:O(2.54),H54-Asp245:OD1(2.96),HA-ILE249:O(2.52),	Pi-Pi-Stacked Phe294(3.79),Pi-Pi-Tshaped-His246(4.79),Alkyl Pro252(4.47),Pi-Alkyl Val202(5.02),Pi-Alkyl Ile249(4.58)	53.042 ± 0.9
Saquinavir	−9.5	H-Ile249:O(2.44),HN-Phe294:O(2.45),H21-Gln110:OE(2.38),HB-Thr292:O(2.49),HD1-Pro293:O(2.52),HA-Phe294:O(3.01),	Pi-anion Asn245(4.56),Pi-Pi-Stacked Phe294 (3.87),Alkyl Ile249(4.45),Alkyl Val202(4.87),Alkyl Pro293(5.28),Pi-Alkyl Ile249(5.01)	73.807 ± 1.0

Figure 1.

Molecular docking interactions between the SARS-CoV-1 main protease and drug compounds. A denoted Epirubicin, B denoted Vapreotida, C denoted Saquinavir.

Molecular dynamics simulation study

The RMSD value of docked complex analyzed from simulated trajectory to find out whether the complex remain stable or complex rigidity. From the Figure 2, it can be observed that, all three chemical drugs had lower pic than 2.5 angstrom during whole simulation time

Figure 2.

Root Mean Square Deviation (RMSD) of the C-alpha atom, (B) Root Mean Square Fluctuation (RMSF) of the drug-protein complex.

However, Epirubicin showed less flexibility till 10 ns and thereafter, fluctuated for the next 5 ns. The Epirubicin and main protease protein from SARS-CoV-2 were stabilized and then reached in a stabilized state. On the other hand, Vapreotida and Saquinavir had similar profile till 10 ns and after then, it exhibited lower RMSD value. It can be denoted from Figure 2, that all the complex stabilized despite have some minor fluctuation but they were way below 2.5 angstrom. On the other hand, root means square fluctuation of the drug-protein complex can help to understand the flexibility in the protein region. From Figure 2 it can be understood that all three drug molecules exhibit a similar level of flexibility over the protein region. Most of the residue did not have a higher peak in RMSF therefore responsible for the complex overall stability. However, Figure 2 indicates that the amino acid; Ser1 (helix), Gly2 (helix), Arg4 (helix), Glu47 (beta turn), Leu50 (beta turn), Glu55 (helix), Arg60 (helix), His64 (helix), Asn72 (beta hairpin), Arg76 (helix), Lys97 (beta turn), Tyr154 (beta hairpin), Pro168 (beta hairpin), Thr169 (beta hairpin), Arg222 (beta turn), Phe223 (beta turn), Thr224 (beta turn), Met235 (helix), Lys236 (helix), Gln244 (helix), Asn277 (beta turn), Gly278 (beta turn), Arg279 (beta turn), Phe294 (helix), Arg298 (helix), Gln299 (helix), Ser301(beta turn), Gly302 (beta turn), Val303 (beta turn), Thr304 (beta turn), Phe305 (beta turn), Gln306 (beta turn) from all the three complex were most flexible than other residue in the system.

The complex compactness and rigidity can be assessed through a radius of gyration. Figure 3 denoted that, Vapreotida and main protease has lower Rg value than Saquinavir and Epirubicin which indicates the degree of compactness of the Vapreotida.

Figure 3.

(A) Radius of Gyration of the complexes was explored to find out the compact nature of the system, (B) Solvent Accessible Surface Area (SASA) was determined to understand the changes in surface area.

Moreover, Epirubicin and protein complex is less compact than the other two docked complex as they had higher Rg value among three complex.

On the other hand, Higher SASA value indicates the high degree of expansion in the protein surface area. The Saquinavir had a slightly higher peak than the other two drug molecules till 10 ns. Subsequently, Saquinavir, Vapreotida, Epirubicin had a similar trend to 30 ns which indicates protein surface area remain same during this simulation time. However, Epirubicin had a lower SASA peak than two complexes after 30 ns which may be the result of complex area expansion. Both Saquinavir and Vapreotida had similar SASA profiles for the rest of the simulation time.

Furthermore, we have calculated MM-PBSA binding energy by taking each trajectory from the simulated protein complex to understand the relation between ligand binding and structural changes (Figure 4). The main protease complex and Epirubicin exhibited lower binding energy than other two complexes. The rest two complex had higher binding energy which indicates better binding with main protease of SARS-CoV-2.

Figure 4.

Molecular mechanics Poisson-Boltzmann surface area calculation of the screened drug complex system.

Discussion

Molecular modelling especially docking and dynamic simulation study has become essential part in modern drug discovery pipeline(Bappy et al., 2020; Elfiky, 2020a; Islam et al., 2019). Moreover, time and lab cost can significantly reduce by narrowing the drug molecule choices in in silico based drug discovery process (Talele et al., 2010). Besides, this combinatorial computational approaches allow the researcher to get the binding insights of the drug molecules to the target protein and rationally design the therapeutic targets (Mullard, 2019).

Among the selected drugs from this study, Epirubicin can have a negative impact on chronic HBV patient through reactivation of the viral genome (Xu et al., 2014). The Epirubicin drug molecules showed −9.5 kCal/mol energy while binding with the main protease of SARS-CoV-2 which is the highest binding energy among the drug molecules. This might be the result of several hydrogen bonding that was found such as Phe140, Glu166, Gln189, Asn142, and Pro168 at the active cavity of the main protease enzyme. Researcher suggests that non-covalent binding at the active site of the target protein may lead to possible inhibition and blockage of the protein (Hasan et al., 2015). On the other hand, the somatostatin analog, Vapreotide (VAP) was used to control tumor growth, progression, and metastasis. Saquinavir is a protease inhibitor of HIV has shown in-vitro effect against SARS and MERS coronaviruses (Feng et al., 2014). Multiple hydrogens and hydrophobic bonds for the Vapreotida and Saquinavir may be responsible for the tight binding and better affinity among the 23 drug molecules (Hossain, Khan, et al., 2016; Hossain, Oany, et al., 2016).

Molecular dynamics study of screened drug molecules were checked to evaluate their stability at atomic scale. The average root mean square deviation (RMSD) of Epirubicin, Vapreotida, Saquinavir are 1.756, 2.008, 1.76 Å which indicates complex stability throughout the entire simulation (Islam et al., 2019). Additionally, the average radius of gyration for these three complexes were found as 22.307, 22.211, 22.301 Å and did not deviate much. Besides, the mean of the surface area of the three complexes was 14005.961, 14419.772, 14514.697 Å2. It was cleared from Figure 3 that simulation analysis of top three therapeutic drug molecules consistent with overall stable nature and complex rigidity (Mahmud et al., 2019).

In this ground, it is clearly seen that Epirubicin, Vapreotida, and Saquinavir may inhibit COVID-19 by synergistic interactions among the 23 potential inhibitors against SARS-CoV-19 main protease and those results pave the way in drug discovery although it has to be further validated by in vitro and in vivo investigations.

Conclusion

Scientists have chosen drug repurposing techniques as alternative options for approved or investigational drugs outside their defined indication for ease of utilization in an emergency with minor associated risk and clue for prompt target drug development and application with minimum cost involvement. Bioinformatics approaches, such as molecular docking and dynamics study with probable drug candidates binding to key target proteins of pathogens for neutralization and inhibition are widely acceptable effective techniques to validate the drug repurposing approach. We have evaluated potential drug candidates against the main protease of the SARS-CoV-2 to understand their binding sites and structural stability. The nature of the binding mode of drug compound found as stable and molecular dynamics simulation study supports their conformational rigidity. As these outcomes solely came from combining computational algorithms, it requires further validation of these findings conducting further in-vitro wet-lab experimentation followed by in-vivo rapid trials for further use of the drugs individually or in combination to reduce mortality during Covid-19 pandemics globally. Moreover, based on these data, further, the experimental study can aid in getting suitable new drug molecules for combating COVID-19.