Molecular docking and dynamic simulations for antiviral compounds against SARS-CoV-2: A computational study.

The aim of this study was to develop an appropriate anti-viral drug against the SARS-CoV-2 virus. An immediately qualifying strategy would be to use existing powerful drugs from various virus treatments. The strategy in virtual screening of antiviral databases for possible therapeutic effect would be to identify promising drug molecules, as there is currently no vaccine or treatment approved against COVID-19. Targeting the main protease (pdb id: 6LU7) is gaining importance in anti-CoV drug design. In this conceptual context, an attempt has been made to suggest an in silico computational relationship between US-FDA approved drugs, plant-derived natural drugs, and Coronavirus main protease (6LU7) protein. The evaluation of results was made based on Glide (Schrödinger) dock score. Out of 62 screened compounds, the best docking scores with the targets were found for compounds: lopinavir, amodiaquine, and theaflavin digallate (TFDG). Molecular dynamic (MD) simulation study was also performed for 20 ns to confirm the stability behaviour of the main protease and inhibitor complexes. The MD simulation study validated the stability of three compounds in the protein binding pocket as potent binders.

Introduction

Coronavirus (CoV) is a genus of the Coronaviridae family named for the crown-like spikes found on their surface. They are a huge family of viruses containing a genome composed of a long RNA strands which is the largest of all RNA viruses, and this genome acts like a messenger RNA when it infects a cell, and directs the synthesis of two long polyproteins that include the machinery that the virus needs to replicate new viruses. A novel betacoronavirus from the subgenus Sarbecovirus has been isolated from human airway epithelial cells [[1], [2], [3]]. Protease inhibitors could likely block a key enzyme that helps viruses replicate and prevents SARS, which is also a coronavirus [4]. There is some information about doctors administering HIV drugs for treating COVID-19 patients, but no solidly positive results. Anti-malarials have also been tested with little evidence for effectiveness. If the virus attacks immunocompromised patients, they will be at risk, and neither HIV retrovirals nor plant-based products are therapeutic [5]. Naturally occurring phyto-chemically based compounds have been shown to exhibit several anti-viral effects including other pharmacological properties [6]. The genome of novel coronavirus (SARS-COV-2) encodes many important proteins for its replication in the host genome viz. the nucleocapsid protein, Spike (S) protein, Envelope (E) protein, Membrane (M) protein and coronavirus main protease, which play crucial roles in gene expression, and cleave polyproteins into replication-related proteins [7,8]. With this in mind, the present in silico study was designed to evaluate the effects of FDA approved anti-viral drugs and plant-based antiviral agent on the COVID-19 main protease viral protein of SARS-COV-2.

Methodology

Protein preparation for docking

The X-ray diffraction-based crystal structure of COVID-19 main protease in complex with an inhibitor N3 with a resolution of 2.16 Å that contains neither carbohydrate polymers nor chain breaks was selected for the study [9]. The complexes bound to the protein receptor molecule were removed. The protein preparation wizard from the Schrödinger module was used to prepare the structure of the main protease on adding the hydrogen atoms, removing the waters beyond 5 Å of the binding site, and the active site grid was generated using the Receptor grid generation application in Glide module. Glide uses a filter search to locate the ligand in the active-site region of the receptor. The shape and properties of the receptor are represented on a grid that provides a more accurate scoring of the ligand poses. The grid of 20 Å was generated over the co-crystallized ligand inhibitor molecule. The docked complexes were superimposed to the original crystal structure to calculate the root mean square deviation (RMSD) using pyMOL.

Ligand preparation

Twenty-four natural plant-based compounds with antiviral property, 22 US FDA approved antiviral drugs, and 16 anti-malarial drugs were identified from the PubMed literature as test ligand molecules against Main protease receptor. LigPrep (Schrödinger) is used to test compounds on assigning chiralities and are converted to 3D structures (Fig. 1 a). Ionization and tautomeric states were generated using the OPLS_2005 force field. For each ligand, 32 stereoisomers were generated.

Fig. 1

a). Main protease (6LU7) and b). Calculation of RMSD and visualization of the Main protease + amodiquine complex with crystal structure 6LU7 using PYMOL.

Virtual screening

Initial screening was done with plant-based and antiviral drugs using a virtual screening workflow with default parameters using the Glide program of Schrödinger. The HTVS mode eliminates most of the stereoisomers and only a few of the isomers that passed after screening were subjected to pass through SP and XP docking modes. The ligands which showed a better affinity towards the main protease would be qualified. Finally, the interactions of selected ligand and protein docked complexes were analyzed by pose viewer.

Molecular dynamic simulations

The protein-ligand complex structure of SARS-CoV-2 main protease and candidate molecules were prepared for MD simulation using VMD. GROMACS-2019 version was used to carry out 20 ns simulations using the OPLS force field. The TIP3P water model was selected for solvating complexes followed by addition of ions to neutralize. Periodic boundary conditions were used. Energy minimization was done with a tolerance of 1000 kJ/mol/nm. Equilibration of the system was done using NVT and NPT ensemble for 100 ps. The trajectories were set to be generated every 2 fs and save every 2ps. The protein-ligand complexes results were then analyzed.

Results and discussion

Virtual screening and molecular docking

For each ligand, 32 conformations have been generated by the LigPrep module and they were further evaluated for Virtual screening using the Glide module. The HTVS mode identified active compounds for screening, and ligands with different conformations were screened and false positives eliminated. The LigPrep program produced results for 1129 ligand conformations of FDA approved antiviral drugs, 459 of anti-malarial drugs, and 110 of plant-based drugs. The stereo conformer per each ligand which showed the lowest energy with the receptor has been selected for further docking study under SP (Standard Precision) mode of Glide. The Protein-ligand interaction of the stable docked Lopinavir, Amodiaquine, and Theaflavin digallate with main protease complex was visualized with ligand interaction. Lower RMSD values of PyMOL superimposed docked complexes with respect to initial structure before docking indicates that our complexes are stable (Fig. 1b). Lopinavir showed interactions of hydrogen bond interactions with Glu 166, GLN 189, GLY 143 and pi -pi stacking with His 41 (Fig. 2 a and b) and was successfully docked with a best docking score of −9.918, glide energy −8.023, and glide e-model −76.898 (Table 1 ). Theaflavin digallate showed interactions of hydrogen bond interactions with THR 26, ASB 142, HIS 41, HIS 163, GLU166, GLN 189 (Fig. 3 a and b) docked with a best docking score of --10.574, glide G-SCORE -10.722 and glide e-model −135.584 (Table 2 ). The ligand amodiaquine shares four hydrogen bond interactions with Leu 141, ARG 187 and salt bridge interaction with Glu 166, whereas pi-pi stacking with His 41 (Fig. 4 a and b) of the active site residues of main protease and the docked complex orientation exhibited a docking score of −7.421, glide g-score of −8.023, and glide e-model −76.898 (Table 3 ). The LIGPLOT 2D interaction showed the interacting residue information. The results produced from docking showed the antimalarial drug and Theaflavin digallate. Hence, from in silico approach, the study suggests that the US FDA approved drugs Lopinavir, Amodiaquine, and Theaflavin digallate may be tested in vivo as potent drugs against SARS-CoV-2 (see Fig. 5).

Fig. 2

a). Docked pose of lopanivir (antiviral drug) molecule with main protease (6LU7) and b). Ligand interaction of lopanivir with 6LU7.

Table 1

Docking scores of US-FDA approved antiviral drugs using GLIDE module.

S.no	Entry Name/ID	Docking score	Glide g-score	Glide e-model
1	Lopinavir	−9.918	−9.918	−101.59
2	Darunavir	−8.843	−8.972	−84.355
3	Amprenavir	−8.655	−8.784	−89.621
4	Rupintrivir	−8.342	−8.342	−94.738
5	Sofosbuvir	−8.324	−8.324	−83.378
6	Adefovirdipivoxil	−8.252	−8.277	−98.829
7	Famciclovir	−7.546	−7.546	−52.595
8	Tecovirimat	−7.546	−7.546	−52.595
9	Darunavir	−7.505	−8.472	−83.862
10	Zidovudine	−7.396	−7.396	−60.709
11	Dolutegravir	−7.279	−7.727	−72.215
12	Entecavir	−7.15	−7.292	−66.857
13	Bictegravir	−7.088	−7.463	−68.443
14	Oseltamivir	−7.037	−7.049	−57.111
15	Emtricitabine	−6.941	−6.941	−56.097
16	Zalcitabine	−6.712	−6.712	−51.947
17	Didanosine	−6.614	−6.933	−47.279
18	Baloxavirmarboxil	−6.304	−6.313	−60.119
19	Emtricitabine	−6.21	−6.21	−46.753
20	Simeprevir	−5.824	−5.824	−66.915
21	Elbasvir	−4.966	−4.966	−51.332
22	Ritiometan	−3.671	−3.671	−49.254

Fig. 3

a). Docked pose of plant based phenol, Theaflavin digallate with main protease (6LU7) and b). Ligand interaction of plant based phenol, Theaflavin digallate with 6LU7.

Table 2

Docking scores of natural plant based molecules using GLIDE.

S. no.	Entry Name	docking score	glide gscore	glide emodel
1	Theaflavin digallate	−10.574	−10.722	−135.584
2	Biorobin	−9.058	−9.087	−97.726
3	Hesperidin	−7.848	−7.848	−87.457
4	Rosmarinic acid	−6.971	−6.971	−73.914
5	Berchemol	−6.793	−6.793	−64.803
6	Baicalin	−6.563	−6.57	−65.605
7	lycorine	−6.522	−6.555	−54.049
8	Chrysin	−6.461	−6.504	−47.835
9	Berberine	−6.429	−6.429	−53.057
10	(−)-Epigallocatechin gallate	−6.142	−6.221	−72.524
11	Hesperetin 7-O-neohesperidoside	−6.124	−6.124	−76.976
12	Clivimine	−6.109	−6.141	−66.176
13	Kouitchenside I	−5.845	−5.858	−70.472
14	Cosmosiin	−5.717	−5.717	−65.508
15	Piceatannol	−5.614	−5.614	−44.629
16	Gnidicin	−5.579	−5.579	−57.887
17	Rosmarinic Acid 3′-O-Beta-D-Glucoside	−5.246	−5.246	−65.572
18	Andrograpanin	−4.941	−4.941	−41.129
19	Andrographiside	−4.647	−4.647	−54.294
20	lycorine	−4.63	−6.366	−50.165
21	Cerevisterol	−4.468	−4.468	−40.912
22	Stigmast-5-en-3-ol	−4.395	−4.395	−37.882
23	Betulonal	−4.245	−4.245	−49.295
24	Andrograpanin acetate	−3.783	−3.783	−42.001

Fig. 4

a). Docked pose of amodiquine (antimalarial drug) molecule with main protease (6LU7) and b). Ligand interaction of amodiquine with 6LU7.

Table 3

Docking scores of anti malarial molecules with main protease (6LU7).

S.no	Entry Name	docking score	glide gscore	glide emodel
1	Amodiaquine	−7.429	−8.023	−76.898
2	Mefloquine	−6.873	−6.876	−53.961
3	Quinine	−6.508	−6.522	−50.551
4	Primaquine	−6.361	−6.361	−59.673
5	Halofantrine	−6.352	−6.354	−68.585
6	Lumefantrine	−6.202	−6.226	−62.308
7	Chloroquine	−6.075	−6.111	−56.559
8	Piperaquine	−5.748	−6.009	−69.174
9	Sulfadoxine	−5.516	−5.668	−52.974
10	Atovaquone	−5.493	−5.5	−47.59
11	Artenimol	−5.178	−5.178	−38.26
12	Atovaquone	−5.121	−5.128	−48.149
13	Artesunate	−4.862	−4.862	−45.28
14	Proguanil	−4.842	−5.252	−38.978
15	Artemether	−4.764	−4.764	−30.727
16	Piperaquine	−4.716	−6.879	−82.308

Fig. 5

Plot of root mean square deviation (RMSD)a, root mean square fluctuation (RMSF)b values, radius of gyration (Rg)c during 20ns MD simulation of SARS-CoV-2 3CL protease in complex with amodiquine.

Molecular dynamic simulation

The molecular dynamic simulations were examined on the basis of Root mean square deviation (RMSD), Root mean square fluctuation (RMSF), and Radius of gyration values as a function of time. The structure variation was calculated by RMSD values of protein-ligand complexes from 0 to 20 ns. The RMSD values steadily increased from 0 to 5 ns, and reached stable state throughout the simulation. The average RMSD values of amodiquine, lopinavir, and theaflavin digallate were found to be 0.25, 0.23, and 0.22 nm, respectively. RMSF provides the fluctuation of each atom in the overall simulation (Fig. 5, Fig. 6, Fig. 7 ). RMSF was calculated for protease with 306 amino acids and three potential drug candidates; the values confirmed that the binding site residues showed less fluctuation. The average RMSF values were 0.15, 0.17, and 0.2 nm for amodiquine, lopinavir, and theaflavin digallate, respectively (Fig. 5, Fig. 6, Fig. 7b). The radius of gyration (Rg) of the protein and ligand complexes was found to be between 2.13 and 2.18 nm initially. The Rg values of the main protease with amaodiquine was stabilized after 5 ns and in case of lopinavir and theaflavin digallate, the values were decreased initially and stabilized from 5 to 20 ns which was an indicator for a stable binding pose (Fig. 5, Fig. 6, Fig. 7 The MD simulation analysis of selective docked complexes confirmed the stability of protein ligand complexes. Dayer et al. (2017) conducted 50ns MD simulations on HIV-1 protease and inhibitors complex and they found that lopinavir was possible drug candidate with a significant binding energy value [10]. Nukoolkarn et al. (2008) performed 2ns MD simulations with main protease and inhibitor complexes (ritronavir and lopinavir) and reported that the ligand-binding site of the SARS-CoV 3CLpro (main protease) was around active sites of H41 and C145 [11]. Similarly, in the present study, we found that the active site of H41 was sharing hydrogen bonds with all three proposed inhibitors.

Fig. 6

Plot of root mean square deviation (RMSD)a, root mean square fluctuation (RMSF)b values, radius of gyration (Rg)c during 20ns MD simulation of SARS-CoV-2 3CL protease in complex with lopanivir.

Fig. 7

Plot of root mean square deviation (RMSD)a, root mean square fluctuation (RMSF)b values, radius of gyration (Rg)c during 20ns MD simulation of SARS-CoV-2 3CL protease in complex with theaflavin digallate.

Conclusion

The tentative agreement of applying drug candidates such as Lopinavir, an antiviral HIV-drug, amodiquine, an antimalarial drug, and Theaflavin digallate, a plant-based phenol derivative, were selected as SARS‐CoV‐2 therapeutics due to their good binding affinity for the active site of main protease. Identification of natural and existing approved drugs is a noteworthy step for early drug discovery against COVID‐19. Based on the observations of docking score, we believed that the phenol derivative and anti-HIV drug, could aid in COVID‐19 drug discovery. The RMSD of three complexes fall between 0.15 and 0.25 nm and inferred that the compounds lopinavir, amodiquine, and theaflavin digallate had undergone good conformational changes while binding, and maintained close affinity with the binding site of the main protease. However, in vitro and in vivo evaluation study is required to repurpose these three drugs against 2019-nCoV.

Ethical statement

The study does not involves the animals.

Declaration of competing interest

The authors declare no competing interests.