"Identification of Nafamostat and VR23 as COVID-19 drug candidates by targeting 3CLpro and PLpro."

The sudden increase in the COVID-19 epidemic affected by novel coronavirus 2019 has jeopardized public health worldwide. Hence the necessities of a drug or therapeutic agent that heal SARS-CoV-2 infections are essential requirements. The viral genome encodes a large Polyprotein, further processed by the main protease/ 3C-like protease (3CLpro) and papain-like proteases (PLpro) into 16 nonstructural proteins to form a viral replication complex. These essential functions of 3CLpro and PLpro in virus duplication make these proteases a promising target for discovering potential therapeutic candidates and possible treatment for SARS-CoV-2 infection. This study aimed to screen a unique set of protease inhibitors library against 3CLpro and PLpro of the SARS-CoV-2. A molecular docking study was performed using PyRx to reveal the binding affinity of the selected ligands and molecular dynamic simulations were executed to assess the three-dimensional stability of protein-ligand complexes. The pharmacodynamics parameters of the inhibitors were predicted using admetSAR. The top two ligands (Nafamostat and VR23) based on docking scores were selected for further studies. Selected ligands showed excellent pharmacokinetic properties with proper absorption, bioavailability and minimal toxicity. Due to the emerging and efficiency of remdesivir and dexamethasone in healing COVID-19 patients, ADMET properties of the selected ligands were thus compared with it. MD Simulation studies up to 100 ns revealed the ligands' stability at the target proteins' binding site residues. Therefore, Nafamostat and VR23 may provide potential treatment options against SARS-CoV-2 infections by potentially inhibiting virus duplication though more research is warranted.

Introduction

In December 2019, a new coronavirus caused an outbreak of the pulmonary disease in Wuhan, the capital of Hubei province in China, and has since spread globally [[67], [68], [79]]. The virus has been named SARS-CoV-2 [20], with 96% genome identical to a bat coronavirus and shares 79.6% sequence identity to SARS-CoV [43,67,79]. This pandemic spread worldwide with more than 104.9 million infections and more than 2.27 million deaths till 3rd February 2021 [https://www.worldometers.info/coronavirus/].

The genome size of coronaviruses is ~30,000 nucleotides in length with a 5′-cap structure and a 3′-poly (A) tail and consist of at least six open reading frames (ORFs) [27,13]. The first ORF (ORF 1a/b) is about two-thirds of the genome length, precisely translates two polyproteins, pp1a and pp1ab and are processed by the main protease, also known as the 3C-like protease (3CLpro) and by one or two papain-like proteases (PLpro), into 16 nonstructural proteins (NSPs) and develop into mature proteins which assist in the viral replication [7,26,24]. PLpro facilitates cleavage at the first three polyproteins sites, whereas CLpro facilitates the cleavage at 11 sites [14,29]. The 3CLpro carries out cleavage at the polyproteins' C-terminal while, N-terminal of the polyproteins are cleaved by PLpro [40]. These NSPs engaged in subgenomic RNAs construction that encodes four main structural proteins, namely envelope (E), membrane (M), spike (S), and nucleocapsid (N) proteins and other accessory proteins [56,57].

The ~306 aa long 3CLpro, a key enzyme for coronavirus replication, is also encoded by the polypeptide and responsible for processing the polypeptide into functional proteins [67,77]. The 3CLpro, also known as Nsp5, is the first to get automatically cleaved from polyproteins to produce mature enzymes, and then it further assists in the cleavage of downstream Nsps at 11 sites to release Nsp4–Nsp16 [71]. 3CLpro directly mediates the maturation of Nsps, which is essential in the life cycle of the virus. Studies have shown that 3CLpro of different coronaviruses and all the 3CLpro cleavage sites on polyprotein 1ab are highly conserved in sequences and 3D structures [72,69]. Since there are no known human proteases with a similar cleaving specificity, 3CLpro inhibitors are less likely to be toxic and, together with its functional importance, have made 3CLpro an attractive target for the design of anti-SARS-CoV-2 drugs [2,69,77].

PLpro is responsible for the cleavages of N-terminus of the replicase polyprotein at three conserved cleavage sites to release Nsp1, Nsp2 and Nsp3, essential for correcting virus replication [21,7]. PLpro was also confirmed to be significant to antagonize the host's innate immunity [12,75,37]. As an essential enzyme in the host's coronavirus replication and infection, PLpro has been a popular and valuable target for identifying potential drugs against novel coronavirus 2019 [42].

Molecular docking has become a promising tool for drug discovery and development. By utilizing this tool, the binding interaction of drug-like molecules inside the target protein's binding pockets can be investigated with other factors like identification of hit molecules, optimization of lead compound and virtual screening. Since COVID-19 is a significant outbreak in 215 nations and there are no effective drugs for COVID-19, it is challenging to cure SARS-CoV-2 infection and control the associated pandemic. Therefore, novel drug design and discovery techniques can be utilized to discover potential therapeutic candidates against SARS-CoV-2.

Materials and methods

Protein structure retrieval

The crystal structures for 3CLpro (PDB ID: 6Y2E, Resolution: 1.75 Å) [77] and PLpro (PDB ID:6W9C, Resolution: 2.70 Å) of SARS-CoV-2 were downloaded from the PDB database (https://www.rcsb.org) in PDB format and were utilized for the in silico studies.

Energy minimization and structure validation

YASARA Energy Minimization Server carried out energy minimization of the crystal structures of 3CLpro and PLpro of the SARS-CoV-2 to obtain an energy-minimized and highly stable protein structure for efficient docking [33]. PROCHECK further validated both the energy minimized crystal structures of SARS-CoV-2 [36].

Binding site prediction

Binding site residues were anticipated through a literature survey for 3CLpro [[32], [58], [61], [70], [77]] and PLpro [4,66] (Table S1).

Ligand selection and ligand file preparation

Protease inhibitors, a class of antiviral drugs, are widely used to treat HIV/AIDS and hepatitis C. Inhibition of serine proteases mitigated SARS-CoV pathogenesis in vivo was also reported. Protease inhibitors prevent viral replication by selectively binding to viral proteases and blocking proteolytic cleavage of protein precursors necessary for the production of infectious viral particles [[17], [44], [52], [63]]. The protease inhibitor library, with a unique collection of 227 small molecule inhibitors used for chemical genomics, high-throughput screening (HTS) and high content screening (HCS), was downloaded from Selleckchem.com (https://www.selleckchem.com/) and saved in sdf file. According to our docking protocol, the MDL MOL format was converted to a pdbqt file by Open Babel [49].

Molecular docking and interaction analysis

Validation of docking protocol is essential in molecular docking to ensure that ligands bind in the correct conformation within the protein's binding site pocket, which is done by validating the size and center of the coordinates of the grid box across the binding pocket [39]. PyRx, a virtual screening software for computational drug discovery, was used to screen the selected ligands (set of protease inhibitors) files against SARS-CoV-2 3CLpro and PLpro [15]. PyRx uses AutoDock 4 and AutoDockVina as docking software implying the Lamarckian Genetic Algorithm and Empirical Free Energy Scoring Function. PyRx was carried out using selected inhibitors' on the energy minimized structure of SARS-CoV-2 3CLpro and PLpro. The macromolecular structure of two proteases (3CLpro & PLpro) and the ligands were prepared and docking was accomplished into the binding site residues inside a grid box with X, Y and Z axis and dimensions adjusted to 16.17 Å × 17.25 Å × 20.14 Å and 49.43 Å × 49.41 Å × 21.06 Å for SARS-CoV-2 3CLpro and PLpro protein respectively. As in PLpro, the catalytic triad Cys111-His272-Asp286 in each of three subunits A, B and C respectively are directed to the center of the protein [66]; the binding pocket for all the three subunits was taken under consideration and was within the validated size and dimensions of the grid box in PyRx docking. The docking protocol was then run at exhaustiveness of 8 and set to only output the lowest energy pose. The interactions between our targeted proteins and the ligands were studied using Ligplot [64], and the pictures were processed and made using Pymol Molecular Visualization Software [38].

Pharmacokinetics studies

The drug-likeness prediction of the selected ligands was carried out by Lipinski filter (http://www.scfbio-iitd.res.in/software/drugdesign/lipinski.jsp), according to which an orally active drug should acquiesce to a minimum of four of the five laid down criteria for drug-likeness, namely: molecular mass, H-bond acceptor, H-bond donor, molar refractive index and LogP [41]. Also, molinspiration (https://molinspiration.com/cgi-bin/properties) was used to calculate the selected inhibitors' molecular properties and bioactivity. The absorption percentage (% AB) was calculated using the formula [[27], [78]]:AB% = 109 – (0.345 x TPSA).

The pharmacodynamics parameters: Absorption, Distribution, Metabolism, Excretion and the Toxicity (ADMET) of the inhibitors were predicted using admetSAR (http://lmmd.ecust.edu.cn:8000/).

Molecular dynamics (MD) simulation

All-atoms MD simulations were performed on the atomic coordinates of the best-docked complex of SARS-CoV-2 target proteins, i.e., 3CLpro and PLpro using Gromacs v5.1.4 with force field CHARMM27 and water model TIP3P [1] and the ligands parameters were defined from Zoete et al. [80]. The simulation box was defined with buffer distance (10 Å) from the centrally placed protein-ligand complex. The prepared system was solvated with water molecules and neutralized with the addition of 0.15 M counter ions (Na+ and Cl−) [31]. The energy minimization process involves 50,000 steps for each steepest descent, followed by conjugant gradients. PBC condition was defined for x, y and z directions [16] and simulations were performed at physiological temperature; 300 K. SHAKE algorithm was applied to constrain all bonding involved, hydrogen and long-range electrostatic forces treated with PME (Particle mesh Ewald). The system was equilibrated in two steps, NVT and NPT, at 300 K for a period of 500 ps. During the simulation, the Berendsen thermostat [9] and Parrinello-Rahman pressure [51] were used to maintain pressure and temperature. LINC algorithm was used to constrain the bonds and angles [23]. The van der Waals interactions are taken care of by LJ potential with a cutoff of 0.10 nm. Using the NPT ensemble, production runs were performed for the period of 100, with time integration. The energy, velocity and trajectory were updated at the time interval of 10 ps. All production runs were done on CUDA enabled Tesla GPU machine (DELL T640 with V100 GPU) and OS Centos 7 [59,53] and the Gromacs utilities were used for the analyses of obtained MD trajectories.

Results and discussions

Energy minimization and model validation

YASARA [33] was used for energy minimization of the crystal structures for SARS-CoV-2 3CLpro (PDB ID: 6Y2E) and PLpro (PDB ID:6W9C) for low energy and high stability of the protein and was validated using PROCHECK [36]. Ramachandran plot analysis showed that 99.6% residues are in favoured, allowed and generously allowed regions and only 0.4% in the disallowed region for the energy minimized 3CLpro protein. In comparison, 100% residues are in favoured, additionally allowed and generously allowed regions for the PLpro protein (Fig. S1).

The top challenging job in computational chemistry is predicting ligands' actual binding conformation into protein binding sites [6]. After successful docking between the selected ligands (set of protease inhibitors) and the binding site of SARS-CoV-2 3CLpro and PLpro, the ligands' docking scores were noted and the best-ranked poses with the lowest docking score were chosen for interaction studies in detail. The mode of interaction of inhibitor compounds with amino acid residues in protein binding sites results in the compounds' binding affinity and potency [76].

The docking study carried out with the set of protease inhibitors as ligands against 3CLpro and PLpro of SARS-CoV-2 gave a formative revelation of molecular interplay. We find that all ligands interact with 3CLpro and PLpro more or less efficiently. Nafamostat and VR23 [7‑chloro-4-(4-((2,4-dinitrophenyl)sulfonyl) piperazin-1-yl) quinoline] showed the lowest binding energy and high efficiency against both the target protein, i.e., 3CLpro and PLpro of SARS-CoV-2. Nafamostat showed binding energy of −9.0 kcal/mol and −9.2 kcal/mol against 3CLpro and PLpro, respectively, while VR23 possessed binding energy of −9.1 kcal/mol against both 3CLpro and PLpro, respectively (Table 1 ). A recent study suggests that nafamostat mesylate can effectively kill SARS-CoV-2 in vitro by targeting spike protein [74]. Earlier reports also suggested that Nafamostat effectively inhibits MERS-CoV besides its extensive use in Japan to treat patients with a dissemination intravascular coagulation (DIC) and acute pancreatitis [73]. COVID-19 manifested DIC characteristics with enhanced fibrinolysis, and nafamostat was anticipated as a promising solution [5]. Nafamostat mesylate was considered a drug candidate for treating DIC induced by the Ebola virus infection [48]. Nafamostat has been used clinically to treat pancreatitis in Japan [3,28] and was approved by the USA Food and Drug Administration as the generic name, Nafamostat mesylate. VR23 is a small molecule that potentially inhibits trypsin-like proteasomes, chymotrypsin-like proteasomes and caspase-like proteasomes. Structurally, VR23 is a novel proteasome inhibitor with desirable anticancer properties that selectively kills cancer via cyclin E–mediated centrosome amplification without any notable ill effects [55]. The best binding poses for 3CLpro and PLpro with the top 2 ligands are shown in Fig. 1 & 2 , respectively.

Table 1

Table showing binding energy values and interactions of the top two ligands (protease inhibitors) with the key residues of 3CLpro and PLpro of the SARS-CoV-2 evaluated by PyRx docking.

Target protein	Inhibitors	Binding energy (kcal/mol)	Key residues interaction	Chain	H-bonds	Bond length (Å)
3CLpro	Nafamostat	−9.0	Ser1	–	N-N4	3.04
			Thr26	–	O-N2	2.93
					N-N2	3.14
			Glu166	–	OE2-N3	2.93
	VR23	−9.1	Ser1	–	N-O3	2.65
			Ser144	–	OG-O5	3.24
			Glu166	–	N-O1	2.90
PLpro	Nafamostat	−9.2	Trp106	C	N2-N	3.21
			Tyr273	B	N5-OH	2.91
			Asp286	C	N1-O	3.13
			Ala288	C	N2-N	3.10
			Asp302	B	N4-OD1	3.26
	VR23	−9.1	Trp106	C	O3-N	2.80
			Asn267	A	O5-ND2	3.13
			Tyr273	A	N5-OH	2.97
			Ala288	C	O3-N	2.93

Fig. 1

The docking results of (A) Nafamostat (Red) and (B) VR23 (aquamarine) inside binding pocket of the 3CLpro (cartoon and blue) of the SARS- CoV-2. Hydrogen bonded interactions are shown as yellow dotted lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

The docking results of (A) Nafamostat (Red) and (B) VR23 (aquamarine) inside the binding pocket of the PLpro (cartoon and magenta) of the SARS- CoV-2. Hydrogen bonded interactions are shown as yellow dotted lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The docking experiment was carried out in PyRx with grid box formation to include all the active site residues for SARS-CoV-2 3CLpro and PLpro. Further bond length, H-bond formation between ligands and target protein were analyzed using Ligplot and interactions were depicted in Fig. 3 A-D. Both the protease inhibitors (ligands) Nafamostat and VR23 did not form any H-bond with the catalytic dyad (Cys145 and His41) of SARS-CoV-2 3CLpro as revealed from the interaction studies, but both Nafamostat and VR23 revealed to form H-bond with the residues Ser1 and Glu166, which provides stability and active conformation of S1 pocket in SARS-CoV-2 3CLpro [67,77]. Apart from these two residues, Nafamostat also forms H-bond with Thr26, while VR23 forms H-bond with Ser144 (Fig. 3 A-B).

Fig. 3

Diagrammatic sketch illustrating the interactions between (A) Nafamostat and 3CLpro protein. (B) VR23 and 3CLpro protein. (C) PLpro protein and Nafamostat. (D) PLpro protein and VR23 by LigPlot. Ligand is shown in purple and: green dashed lines indicate hydrogen bonds with distance in angstrom (Å), spoked red arcs indicate hydrophobic contacts, atoms are shown in black for carbon, blue for nitrogen, red represents oxygen, green represents fluorine and yellow represents sulfur. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

On the other hand, the 2D interaction of the selected protease inhibitors with SARS-CoV-2 PLpro revealed Nafamostat's interaction and H-bond formation with the residue Asp286 of the catalytic triad Cys111-His272-Asp286 of subunit C apart from other vital interactions with the binding site residues across other subunits (Fig. 3C). VR23 did not form an H-bond with the catalytic triad Cys111-His272-Asp286 of any of the subunits of PLpro. However, VR23 was highly efficient in interacting with other binding site residues such as Trp106 and Ala288 of subunit C and Asn267 and Tyr273 of subunit A respectively, of PLpro (Fig. 3D).

Lipinski's rule of 5 (http://www.scfbio-iitd.res.in/software/drugdesign/lipinski.jsp) predicted the selected molecular properties of the ligands from the SMILES imported from PubChem (https://pubchem.ncbi.nlm.nih.gov/). Under Lipinski rule of 5, the calculation for drug-likeness was done with criteria namely: (M.W. <500 Dalton), high lipophilicity expressed as LOGP (LOGP<5), hydrogen bond donors (HBD <5), hydrogen bond acceptors (HBA <10) and molar refractivity should be between 40 and 130 indicated good absorption and permeation across the cell membrane. For an effective drug to be taken orally, it should satisfy at least four criteria out of the 5, thus distinguishing a molecule from a drug or a non-drug [[41]; http://www.scfbio-iitd.res.in/software/drugdesign/lipinski.jsp]. The selected inhibitors Nafamostat and VR23 both satisfied Lipinski's rule of 5 with LOGP value (<5), evaluating them with good permeability across the cell membrane [45] and establishing them as pharmacologically active (Table S2). As evaluated by molinspiration, the molecular property showed excellent results for both the selected ligands (protease inhibitors) to evaluate a single violation for both the inhibitors in Lipinski's rule of 5 (Table 2 ). The AB% (>50%) of a drug is a good sign for its excellent bioavailability, distribution and circulation by oral route [[18], [22]]. Both Nafamostat and VR23 possessed AB% > 50% (Table 2).

Table 2

Physicochemical proprieties prediction of the selected ligands (Nafamostat and VR23) using Molinspiration software.

Inhibitors	MW (g/mol) (<500)	MiLogP (≤ 5)	nON (acceptors) (<10)	nOHNH (donors) (<5)	Volume (A3)	TPSA (A2)	AB% (>50%)	Lipinski's violations ≤ 1
Nafamostat	347.38	2.16	7	7	306.84	140.59	60.49	1
VR23	477.89	3.58	11	0	368.86	145.16	58.91	1

The selected ligands' bioactivity was checked through molinspiration and calculated the activity against GPCR ligand, ion channel modulator, a kinase inhibitor, nuclear receptor ligand, protease inhibitor and enzyme inhibitor [45]. The bioactivity values were interpreted as: active (bioactivity score > 0), moderately active (bioactivity score: −5.0 – 0.0) and inactive (bioactivity score < −5.0) [62]. The bioactivity score of the selected inhibitors towards GPCR ligand, ion channel modulator, nuclear receptor ligand, kinase, protease and enzyme inhibitions indicated both Nafamostat and VR23 as active protease and enzyme inhibitors. The bioactivity prediction by molinspiration is tabulated in Table 3 .

Table 3

Bioactivity prediction of the Nafamostat and VR23 through Molinspiration software.

Inhibitors	GPCR ligand	Ion channel modulator	Kinase Inhibitor	Nuclear receptor ligand	Protease Inhibitor	Enzyme Inhibitor
Nafamostat	0.28	0.14	−0.03	−0.16	0.57	0.19
VR23	−0.05	−0.15	−0.04	−0.32	0.07	0.05

The pharmacodynamic study for the selected ligands was carried out through admetSAR to understand the drug's action inside a host's body. The ADMET properties of the selected ligands revealed through admetSAR (http://lmmd.ecust.edu.cn:8000/) are presented in Table 4 . The ADMET study in this work focused on the parameters such as solubility (LogS), human intestinal absorption (HIA), Caco-2 permeability, cytochrome substrate/inhibitor, human ether a go-go gene inhibition, AMES toxicity, carcinogens and acute rat toxicity (LD50). The result shows that Nafamostat and VR23 exhibit an ability to cross the blood-brain barrier (BBB), with a probability of 0.92 and 0.68. Besides, an excellent human intestinal absorption for selected ligands and a moderate ability to penetrate human colon adenocarcinoma (Caco-2+), with Caco-2 permeability of 0.53 cm/s and 0.57 cm/s for Nafamostat and VR23, respectively, were predicted (Table 4). These results confirm the high drug absorption (AB%) as predicted by Molinspiration [45] and further supported by the excellent values of HIA of 1.00 and 0.99 for Nafamostat and VR23 respectively in admetSAR (Table 4). LogS values of Nafamostat and VR23 were found to be −4.07 and −4.29, respectively, and were considered moderately soluble (Table 4).

Table 4

Table indicating pharmacodynamic profile for the selected ligands along with Remdesivir and Dexamethasone as control by admetSAR.

Inhibitors	Log S (>−4)	Blood Brain Barrier (BBB)	Human Intestinal Absorption (HIA)	Caco2 Permeability (cm/s)	CYP substrate / Inhibitor	Human ether a go-go gene	AMES toxicity	Carcinogenicity	LD50 (Rat Acute Toxicity) (mol/ kg)
Nafamostat	−4.07	0.92	0.99	0.53	Non-substrate/ Non-inhibitor	weak inhibitor	nontoxic	Non-carcinogen	2.56
VR23	−4.29	0.68	1.00	0.57	Substrate/ Inhibitor	weak inhibitor	nontoxic	Non-carcinogen	2.52
Remdesivir	−3.47	0.74	0.88	−0.14	Substrate/ Non-inhibitor	weak inhibitor	nontoxic	Non-carcinogen	2.71
Dexamethasone	−3.70	0.97	0.99	1.08	Substrate/ Non-inhibitor	weak inhibitor	nontoxic	Non-carcinogen	2.14

In terms of metabolism, Nafamostat was a non-substrate/non-inhibitor of cytochrome P450 (CYP 450), indicating its proper metabolism by CYP450, while VR23 was a substrate for CYP450 3A4 and inhibitor of CYP450 (Table 4).

Toxicity analysis predicted that both the selected protease inhibitors were non-AMES toxic and non-carcinogens. The Nafamostat and VR23 were also weak inhibitors of human ether, a go-go gene, a regulatory potassium channel that leads to long QT syndrome [60,[10], [11]]. Comparing the predicted LD50 doses, a compound with a lower dose is more lethal than the compound having a higher LD50 [47]. From our observation in admetSAR, we found that both the selected protease inhibitors have optimal LD50 doses indicating they are nonlethal and possess excellent pharmacodynamic properties (Table 4).

Due to the emerging and effectiveness of remdesivir and dexamethasone in treating COVID-19 patients in clinical trials [8,65,19,[25], [30]], the ADMET properties of the selected ligands were thus compared with it. All the selected ligands showed immensely satisfying results, with some of the values even better than the control drugs, as shown in Table 4 .

Molecular dynamics (MD) simulation studies

Molecular docking and pharmacokinetic analyses resulted in selecting two possible protease inhibitors, Nafamostat and VR23, against the SARS-CoV-2 protease proteins, 3CLpro and PLpro. In order to examine the molecular stability of the inhibitor's interactions with viral proteases, the best docking poses were selected and all-atoms MD simulations were carried out for the period of 100 ns at the physiological temperature. The changes in the structural features underlying the protease-inhibitors interactions were analyzed through investigating the plots of RMSD, Rg, SASA, RMSF, H-bonds and potential energies, as shown in Fig. 4 and 5 . The obtained backbone Cα-RMSD plot of inhibitors, Nafamostat and VR23 complexed with 3CLpro, during the simulation is depicted in Fig. 4 A. In this Fig. 4A, we can see that the RMSD plot of the 3CLpro-Nafamostat complex reaches the steady equilibrium at ~30 ns. The initial rising in RMSD of ~0.15 nm during 0–25 ns indicates the structural perturbation to well-fit the inhibitor, which is gradually dropped at ~30 ns. The complex structure remains stable with RMSD 0.19±0.02 nm till the end of simulation at 100 ns. The RMSD plot of 3CLpro-VR23 shows that it also attains equilibrium at ~30 ns, after the initial deviation of RMSD ~0.10 ns; however, the structure is observed a stable RMSD value of 0.21±0.03 nm. The RMSD plot of the PLpro-Nafamostat complex suggested relatively stable conformational dynamics of protease-inhibitor, showing a steady equilibrium at RMSD 0.14±0.02 nm can be seen throughout the simulation time (Fig. 5 A). However, the binding of VR23 with PLpro shows the sharp drift of 2.0 nm, optimized at ~40 ns. This equilibrium state is prolonged up to 80 ns and the further increase in RMSD can be seen for the last 20 ns. The result shows that the structure of the PLpro-VR23 complex propagated with an average RMSD 1.89±0.60 nm. These findings indicate a stable global conformational dynamic of Nafamostat with both proteases, 3CLpro and PLpro, and VR23 with 3CLpro. Whereas the binding of VR23 with PLpro underwent significant conformational changes.

Fig. 4

The 100 ns molecular dynamic simulation results of two protein-ligand complexes (3CLpro- Nafamostat (black) and 3CLpro-VR23 (red)). (A) RMSD values of backbone atoms. (B) Rg of backbone atoms. (C) SASA of the ligand. (D) RMSF values of the chain. (E) Total number of H-bonds between ligands and 3CLpro. (F) Gromacs/ Potential energies values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

The 100 ns molecular dynamic simulation results of two protein-ligand complexes (PLpro- Nafamostat (blue) and PLpro-VR23 (pink)). (A) RMSD values of backbone atoms. (B) Rg of backbone atoms. (C) SASA of the ligand. (D) RMSF values of the chain. (E) Total number of H-bonds between ligands and PLpro.(F) Gromacs/ Potential energies values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Further, protease-inhibitors' conformation stability was determined by the radius of gyration (Rg), which defines the structural compactness of protein. The average Rg values for 3CLpro-Nafamostat, 3CLpro-VR23, PLpro-Nafamostat and PLpro-VR23 complexes were calculated as 2.19, 2.20, 2.32 and 3.32 nm, respectively (Fig. 4 B & 5 B). The Rg plot shows a minor increase in Rg values at 0–10 and 20–25 ns, which gradually settles around 30 ns in the Nafamostat docked structure with 3CLpro (Fig. 4 B). The initial deviation in Rg values may be considered the time taken to optimize Nafamostat's well-fitting at the binding pocket of 3CLpro. After that, no structural drifts were observed in the Rg trajectory of 3CLpro-Nafamostat, suggesting the complex's stability for the remaining period of simulation 30–100 ns. Protease-inhibitor interactions. The Rg plot of PLpro-Nafamostat quickly attains the stable equilibrium and is observed consistently throughout the simulation time (Fig. 5 B). The slight deviation in the Rg plot of 3CLpro-VR23 can be seen during the 0–30 ns, which attained equilibrium at ~30 ns and the structural integrity of the complex is observed stable for the simulation time 30–100 ns. Remarkably, the Rg plot of PLpro-VR23 shows a sharp drift of ~0.75 nm at ~15 ns, which settles gradually around ~40 ns and a stable Rg trajectory is seen up to 100 ns. The substantial jump in Rg values may be due to the significant conformational changes before accommodating the inhibitor at the binding pocket of PLpro. Thus, the significant deviation in Rg values and a slow rate to achieve equilibrium suggested a less stable conformational dynamics of PLpro-VR23 than the binding of VR23 with 3CLpro. Moreover, a better equilibrated and stabilized structure of Nafamostat with both 3CLpro and PLpro, respectively.

The changes in protease-inhibitor complexes' structural features have also been analyzed by computing the solvent-accessible surface area (SASA). The solvent environment around the protein act as a driving force for maintaining the protein folds. SASA is one of the fundamental properties of a protein that maintains the practical orientation of protein and accompanies protein-ligand interactions [50,54]. Thus, it helps to evaluate protein-inhibitor complexes' structural stability under the solvent environment [10]. Fig. 4C shows that the complex structure of 3CLpro-Nafamostat and 3CLpro-VR23 have average SASA values 180.82 ± 1.12 and 189.49 ± 1.14 nm2, respectively. The binding of Nafamostat and VR23 with PLpro has the average SASA values 177.41± 1.15 and 356.70 ± 1.62 nm2, respectively (Fig. 5 C). Thus, comparing the binding of inhibitors with proteases, it is observed that Nafamostat has stable molecular interaction with 3CLpro and PLpro. Eventually, the binding of VR23 with 3CLpro also provides evidence of stable molecular interactions. However, the molecular interaction of VR23 with PLpro indicates the orientational change in the protein surface; thus, we observed a higher value for the accessible area of PLpro-VR23.

To examine protease-inhibitor interactions' dynamic progression, we calculated the root mean square fluctuation (RMSF), which provides a clue about the average positional fluctuations of each amino acid residue [35,46]. Usually, the higher atomistic fluctuations are possessed by the N-and -C terminal residues and loops of protein; however, larger loops exist, interconnecting the stable conformation of α-helices and β-sheets may deviate the RMSF plots [34]. Fig. 4 D shows that the binding of inhibitors with protease 3CLpro complies with a stable molecular interaction and displays almost completely overlapping RMSF plots for 3CLpro-Nafamostat and 3CLpro-VR23. The much higher deviation in the RMSF plot of PLpro-VR23 can be observed compared to PLpro-Nafamostat (Fig. 5 D). The residues of 3CLpro belonging to the stable secondary structure of α-helices and β-sheets display average atomic fluctuations <0.25 nm and no significant perturbation is observed for the region belonging to loops. Indeed, the binding pocket residues (Ser1, Glu166, Thr26 & Ser144) having lower atomic fluctuation, which evident in the stable molecular interaction of inhibitors with 3CLpro. We also find the consistent molecular interaction of Nafamostat with PLpro(Fig. 5 D). However, the RMSF plot of PLpro -VR23 shows the overall higher atomic fluctuations for all residues. Interestingly, the lower atomic fluctuation for residues (Trp106, Ala288, Asn 267 &Tyr273), involved in the interaction with VR23, clearly indicates that ligand remains occupied at the binding pocket during the simulation and higher RMSF is observed due to the conformational shifting of protease, PLpro.

More importantly, we computed the time evolution plot of hydrogen bonds' occupancy (H-bonds) between target proteases and inhibitors. H-bonds play a crucial role in maintaining the shape and stability of protein structure and govern the protein-ligand molecular interactions. Fig. 4 E shows that the maximum propensity of three H-bonds between 3CLpro and VR23. Among them, only one H-bond remains consistent and two appeared and disappeared transiently. The molecular interaction of 3CLpro with Nafamostat displays two H-bonds, which can be seen up to ~65 ns. However, one H-bond is lost during the last ~35 ns of simulation. The H-bond interactions of inhibitors with PLpro is shown in Fig. 5 E. Fig. 5E shows Nafamostat's structure is stabilized at the binding pocket of PLpro through four H-bonds, out of which two are observed consistently during the simulation.

On the other hand, we find the propensity of single H-bond interaction between VR23 and PLpro, which appears intermittently. Thus, H-bond interactions, along with structural parameter analyses described in Fig. 5A-D, provide clear evidence of the loosely bound structure of VR23 with PLpro. Finally, we measured the potential energy of inhibitor bound complexes of 3CLpro and PLpro. Fig. 4 F shows that the structure of 3CLpro-Nafamostat and 3CLpro-VR23 have potential energy values −7.60 × 105 and −7.87 × 105 kJ/mol, respectively. Whereas the structure of PLpro- Nafamostat and PLpro-VR23 shows the potential energy values −1.30 × 105 and −1.03 × 105 kJ/mol (Fig. 5 F). Thus, the obtained higher potential energy values for inhibitors bound complexes of 3CLpro compared to complexes with PLpro suggested the higher molecular affinity of Nafamostat and VR23 towards 3CLpro. This result also indicates Nafamostat's higher molecular affinity with PLpro than VR23, showing approximately 2.7 × 105 kJ/mol higher potential energy.

The MD simulation analyses collectively revealed a higher structural stability of Nafamostat with both proteases, 3CLpro and PLpro. Moreover, VR23 shows preferably stable molecular interaction towards 3CLpro as compared to PLpro. This study provides the structural insights of the molecular interaction of identified inhibitors, Nafamostat and VR23, against the two significant proteases of SARS-CoV-2, 3CLpro and PLpro, which can be explored as a potential lead in the development of anti-CoV drug therapy to control the pandemic of COVID-19.

Conclusion

Our study revealed that Nafamostat and VR23 interact with the critical binding site residues of 3CLpro and PLpro of SARS-CoV-2 and probably inhibit the growth of SARS-CoV-2 by targeting 3CLpro and PLpro. Pharmacokinetics properties matching the required criteria for drug-likeness and Nafamostat and VR23 may be considered a potential drug against SARS-CoV-2. Protein-ligand complexes for 3CLpro and PLpro were structurally stable throughout the 100 ns simulation period concerning distance and fluctuation dynamics. However, further experimental studies are needed to check the possible preclinical and clinical efficacy of Nafamostat and VR23 to prevent and treat SARS-CoV-2.

Author contribution statements

DK, RDS, AP and DB carried out the experiment. DB, AP and DK wrote the manuscript. DK, AP and DB contributed to the analysis of the results. DK supervised the project and conceived the original idea.

Declaration of Competing Interest

The authors declare that no conflict of interest exists.