Discovering potential inhibitors against SARS-CoV-2 by targeting Nsp13 Helicase

Abstract The rise in the incidence of COVID-19 as a result of SARS-CoV-2 infection has threatened public health globally. Till now, there have been no proper prophylactics available to fight COVID-19, necessitating the advancement and evolution of effective curative against SARS-CoV-2. This study aimed at the nonstructural protein 13 (nsp13) helicase as a promising target for drug development against COVID-19. A unique collection of nucleoside analogs was screened against the SARS-CoV-2 helicase protein, for which a molecular docking experiment was executed to depict the selected ligand’s binding affinity with the SARS-CoV-2 helicase proteins. Simultaneously, molecular dynamic simulations were performed to examine the protein’s binding site’s conformational stability, flexibility, and interaction with the ligands. Key nucleoside ligands were selected for pharmacokinetic analysis based on their docking scores. Selected ligands (cordycepin and pritelivir) showed excellent pharmacokinetics and were well stabilized at the proteins’ binding site throughout the MD simulation. We have also performed binding free energy analysis or the binding characteristics of ligands with Nsp13 by using MM-PBSA and MM-GBSA. Free energy calculation by MM-PBSA and MM-GBSA analysis suggests that pritelivir may work as viable therapeutics for efficient drug advancement against SARS-CoV-2 Nsp13 helicase, potentially arresting the SARS-CoV-2 replication. Communicated by Ramaswamy H. Sarma


Introduction
Till July 2021, the COVID-19 outbreak has claimed more than 4,229,249 lives, establishing it as a global threat to the human population, jeopardizing social and economic life worldwide. It all began in the Wuhan provinces of China in December 2019, where the first known human infection was reported .
Though several potential drug and vaccine candidates are already approved for emergency use, their potency and efficacy data are still in a preliminary stage to recognize them as a guaranteed treatment for COVID-19. Thus, repurposing drugs can be a short time alternative and cost-effective compared to de novo drug discovery to control the outbreak.
Coronavirus (CoV) is an enveloped, single-stranded, positive-sense RNA virus, has a 5 0 -capping and 3 0 -polyadenylated tail. It belongs to the family Coronaviridae, which can be further subdivided into a, b, gamma, and delta coronavirus (Fehr et al., 2015;Perlman et al., 2009). The closely related SARS-CoV and MERS-CoV, which shares a close genomic relationship with the SARS-CoV-2 with approx 79.5% and 54% sequence similarity, respectively, gave us the right level of understanding of the virus, the immunity response (host) to the infection, and potential therapeutic targets . These studies revealed lucrative approaches to target key enzymes such as protease, polymerase, and helicases (Katsiki et al., 2020). Thus, recognizing the coronavirus genus conserved sequences became a work of utmost importance. The nonstructural protein 13 (Nsp13) or helicase of SARS-CoV-2 is requisite for RNA replication and an attractive drug target for antivirals development.
Comparing the protein sequence of SARS-CoV and SARS-CoV-2 has revealed the helicase gene (nsp13) to be highly conserved. Helicase protein is also well conserved among a, b, and gamma coronaviruses (Adedeji et al., 2014). Current studies suggested the conserved amino acid residues belonging to the active site domain of NTPase, could be explored as a promising molecular target. In another study on the Nsp13 of 2019-nCoV, it is well demonstrated that helicase of SARS-CoV-2 catalyzes the unwinding of nucleic acid into single strands from 5 0 -3 0 direction in an NTP-dependent manner (Jia et al., 2019;Singleton et al., 2007). Helicase of SARS-CoV-2 was also found to target the host's innate immune signaling proteins (Gordon et al., 2020). Thus, the inhibition of helicase protein may directly lead to disrupting other proteins like protease (Gordon et al., 2020). Considering all these vital roles, we have specifically targeted the critically important Nsp13 helicase enzyme to identify some promising inhibitors that can act as antivirals and help eliminate the deadly SARS-CoV-2 in the future.
Recent investigations revealed that some small molecular weight compounds could restrain and inhibit NTPase activity and its associated binding site for ATP (Mirza et al., 2020). Other studies also showed nucleoside analogs' potential to inhibit NTPase/helicase protein (Pruijssers et al., 2019;Te et al., 2007;Yedavalli et al., 2008;Zhang et al., 2003). With this in mind, the present study emphasizes screening a unique collection of nucleoside analogs against Nsp13 of the SARS-COV-2 in molecular docking and dynamic settings.

Protein target selection and preparation
The Crystal Structure of Nsp13 protein (PDB ID: 6ZSL) (Newman et al., 2020) was retrieved from the Protein Data Bank (http://www.pdb.org). Water molecules, hetero atoms, and co-crystallized ligands (PO 4 and ZN) were removed and subsequently, the protein was energy minimized by the YASARA Energy Minimization Server (Krieger et al., 2009). The docking-ready protein structure was first converted into a pdbqt file and then was subjected to the Graphical User Interface program "Auto Dock Tools (ADT)," where the chain A of the finalized protein (6ZSL) was prepared (Huey et al., 2008).

Prediction of the binding site, ligand selection and ligand file preparation
Binding sites were selected based on the previous studies (Mirza et al., 2020;Pal et al., 2020). A unique collection of 135 nucleoside analogs from the Selleckchem database (https://www.selleckchem.com/) was downloaded in sdf file (MDL MOL format) along with the two viral helicase-primase complex inhibitors Amenamevir and Pritelivir (Kawashima et al., 2017;Shiraki et al., 2017;Wald et al., 2014;Wozniak et al., 2013). The ligands were selected based on the recent studies evaluating the Nsp13 of 2019-nCoV possessing the helicase and NTPase activity in the presence of NTP (Shu et al., 2020). Previous studies put forward a stronghold on nucleoside analogs as they change the genetic foundation of viruses by activating or inhibiting the unwinding process of NTPase/helicase, thus leading to a decline in the viral load with each replication cycle (Pruijssers et al., 2019;Te et al., 2007;Yedavalliet al., 2008;Zhang et al., 2003). The PubChem database detailed all the compounds' chemical information, thus confirmed and saved in smi format converted to pdbqt and mol file by Open Babel (Boyle et al., 2011) for docking studies.

Molecular docking and interaction analysis
For docking purposes, two software, PyRx and iGEMDOCK (Dallakyan et al., 2015;Hsu et al., 2011), applying two separate algorithms were utilized. Using the PyRx platform, which provides the graphical interface for the docking program AutoDock Vina (Trott et al., 2010), the macromolecular structure of SARS-CoV-2 Nsp13 protein and ligands were processed, the grid center was set at X¼ À15.85, Y ¼ 6.79, and Z¼ À76.03 with the grid box parameters of 23.28 Å Â 20.95 Å Â 13.26 Å for X, Y, and Z co-ordinates respectively and then, the docking was executed. GEMDOCK utilizes the evolutionary genetic method (Yang and Chen, 2004), providing the interface for protein-ligand interaction and dock through its iGEMDOCK tool. The docking protocol in the iGEMDOCK tool was set to a "Stable docking" procedure with a population size of 80 generations (N ¼ 800) and ten solutions were accomplished (Hsu et al., 2011). Interactions using Ligplot (Wallace et al., 1995) and visualization through Pymol Molecular Visualization Software (The PyMOL Molecular Graphics System, Version 2.0 Schr€ odinger, LLC) were studied between the target protein and the inhibitors.

Molecular dynamics simulation
Using the biosimulation package Amber16 (Case et al., 2005), all atoms molecular simulation was applied on the coordinates of Nsp13 and docked complexes with ligands, Cordycepin, and Pritelivir. The force field was selected as ff14SB, and the water model TIP3P was used to solvate all three systems (Jorgensen et al., 1983;Maier et al., 2015). With the help of the Antechamber tool, the ligands parameters were defined. AM1-BCC was used to calculate charges for ligands parameters (Wang et al., 2006), and the force field was taken as GAFF (Jakalian et al., 2002;Souda da Silva & Vranken, 2012). The Amber16 tool, tleap, defines the simulation boxes of dimension octahedral, having a buffer distance of 10 Ð. The prepared system was neutralized with counter ions, 0.15 M of Na þ and Cl - (Joung & Cheatham, 2008). SHAKE algorithm was applied to restrain the bonds of hydrogen atoms. The long-range electrostatic was handled with Particle mesh Ewald (PME) (Darden et al., 1993). The temperature and pressure during the simulation were maintained by the Langevin thermostat and Berendson's barostat, respectively (Berendsen et al., 1984). The energy minimization was performed in two phases, 5000 steps of each, using the steepest gradient and conjugate gradient algorithm. After the minimization, the ensemble process NVT and NPT were applied to equilibrate prepared systems as described else (Mishra et al., 2020). The pmemd.cuda was used for the production run of 100 ns on the NVT ensemble, with time step 2 fs. All files, trajectories, velocity and energy, were saved at a gap of every ten ps. The Amber16 tool, cpptraj, was used to analyze the simulation trajectories (Roe & Cheatham, 2013).

Binding free energies analysis
The binding free energies of ligands with Nsp13 were calculated using MM-PBSA, which estimates binding free energy components applying the combination of molecular mechanics and continuum solvent models, which accounts for the structural stability of ligands at the binding pocket of the protein. This approach considers the solute molecule as a dielectric object whose shape is determined by the atomic coordinates and their radii. The Poisson-Boltzmann (PB) equation was applied to the electrostatic potential term. The nonpolar part of solvation-free energy G NPOL was calculated through the solvent-accessible surface area (SASA) of the solute. Another implicit solvation model, molecular mechanics generalized Born surface area (MM-GBSA), was also applied to calculate the binding free energy of the ligand-bound complexes of Nsp13, which calculates the polar or electrostatics (G SOL-GB ) based on generalized Born expression (Genheden & Ryde, 2015;Wang et al., 2019). Using Amber tool mmpbsa.py (Miller et al., 2012), MM-PB(GB)SA can estimate the binding free energy (DG bind ) as, DE MM denotes the nonbonded and bonded interaction energies. DG Solv defines the solvation free energy (polar and nonpolar contributions of solvation) and the polar solvation term calculated using the Generalized-Born model or Poisson-Boltzmann and the nonpolar term estimated from SASA. Excluding the entropy term (TDS), which is computationally expensive, the above equation for the binding free energy can be approximately written as, Where DE MM is the change in the average molecular mechanics' interaction energy (gas-phase) upon ligand binding computed as the sum of the changes in the bonded and nonbonded (electrostatics and Van der Waals) interactions upon ligand binding ( , DG solv is the change in solvation free energy upon ligand binding. For binding free energy calculation using MM-GBSA, no dielectric constant has been used for protein interior (i.e. in vacuum or dielectric constant 1), while a dielectric constant of 4 has been used in MM-PBSA methodology.

Protein binding site prediction
The active/binding site residues of SARS-CoV-2 Nsp13 were predicted from the literature (Mirza et al., 2020;Pal et al., 2020). Large pocket or active site selection in the docking process increases the chance of ligand binding randomly ( Figure 1 and Supplementary Table S1).

Molecular docking and interaction analysis
The most challenging computational chemistry task is to foretell the substantial interacting arrangement of inhibitors into protein active sites (Azam et al., 2014). The docking analysis carried out with the 135 nucleoside analogs from Selleckchem.com, Bioactive Compound Expert and the two helicase-primase inhibitors Amenamevir and Pritelivir showed some excellent results. The active/binding site residues of SARS-CoV-2 Nsp13 were retrieved from the works of literature ( Figure 1 and Supplementary Table S1). Both AutoDock Vina and GEMDOCK results show Cordycepin and Pritelivir are the best ligands against Nsp13 of SARS-CoV-2. Using AutoDock Vina, the docking score of ligands, Cordycepin and Pritelivir, obtained À9.0 kcal/mol and À8.8 kcal/mol, respectively (Table 1), whereas GEMDOCK shows binding scores, À112.54 kcal/mol and À126.61 kcal/mol, respectively (Table  2). Cordycepin (3 0 -deoxyadenosine) is a natural adenosine analog and a biologically active component isolated from Cordyceps militaris and is also found in other Cordyceps species Jeong et al., 2011;Li et al., 2015;Zhou et al., 2008). Intracellular targets for cordycepin include nucleic acid (DNA/RNA) and apoptosis (Tuli et al., 2013) and possesses anti-inflammatory, antioxidant, immunomodulatory, and anticancer activities (Soltani et al., 2019). While, Pritelivir (BAY 57-1293) showed antiviral activity performed on animal models of herpes simplex virus disease (Wald et al., 2014;Wozniak et al., 2013).
Interaction studies showed the H-bond formation of the selected ligands cordycepin and pritelivir with the Nsp13 binding site residues, signifying their efficient binding to Nsp13, probably restricting the helicase activity and the unwinding of DNA. Ligplot analysis of the nucleoside analog, cordycepin revealed the strong H-bond formation with the binding site residues Lys288, Ser289, Gln404 with Nsp13 protein. Besides the binding residues of Nsp13 of SARS-CoV-2, the cordycepin also formed H-bonds with other residues. Simultaneously, the viral helicase-primase inhibitor pritelivir showed H-bond interactions with the binding site residues Glu375 and Gln404. All the interaction figures were visualized in PyMol (The PyMOL Molecular Graphics System, Version 2.0 Schr€ odinger, LLC). The best docking poses for selected ligands are shown in Figure 2(a,b). The 2D interplay of the top ligands with the Nsp13 protein using Ligplot is shown in Figure 3(a,b).

In silico pharmacokinetics studies
Prediction and calculation for drug-likeness through Lipinski rule of five was performed with standards, viz: Molecular Weight (MW <500 Dalton), lipophilicity (LogP <5), hydrogen bond acceptors (HBA <10), hydrogen bond donors (HBD <5) and molar refractivity range within 40-130 showed proper absorption or permeation across the cell membrane. A drug should fulfill four criteria to be considered active and taken orally (Lipinski, 2004). Molsoft (http://www.molsoft.com/ mpropdesc.html) applies fragment-based contributions and, based on molecular features, renders a drug score for the ligands under consideration. Potential drug candidates should possess a drug score greater than zero (>0) (Chandrasekaran & Kumar, 2016). Selected ligand pritelivir satisfied Lipinski's rule of 5 with impressive drug scores of 0.98. Simultaneously, cordycepin scored 0.49, which is quite good, though not too impressive as the former. The obtained results show excellent physicochemical properties, suggesting them as pharmacologically active compounds (Table 3 and Table S2; Figure 4(a,b)).
Furthermore, SwissADME revealed that cordycepin and pritelivir have good oral bioavailability with a score of 0.55 (Table 3). The results from the five different methods: Lipinski (Pfizer) filter, Ghose filter, Veber (GSK) filter, Egan (Pharmacia) filter and Muegge (Bayer) filter are summarized in Table 3. These results were in good agreement for the selected ligands with those obtained in the "pink area" of the bioavailability radar plot ( Figure 5(a,b)).
This plot gives a precise figure of the drug-likeness criterion of a bioactive drug to be taken orally. A hexagonal radar plot is conferred with the vertices representing a single criterion for an orally bioactive drug (Daina et al., 2017). The "pink area" inside the hexagon denotes the optimum dimension for lipophilicity (À0.7 to þ5.0), MW (150 to 500 g/mol), polarity (TPSA from 20 to 130 A2), solubility (<6), saturation (carbons in the sp3 hybridization not less than 0.25) and flexibility (<9 rotatable bonds).
Complimenting the above drug-likeness prediction of the small molecules, PAINS, Brenk and Lead likeness screening evaluated by SwissADME was carried out. PAINS (Baell JB et al., 2010) is a high throughput screening model that identifies or predicts compounds highly engaged in various biological assays and is frequently published in the literature. No PAINS alerts were evaluated in SwissADME for the selected ligands (Table 3). Brenk et al. (2008) considered smaller and less hydrophobic compounds and not the "Lipinski's rule of five" evaluated compounds for widening hit optimization opportunities. Lead likeness tests (Teague et al., 1999) in high-throughput screening are anticipated to specify leads with substantial efficiency and enable the detection and development of novel interactions to discover and develop a new drug. None of the molecules flouted the lead likeness test (Table 3).
The ADMET properties through admetSAR of the selected ligands are shown in Table 4. The ADMET analysis concentrated on the properties such as solubility (LogS), human intestinal absorption (HIA), Caco-2 permeability, cytochrome substrate/ inhibitor, AMES toxicity, carcinogens and acute rat toxicity (LD50) for understanding the action of the selected ligands inside the host's body. The selected ligands possessed excellent HIA and Caco-2 permeability values and are expected to be highly absorbed upon oral administration (Table 4) (Sood et al., 2018). Cordycepin showed a better absorption rate with the highest HIA of 0.99 (Table 4). All the selected ligands exhibited LogS values higher than À4 (> À4) and were considered optimally soluble (Table 4) (Mokhnache et al., 2019).
In terms of metabolism, cordycepin was a non-substrate/ non-inhibitor of cytochrome P450 (CYP 450). At the same   (Table 4). Inhibition of CYP450 will lead to the drug's improper metabolization by the CYP 450 enzyme and, consequently, hamper the biological transformation of the drug inside the host's body (Nisha et al., 2016). Moreover, no AMES toxicity and carcinogenic properties were related to our selected ligands ( Table 4). The selected ligands, cordycepin and pritelivir, showed optimal LD50 doses in admetSAR (Table 4). It was indicative of their nonlethal character, as a compound with a higher LD50 dose is generally considered to be less lethal than a compound with a lower dose .
Together, all the selected ligands' ADMET parameters were compared with the controls, i.e. remdesivir and dexamethasone (Table 4). Remdesivir (GS-5734), a nucleoside analog, was used for treating COVID-19 in multicentre trials for its inhibition of pathogenic human as well as animal coronaviruses, including in vitro inhibition of SARS-CoV-2 in animal models (Beigel et al., 2020;Goldman et al., 2020;. Whereas dexamethasone, an inexpensive and widely ready to use steroid, lower the death rates among critically COVID-19 patients with good results in the preliminary clinical trials (Horby et al., 2020;Johnson et al., 2020). Our selected ligands showed excellent values and results compared with the control, which agrees that our selected ligands, cordycepin and pritelivir, possess excellent pharmacodynamic properties.
Mutations had already imposed resistance to nucleoside analogs on the viral RdRp in coronavirus along with other RNA viruses (    effective anti-CoV nucleoside analogs could restrict viral resistance and enhance potency (Pruijssers et al., 2019). Helicase inhibitors, protease inhibitors, monoclonal antibodies, viral entry inhibitors, and DEDDh family exoribonuclease inhibitors are compounds for effective treatment combinations (Huang et al., 2016;Liang et al., 2018;Zumla et al., 2016). Thus, nucleoside analog cordycepin is evaluated as a potential SARS-CoV-2 Nsp13 inhibitor and can also be used in combination with the viral helicase-primase inhibitor pritelivir better efficacy and treatment.

Conformation stability of protein-ligand complex
The efficiency of a drug molecule depends on the effective binding of ligands at the binding site of the target protein  . The bioavailability radar of the selected ligands was evaluated using the swissADME web tool. Cordycepin (a) and Pritelivir (b). The pink area represents the optimal range for each property (lipophilicity: XLOGP3 between À0.7 and þ5.0, size: MW between 150 and 500 g/mol, polarity: TPSA between 20 and 130 Å 2 , solubility: log S not higher than 6, saturation: fraction of carbons in the sp 3 hybridization not less than 0.25, and flexibility: no more than nine rotatable bonds. Cordycepin is predicted best orally bioavailable ligand than pritelivir and fludarabine phosphate.  Majewski et al., 2019;Mishra et al., 2018;Mishra et al., 2020;. Hence, to investigate the conformational stability of protein-ligand complexes, MD simulations were done in the aqueous environment for 100 ns at 300 K using the parameters as described . To evaluate the structural order parameters of Nsp13 binding with antiviral drug molecules, we examine the RMSD, Rg, SASA, and RMSF (Prakash & Luthra, 2012). The time evolution plot for root means square deviation (RMSD) of Ca-atoms backbone is shown in Figure  6(a), indicating that the conformational dynamics of Nsp13 settles quickly around five ns and remains consistent up to 30 ns. We can see the drift of $1.0 Å at $30 ns which gradually dropped to 40 ns. The successive increase in RMSD can be noticed during 55-75 ns, and it settles at 80 ns. Similarly, the slight rise in RMSD is noted at 95-100 ns. Remarkably the binding of Pritelivir results in the more stable structure dynamics of Nsp13. The RMSD plot of the Nsp13-Pritelivir complex shows that it quickly attains equilibrium at $5 ns, and a pretty steady trajectory can be seen for the period of 5-100 ns. In contrast, the RMSD trajectory of the Nsp13-cordycepin complex shows higher conformational fluctuation with the drift of $1.5-2.5 Å during the initial 0-45 ns. The trajectory settles at $45 ns, and a comparatively stable structure is observed up to 70 ns. With the slight drift of $1.0 Å at $75 ns, the structure stabilized at an RMSD value of $5.0 Å. Thus, RMSD results indicate the more stable conformational dynamics of the Nsp13-Pritelivir complex compared to the Nsp13 and Nsp13-cordycepin complex. We also measured the ligand RMSD, which shows the initial perturbation required to spatially well settled at the active site of Nsp13. The trajectory of both ligands at $30 and $50ns; however, the trajectory of Pritelivir shows relatively stable equilibrium during $65-100 ns, whereas cordycepin achieved equilibrium at $30 ns, but the slight perturbation in trajectory can be seen up to the end of simulation at 100 ns (Supplementary Figure S1). To examine the conformational dynamics of protein-ligand complexes, we computed the radius of gyration (Rg) plots that described the protein's structural compactness (Figure 6(b)). Results show that the Rg trajectories of Nsp13 and Nsp13-Pritelivir complex quickly achieves equilibrium during the initial $5.0 ns and remain consistent with the marginal perturbation at $70-80 ns, respectively. In comparison, the trajectory of Nsp13 complexed with Cordycepin shows the sharp drifts of R g approximately 2.0-3.0 Å during the initial 0-25 ns. A continuous dropdown in the plot can be seen during the $25-50 ns. It obtained equilibrium at $50 ns that remains consistent up to 85 ns; after that, with a bit of dropdown, the simulation ends with the average R g value of 35.89 ± 0.02 Å. The average change in Rg value of Nsp13-Pritelivir (32.47 ± 0.01 Å) and Nsp13 (34.83 ± 0.02 Å) suggest intact structural integrity of Nsp13 complexed with Pritelivir as compared to Cordycepin. Furthermore, to investigate the molecular binding stabilities of ligands at the active site Nsp13, we measured the average distance of ligands from the center of the binding pocket. Figure 6(c) shows that during the initial 0-30 ns, the average distance of Pritelivir from the center of the active site remains around $3.0 Å. Drifts in trajectory seen at $35 ns which may be due to spatial adjustment of ligand at the active site of Nsp13; however, the center of ligand remains consistent at $5.0 Å, which suggest the stable molecular interaction between Pritelivir-Nsp13 for the remaining period of simulation $40-100 ns. The average distance of Cordycepin was also observed around $4.0 Å during the initial 0-25 ns. However, it shows a sharp drift at $25 ns which is stabilized at $30 ns, and a consistent interaction observed up to $60 ns. After that, the sudden drop in trajectory followed by an abrupt increase in the distance during $75-100 ns indicates that Cordycepin moves out from the active site of Nsp13. This may be one of the reasons we observed inconsistent trajectory of RMSD and Rg for the Nsp13-Cordycepin complex.
Next, to examine the involvement of binding pocket residues in the spatial stability of ligands, RMSF was calculated, which defines the average atomic fluctuation of amino-acid residues. Usually, the residues belonging to stable conformations of a-helices and b-sheets indicate the lower RMSF value. The residues at the N-and C terminal of protein and loops show high RMSF values. The secondary structure of Nsp13 consists of eleven a-helices and twenty-three b-sheets which shows the stable atomic fluctuation; however, the loops interconnecting the stable structure of a-helices and b-sheets, having residues length >10 show higher atomic fluctuation. Figure 6(d) shows high mobility of residues at Nand C-terminals, along with the loops (Val49-Gln62, Leu227-Gln243, Val247-Val266, Asp315-Glu353, Asp401-Phe422, Lys430-Ala448, Val452-Gln470, Met474-Asn516, Val521-Tyr543, Ile575-Asn596). The average atomic fluctuations of residues < 1.5 Å belong to the stable conformations of a-helices and b-sheets. The active site of Nsp13 is composed of a4 and a5, capped with small loops: Gln281-Ser288 and Tyr299-Arg303. The four parallel b-sheets: b17-b20 and long loop: Aps401-Phe422 formed the shallow hydrophobic cavity. This is why we observed the lower average atomic fluctuations of residues: Aps401-Phe422, even though belonging to loop. Moreover, the progressive increase in mobility of the residues belongs to five loops between Lys430-Asn596 at the C-terminal of Nsp13. Whereas the residues actively involve molecular binding of ligands shows significantly low RMSF for Nsp13-Pritelivir compared to Nsp13-cordycepin. Indeed, the RMSF plot of Nsp13-Pritelivir shows fewer average atomic fluctuations, which are evident in the stable binding of Pritelivir. The secondary structure analysis applying DSSP also suggests no structural distortion during the simulation, as shown in Supplementary Figure S2.

H-bond occupancy
Hydrogen bonds (H-bonds) interactions play a crucial role in the molecular stability of protein structure and the stable spatial binding of ligands at the active site of proteins. Figure 7 shows the H-bonds occupancy of ligands, pritelivir, and Cordycepin with Nsp13. Pritelivir shows the maximum occupancy of four H-bonds; however, only two H-bonds were observed consistently throughout the simulation time (Figure 7(a)). The binding of Cordycepin at the active site of Nsp13 shows the maximum possibility of six H-bonds interaction observed during the 0-20 ns. However, the four Hbonds were disrupted gradually at $20-30 ns (Figure 7(b)). We can see the sudden gain in H-bonds at 30-35 ns which again dropped to three H-bonds and remain consistent up to 60 ns. The occurrence of intermittently loss and gain of Hbonds was noted during 60-100 ns. Thus, these results suggest that both ligands, Pritelivir and Cordycepin, stabilized through the H-bond interaction at the binding pocket of Nsp13.

Essential dynamics
The dynamics of Nsp-13, Nsp13-Pritelivir, and Nsp13-cordycepin were analyzed through the principal component analysis (PCA) in the conformational space. PCA characterizes the essential dynamics (ED) governing the conformational changes in the protein during MD simulations. This is calculated by projecting the first two principal components, PC1, and PC2, for the backbone atoms of protein to adequately describe conformational redistributions induced by binding events, as plotted in Figure 8(a,b). From this Figure 8, it is apparent that the collective motion of Nsp-13 navigates the ample conformational phase space represented by the broad distribution of components along the PC1 and PC2. The highly scattered plot in the broad region of phase space indicates the flexible conformational dynamics of Nsp13-cordycepin. However, the confined collective motion of Nsp13- Pritelivir to a small conformational phase space confirms that the overall decreased in the structural flexibility of Nsp-13 on the binding of Pritelivir. Thus, the coherent change in the ED of Nsp13-Pritelivir again provides the elegance evidence of stable molecular interaction of Pritelivir with Nsp-13 compared to Cordycepin.

Binding free energy estimation: MM-PBSA and MM-GBSA
Finally, for the quantitative assessment of the binding affinity of ligands Pritelivir and Cordycepin to Nsp13, we computed the binding free energy (DG effective ) using Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) as shown in Table 5 (Table 6 and Supplementary Figure S4). Thus, the molecular docking, followed by MD simulation and binding free energy estimation using MM-PBSA and MM-GBSA analyses, suggests that Pritelivir can be explored as a promising lead molecule for developing inhibitors against Nsp13 in the therapy COVID-19. It should be noted here that the weak binding of Cordycepin is also because it escapes from the binding pocket near the end of the simulation (Figure 6(c)).

Conclusion
The present study identified two molecules: Cordycepin and pritelivir (BAY 57-1293), targeting the Nsp13 of SARS-CoV-2 and limiting viral replication inside the host. The molecules were selected from an array of 135 nucleoside analogs and were found to have interacting solid efficiency and stability at the inhibition site. H-bond formation of the ligands with the Nsp13 binding site residues indicated their high tendency towards the Nsp13 and was stable throughout the 100 ns MD simulation period. This indicated that our selected  ligands might act in the abnormal activity of NTPase and helicase and consequently, the unwinding of the RNA helices will be hampered, leading to low viral replication and viral load. The physicochemical and ADMET studies evaluated the above-selected inhibitors as highly potent against SARS-CoV-2 infection and are probable drug candidates. The binding free energy calculation using MM-PBSA and MM-GBSA analysis revealed that pritelivir possesses better binding free energy and molecular affinity than cordycepin. Hence, we proposed that pritelivir might be used as lead molecules in COVID-19 therapy. However, in vivo work needs to be carried out to establish the effectiveness of these potential ligands against SARS-CoV-2 and their safety and use.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
The lab is supported by the Science