Identification of Papain-Like Protease inhibitors of SARS CoV-2 through HTVS, Molecular docking, MMGBSA and Molecular dynamics approach.

Coronaviruses (CoVs) are a large group of enveloped positive sense single-stranded RNA viruses that can cause disease to humans. These are zoonotic having potential to cause large-scale outbreaks of infections widely causing morbidity and mortality. Papain-Like Proteases (PLpro) is a cysteine protease, essential for viral replication and proliferation, as a highly conserved enzyme it cleaves peptide linkage between Nsp1, Nsp2, Nsp3, and Nsp4. As a valid therapeutic targets it stops viral reproduction and boosts host immune response thereby halting further spread of infection. In the purpose of identifying inhibitors targeting Papain-Like Proteases (PLpro) we initiated a high throughput virtual screening (HTVS) protocol using a SuperNatural Database. The XP docking results revealed that two compounds SN00334175 and SN00162745 exhibited docking score of -10.58 kcal/mol and -9.93 kcal/mol respectively. The Further PRIME MMGB-SA studies revealed Van der Waal energy and hydrophobic energy terms as major contributors for total binding free energy. The 100 ns molecular dynamics simulation of SN00334175/7JN2 and SN00162745/7JN2 revealed that these complexes were stabilized with ligand binding forming interactions with Gly266, Asn267, Tyr268, Tyr273, Thr301 and Asp302, Lys157, Leu162, Asp164, Arg166, Glu167, Pro248, Tyr264.

Introduction

Coronaviruses (CoVs) are a large group of enveloped positive sense single-stranded RNA viruses that can cause disease to humans (V'kovski et al., 2021). These are zoonotic having potential to cause large-scale outbreaks of infections widely causing morbidity and mortality (Kumar et al., 2021). The 2019 -2020 viral pneumonia epidemics in Wuhan are caused by a newly found Coronavirus disease (COVID-19) (Lakshmi Priyadarsini and Suresh, 2020). It is primarily transmitted through droplets produced by infected peoples when they are coughing, sneezing, or exhaling. The symptoms of recently emerged coronavirus COVID-19 manifests common cold, fever to serious respiratory disorders like severe acute respiratory syndrome (SARS) (Adithya et al., 2021). Currently, 3.75 million deaths with 174 million infected cases reported worldwide. India is facing daily rise in coronavirus cases with an active 29 million cases and 0.35 million deaths. In the present situation of viral outbreak which is spreading like a wildfire with a limited therapeutics as arsenal in hand, made us to rapidly identify new therapeutic agents with clinical use (Umakanthan et al., 2020). Since earlier days being as a basis and source of current medicines the natural products are widely acceptable chemical entities for therapeutic use. The pharmacological and therapeutic benefits with minimal toxicity of natural compounds play an important role in development of more effective drugs. In this study, we want to highlight yet untapped potential of natural compounds as antiviral agents from the plant source. The viral structure was built with major structural proteins like membrane (M) protein, nucleocapsid (N) proteins, spike (S) protein, and the envelop (E) protein, protease, hemagglutinin esterases, and helicases with other proteins which play key role in viral entry, survival and replication. The absence of proteins in human counterparts made them as attractive targets (Raj, 2021; Yadav et al., 2021). Papain-Like Proteases (PLpro) is a cysteine protease, essential for viral replication and proliferation and as a domain of NSP3 it is crucial component of replicase-transcriptase complex (RTC). The highly conserved enzyme cleaves peptide linkage between Nsp1 and Nsp2, Nsp2 and Nsp3, and Nsp3 and Nsp4 (Osipiuk et al., 2021). Proteases play a vital role in replication and maturation of viral proteins and PLpro are valid therapeutic targets which by inhibiting stops viral reproduction and boosts host immune response thereby halting further spread of infection (Lim et al., 2000; Harcourt et al., 2004; Lei et al., 2018; Mielech et al., 2014). As part of our continuous study using in silico and wet lab methodologies to find active compounds for various biological activities, (Rajagopal et al., 2017) we have designed and evaluated various compounds for their biological activities such as Anti-cancer, anti SARS CoV-2 etc (Kalirajan et al., 2012a, 2012b, 2018, 2019a, 2019b, 2020b, 2020c). In Schrödinger suite LLC, there are various modules like Glide, Qikprop, Prime, Desmond etc. The present study aims to identify inhibitors targeting Papain-Like Proteases (PLpro) through HTVS protocol using SuperNatural Database comprising of 325000 natural compounds. Sequential docking of HTVS, SP and XP modes were used to identify the binding modes. Further MMGB-SA post docking minimization was done and molecular dynamics simulations to identify natural compounds binding patterns in the PLpro active site (7JN2.pdb). The flow chart that follows explains the above-mentioned investigations (Figure 1 ).

Figure 1

Flow Chart illustrating step wise protocol for the identification of SARS CoV-2 inhibitors targeting Papain-Like Proteases (PLpro).

Materials and methods

High-throughput virtual screening and Molecular docking

The 3D x-ray crystal structure of Papain-Like Protease co-crystallized with naphthalen-1-yl- ethyl benzamide derivative (7JN2.pdb) was retrieved from protein data bank and was then further prepared using the protein preparation wizard of Schrödinger suite 2020-1(Sastry et al., 2013). The protein was prepared by eliminating crystal waters and bond orders adjusted with addition of hydrogen's (Sunitha Sukumaran et al., 2020). At pH 7.0, Prime was used to add missing side chains and loops, and then protonation and tautomeric states for acidic and basic residues were generated. The OPLS3e (Optimized Potentials for Liquid Simulations) molecular force field (Roos K et al., 2019) was used to minimize the protein, with the RMSD of crystallographic heavy atoms kept at 0.30. Å Grid box was generated (x = 52.16; y = 30.59; z = -0.62) at the centroid of active site keeping the Van der Waals scaling of 0.8 for the receptor with 0.15 as the partial charge cut-off. The docking protocol was validated using re-docking of Co-crystal ligand followed by calculating the RMSD difference between initial energy minimized pose and XP docked poses of Co-crystal ligand. The overlay of energy minimized initial pose and XP docked pose with RMSD difference was illustrated in Figure 2 . The 3D conformers of 3, 25, 000 natural compounds were obtained from SuperNatural Database (Dunkel et al., 2006). A virtual workflow was initiated with preparation of ligands using Ligprep and prefilter option to exclude redundant ligands. The ligands were sequentially docked into the catalytic pocket of Papain-Like Protease (7JN2.pdb) using HTVS, (High Throughput Virtual Screening), SP (standard precision) and XP (extra precision) modes keeping default parameters. The Best docking pose was selected based on glide g-score, glide energy values and hydrogen bond analysis.

Figure 2

Overlay of initial pose (green) and XP docked poses (gray) with RMSD of 1.8 Å.

MMGB-SA free energy calculation studies

The contributions of enthalpy and entropy related components towards binding of ligand-protein complex was calculated using Prime MMGB-SA approach (Shivakumar D et al., 2010) which integrates Generalized-Born/Surface Area (GB/SA) continuum solvent model (Naresh et al., 2020). The contributions from molecular mechanics energies, polar solvation and a nonpolar solvation terms were estimated (kcal/mol) using equation:

ΔG bind = Calculated binding free energy of complex

G complex = Binding free energy of minimized complex

G protein = Binding free energy of receptor

G ligand = Binding free energy of unbound ligand

Molecular Dynamics (MD) simulation (Kalirajan et al., 2017, Rajagopal et al., 2020a, Rajagopal et al., 2020d)

In order to study binding behavior of high ranked compounds at atomic level and to understand the molecular interaction analysis we performed molecular dynamics simulation using Desmond module of Schrödinger 2020-1, LLC, New York, NY (Bowers et al., 2006; Jorgensen et al., 1983). The complexs of SN00162745/7JN2 and SN00334175/7JN2 were solvated with TIP3P water model (Essmann et al., 1995) in an orthorhombic periodic boundary conditions having dimensions of 10 Å buffer region between protein atoms and box edges. The solvated system was neutralized by adding 0.15M NaCl counter ions. Then the system was minimized using default OPLS3e force field (Harder et al., 2016) parameters. The long range electrostatic interactions were calculated at a tolerance of 1e-09 using smooth particle mesh Ewald method. The short-range Van der Waals and Coulomb interactions were calculated at cut-off radius of 9.0 Å. A total of 100 ns MD simulation was performed with a time step of 2 fs under an isothermal-isobaric ensemble (NPT) at a temperature of 300 K and pressure of 1 bar. Nose-Hoover chain thermostat (Martyna et al., 1992) and Martyna-Tobias-Klein barostat (Martyna et al., 1994) methods were ensembled at 100 and 200ps respectively. A multiple time-step algorithms RESPA (REference System Propagator Algorithm) was used at 2, 2 and 6 fs for bonded, short-range non-bonded and long range electrostatic forces respectively. The data was collected at every 100ps and the obtained trajectories were analysed.

In silico predicted ADMET properties of the top hits

Prediction of physically and pharmacokinetically significant descriptors for the top hits (Compounds 1-28) were predicted employing various tools such as QikProp module of Schrodinger suite, ChemAxon tools and DataWarrior. Properties like Molecular weight, Total Polar Surface Area (TPSA), Hydrogen bond acceptor and donor count, log P, log D, log S, Molar volume, Dissociation constant (KD), No. of violations in lipinski's rule of five, Van der Waals Volume etc were given in Supplementary Table S1.

Results and discussion

Molecular docking and binding free energy calculation

A structural based virtual screening work flow from Schrodinger suite was performed against a library of 3, 25, 000 compounds from SuperNatural Database. The virtual screening work flow was carried out with initial prefilters which will screen out ligands based on Lipinski's rule and excluding ligands with reactive functional groups. The following sequential docking protocol was performed at 3 accuracy stages, HTVS, SP (standard precision) and XP (extra precision) docking keeping default constraints. The binding poses and hydrogen bond formation for the top scored hits after XP docking was visually analyzed and a final 28 compounds were identified. Post docking minimization was done employing MMGB-SA of PRIME module. The virtual screening protocol identified wide variation in scaffold topology like pyranyl hexanoic acid, pyran trihydroxy benzoates, 2-phenyl chromene ring linked with sugar moieties, tetrahydroxy hexanal linked with pyran, dicarbamimidamido cyclohexyl dihydrogen phosphates linked with sugars and benzyloxy benzoates linked with pyranoacetates etc. The glide scores and PRIME MMGB-SA energy values were given in Tables 1 and 2 .

Table 1

Glide XP docking values (kcal/mol) for the obtained hits in the active site of Papain-Like Protease (7JN2.pdb)

S. No	Compound Code	Glide gscore	Glide evdw	Glide ecoul	Glide energy	Glide emodel
1.	SN00216711	-9.742	-32.139	-15.09	-47.229	-48.899
2.	SN00007464	-9.805	-17.465	-24.569	-42.035	-56.702
3.	SN00334175	-10.581	-45.019	-16.152	-61.171	-83.343
4.	SN00216710	-9.506	-21.625	-26.927	-48.551	-49.238
5.	SN00293542	-9.473	-10.005	-25.296	-35.301	-54.168
6.	SN00168969	-9.397	-38.426	-19.934	-58.36	-88.319
7.	SN00213181	-9.387	-13.771	-28.123	-41.894	-59.98
8.	SN00213037	-9.33	-32.933	-23.816	-56.749	-76.776
9.	SN00249174	-9.052	-32.021	-35.305	-67.327	-86.228
10.	SN00162745	-9.939	-35.354	-16.911	-52.265	-60.996
11.	SN00041592	-8.906	-14.207	-24.167	-38.374	-55.815
12.	SN00330810	-8.83	-45.568	-9.492	-55.059	-83.974
13.	SN00403420	-8.783	-35.484	-23.999	-59.483	-81.909
14.	SN00379716	-8.598	-46.416	-14.061	-60.477	-89.631
15.	SN00165563	-8.333	-41.955	-20.633	-62.588	-83.089
16.	SN00338961	-8.206	-44.829	-18.619	-63.448	-64.351
17.	SN00350811	-7.872	-37.225	-18.056	-55.281	-72.941
18.	SN00296151	-7.613	-16.738	-24.797	-41.535	-52.462
19.	SN00215944	-7.457	-28.629	-21.661	-50.29	-73.531
20.	SN00143458	-7.455	-26.779	-14.116	-40.895	-60.705
21.	SN00162335	-7.227	-47.689	-12.542	-60.231	-79.355
22.	SN00382835	-7.122	-45.247	-10.131	-55.378	-79.424
23.	SN00332128	-6.972	-45.636	-11.765	-57.4	-60.013
24.	SN00040401	-6.936	-29.496	-21.513	-51.009	-69.388
25.	SN00347999	-6.836	-25.266	-22.44	-47.705	-60.952
26.	SN00038781	-6.769	-36.739	-14.896	-51.635	-74.733
27.	SN00175930	-6.679	-40.803	-13.906	-54.709	-78.135
28.	SN00334894	-6.554	-30.644	-19.429	-50.073	-34.005

Table 2

Prime MMGB-SA binding free energy values (kcal/mol) for the obtained hits in the active site of Papain-Like Protease (7JN2.pdb)

S. No	Compound Code	1MMGBSA Δ Bind	2MMGBSA ΔCoul	3MMGBSA ΔHbond	4MMGBSA ΔLipo	5MMGBSA ΔvdW
1.	SN00216711	-39.043	-6.822	-2.485	-16.307	-42.343
2.	SN00007464	-19.537	-12.244	-4.446	-9.9376	-18.745
3.	SN00334175	-58.344	13.603	0.218	-20.884	-57.851
4.	SN00216710	-39.499	-2.427	-1.416	-19.259	-38.890
5.	SN00293542	-16.242	6.6227	-3.175	-11.430	-26.193
6.	SN00168969	-46.517	9.8302	-1.121	-19.897	-45.915
7.	SN00213181	-29.152	0.8306	-1.851	-6.455	-23.522
8.	SN00213037	-37.953	11.684	-0.691	-7.708	-43.331
9.	SN00249174	-39.916	-17.717	-2.891	-3.213	-27.760
10.	SN00162745	-63.87	-25.085	-2.994	-16.865	-37.245
11.	SN00041592	-32.035	-11.768	-2.790	-14.923	-35.033
12.	SN00330810	-40.522	-8.150	-1.471	-16.005	-41.747
13.	SN00403420	-58.27	32.513	-1.892	-15.811	-43.606
14.	SN00379716	-63.729	-17.676	-1.647	-19.963	-48.175
15.	SN00165563	-54.368	13.583	-1.457	-17.197	-57.794
16.	SN00338961	-48.950	-9.735	-2.611	-17.182	-40.746
17.	SN00350811	-41.55	-3.480	-1.819	-15.571	-48.863
18.	SN00296151	-24.215	-9.008	-0.781	-14.584	-37.339
19.	SN00215944	-52.131	-43.926	-0.185	-5.7200	-19.320
20.	SN00143458	-36.163	21.586	-0.444	-16.304	-35.707
21.	SN00162335	-59.915	22.470	-0.378	-22.370	-49.717
22.	SN00382835	-54.554	-7.905	-1.924	-30.405	-56.870
23.	SN00332128	-53.936	13.051	-0.920	-24.397	-50.087
24.	SN00040401	-45.478	-11.242	-2.986	-14.937	-48.044
25.	SN00347999	-23.027	10.675	-0.021	-8.514	-40.120
26.	SN00038781	-40.170	56.419	1.947	-22.992	-53.651
27.	SN00175930	-48.972	4.066	0.7442	-20.836	-37.671
28.	SN00334894	-11.142	14.1967	-4.218	-7.876	-34.546

Binding free energy.

Coulombic energy.

Hydrogen bond energy.

Lipophilic energy.

Van der Waal energy.

The 3D interactions of top 10 hits were illustrated in Figure 3 . The XP docked poses of selected first 10 hits exhibited hydrogen bonding (5-10) and hydrophobic interactions mainly with catalytic residues Lys157, Asp164, Arg166, Glu167, Pro248, Tyr264, Gly266, Asn267, Tyr268, Gln269, Tyr273 and Thr301 (Table 3 ). The glide score (Table 1) ranges between -10.58 kcal/mol to -6.55 kcal/mol. The highest glide score was observed for the hits SN00334175 (-10.58 kcal/mol) and SN00162745 (-9.93 kcal/mol).

Figure 3

3D interaction diagram of top ten hits in the catalytic pocket of Papain-Like Protease (7JN2.pdb). (a) SN00216711; (b) SN00007464; (c) SN00334175; (d) SN00216710; (e) SN00293542; (f) SN00168969; (g) SN00213181; (h) SN00213037; (i) SN00249174; (j) SN00162745

Table 3

Number of hydrogen bonds and interacting residues for the obtained hits in the active site of Papain-Like Protease (7JN2.pdb)

S. No	Compound Code	No. of hydrogen bonds	Interacting residues
1.	SN00216711	5	Asp164, Arg166, Tyr268, Gln269
2.	SN00007464	6	Asp164, Arg166, Glu167, Tyr264
3.	SN00334175	5	Arg166, Tyr264, Gly266, Tyr268, Gln269
4.	SN00216710	7	Asp164, Arg166, Pro248, Gly266, Gln269
5.	SN00293542	8	Asp164, Arg166, Glu167, Gln269, Thr301
6.	SN00168969	5	Asp164, Arg166, Pro248, Tyr273
7.	SN00213181	10	Asp164, Arg166, Glu167, Gln269
8.	SN00213037	6	Asp164, Arg166, Pro248, Gln269, Tyr273
9.	SN00249174	6	Lys157, Arg166, Glu167, Tyr264, Tyr268, Gln269, Thr301
10.	SN00162745	7	Pro248, Tyr264, Gly266, Asn267, Tyr268, Tyr273, Thr301

From the obtained hits the compound SN00334175 (-10.58 Kcal/mol) (Figure 3c) exhibited five hydrogen bonds with the active site residues of Papain-Like Protease, carbonyl oxygen formed hydrogen bond with backbone NH of Gln269 (-C=O…HN-; 2.02Å). The two carbonyl oxygens of pyran acetates formed two hydrogen bonds with side chain NH and OH of Arg166 (-C=O…HN-; 2.15Å) and Tyr268 (-C=O…HO-; 2.68Å) respectively. The hydroxyl groups on pyran and benzene rings formed two hydrogen bonds with Gly266 and Tyr264 (-OH…OH) with 1.86 and 2.00Å respectively. Another hit compound SN00162745 (-9.39 Kcal/mol) (Figure 3 j) formed seven hydrogen bonds, the p-hydroxyl group on 2-phenyl ring of chromene ring accepted and donated electrons forming two hydrogen bonds with Tyr273 (1.91Å) and Thr301 (2.07Å). The hydroxyl groups on pyranose and furanose rings formed five hydrogen bonds with Asn267, Gly266 and Pro248. Two hydrophobic interactions by Tyr268 and Tyr268 stabilized the chromene ring. Post docking minimization for the top ranked poses of the selected hits exhibited binding free energy (ΔBind) ranging from -11.14 to -63.87 Kcal/mol. As mentioned in Table 2 the major Van der Waal energy terms (ΔVdW) -15.79 to -57.79 Kcal/mol favor total binding energy with moderately favored hydrophobic energy terms (ΔLipo) -3.21 to -30.40 Kcal/mol.

In silico predicted ADMET properties

The in-silico ADMET properties for the top 28 hits were predicted using various software tools. Molecular weight of the compounds were observed in the range 358 to 990.8 g.mol−1. Total polar surface area (TPSA), which is surface sum of polar atoms was observed between 157 and 506 Å2. The lipophilicity Log P values of the compounds are in the range of -11.2 to 5.2 and Log D values ranges between -13.8 to 6.1. The molar volume and Van der waal volume were represented in the Supplementary Table S1. The obtained ADMET properties are within the recommended values with fewer exceptions.

Molecular dynamics simulation

A 100ns molecular dynamics simulation for the docked pose of SN00334175/7JN2 complex revealed that the RMSD of protein Cα atoms (Figure 4 a) were stabilized after ligand binding showing minute fluctuations from 2.1 to 2.7 Å. The ligand RMSD (Figure 4 a) exhibited major fluctuations with 5.9-10.2 Å from 35 to 70 ns which then stabilized showing 6.5 to 8.0Å till 100 ns. The RMSF (Figure 4b) of protein looks stable except amino acids from 180 to 240 showed higher fluctuations till 4.8 Å which are present in loops. From the ligand interaction fraction (Figure 4c) the compound was positioned in the active pocket by forming hydrogen bonding, hydrophobic and water bridged interactions with Lys157, Leu162, Asp164, Arg166, Glu167, Pro248, Tyr264, Tyr268, Gln269, Tyr273 and Thr301. From the 2D ligand interaction diagram (Figure 4d) of the 100ns simulation, the carbonyl oxygen of benzoxy benzoate group exhibited two weak hydrogen bonds with Gln269 and Tyr264 at 20% and 21% of total simulation trajectory. Asp164 formed extensive water bridges interactions with triacetyl pyran ring. The phenyl hydroxyl group formed another water bridged interaction with Thr301. For the other complex a 100 ns molecular dynamics simulation for the docked pose of SN00162745/7JN2 complex revealed that RMSD of the protein Cα atoms (Figure 5 a) were fluctuated from 1.4 to 3.2 Å. The ligand RMSD (Figure 5a) exhibited fluctuations with 1.2-2.8 Å. The RMSF (Figure 5b) of protein looks stable after 50 ns except amino acids from 180, and 230 showed higher fluctuations till 4.0 Å which are present in loops. From the ligand interaction fraction (Figure 5c) the compound was positioned in the active pocket by forming hydrogen bonding, hydrophobic and water bridged interactions with Pro248, Ala249, Gly266, Tyr273 and Asp302. From the 2D ligand interaction diagram of the 100 ns simulation (Figure 5d), the p-hydroxyl group on the 2-phenyl ring formed strong hydrogen bonds Tyr273 and Asp302 with 97% and 94% of simulation trajectory. Whereas, the m-hydroxyl exhibited strong hydrogen bond interaction with 98% of simulation trajectory. The hydroxyl group on the pyranose ring formed water mediated hydrogen bond with Pro248. Two hydrophobic interactions were observed with Tyr264 and Tyr268. At last, the chromone hydroxy group also involved 82% with Gly266 amino acid residue in hydrogen bond formation.

Figure 4

(a) RMSD graph for the 100ns simulation trajectory of SN00334175/7JN2 complex (b) RMSF graph for the 100ns simulation trajectory of SN00334175/7JN2 complex (c) Ligand interaction fraction for the 100ns simulation trajectory of SN00334175/7JN2 complex (d) 2D interaction diagram for the 100ns simulation trajectory of SN00334175/7JN2 complex

Figure 5

(a) RMSD graph for the 100ns simulation trajectory of SN00162745/7JN2 complex (b) RMSF graph for the 100ns simulation trajectory of SN00162745/7JN2 complex (c) Ligand interaction fraction for the 100ns simulation trajectory of SN00162745/7JN2 complex (d) 2D interaction diagram for the 100ns simulation trajectory of SN00162745/7JN2 complex.

Conclusion

Currently, 3.75 million deaths with 174 million infected cases reported worldwide. India is facing daily rise in coronavirus cases with an active 29 million cases and 0.35 million deaths. The present study aims to identify inhibitors targeting Papain-Like Proteases (PLpro) which plays vital role in viral replication and proliferation. A high throughput virtual screening (HTVS) protocol was initiated with sequential docking protocol performed at 3 accuracy levels, HTVS, SP (standard precision) and XP (extra precision) modes of docking using a SuperNatural Database comprising of 325000 natural compounds. Further post docking minimization was performed identifying top 28 hit compounds. The XP docking results revealed that two compounds SN00334175 and SN00162745 exhibited docking score of -10.58 kcal/mol and -9.93 kcal/mol. The PRIME MMGB-SA studies indicated that total binding energy was majorly contributed from Van der Waal energy (ΔVdW) -15.79 to -57.79 kcal/mol and hydrophobic energy terms (ΔLipo) -3.21 to -30.40 kcal/mol. The 100 ns molecular dynamics simulations for the complexes SN00334175 and SN00162745 revealed that they were stabilized in the catalytic pocket forming hydrogen bonding, hydrophobic and water bridged interactions. Especially from molecular dynamics simulation the compound SN00162745 comprising 2-phenyl-chromene-4-one ring exhibited strong hydrogen bonds with Gly266, Asp302 and Tyr273 with 82%, 99% and 94% of total 100ns simulation trajectory. This study highlighted the importance of a glycosylated flavonoids as a starting point for designing Papain-Like Proteases (PLpro) inhibitors where chemical modifications with improved pharmacokinetic properties could be investigated clinically contributing development of potential anti SARS-CoV 2 agents.

Declaration of Competing Interest

The Authors Stated That there is no conflict of Interest.