Identification of Persuasive Antiviral Natural Compounds for COVID-19 by Targeting Endoribonuclease NSP15: A Structural-Bioinformatics Approach

SARS-CoV-2 is a positive-stranded RNA virus that bundles its genomic material as messenger-sense RNA in infectious virions and replicates these genomes through RNA intermediates. Several virus-encoded nonstructural proteins play a key role during the viral life cycle. Endoribonuclease NSP15 is vital for the replication and life cycle of the virus, and is thus considered a compelling druggable target. Here, we performed a combination of multiscoring virtual screening and molecular docking of a library of 1624 natural compounds (Nuclei of Bioassays, Ecophysiology and Biosynthesis of Natural Products (NuBBE) database) on the active sites of NSP15 (PDB:6VWW). After sequential high-throughput screening by LibDock and GOLD, docking optimization by CDOCKER, and final scoring by calculating binding energies, top-ranked compounds NuBBE-1970 and NuBBE-242 were further investigated via an indepth molecular-docking and molecular-dynamics simulation of 60 ns, which revealed that the binding of these two compounds with active site residues of NSP15 was sufficiently strong and stable. The findings strongly suggest that further optimization and clinical investigations of these potent compounds may lead to effective SARS-CoV-2 treatment.


Introduction
Viral diseases are a lethal threat to humans. New viruses emerge, and all major acute respiratory syndromes are affected [1]. SARS-CoV and MERS-CoV, among these, have high fatality rates. SARS-CoV-2 is now undergoing fast dissemination in a pandemic situation. The mortality rate of COVID-19 in the affected population is very high compared to those of SARS-CoV and MERS-CoV; it increased from 2% (January 2020) to more than 7% (May 2020), but is now globally decreasing (at present, 3.28%) [2]. Animal-to-human transmissions were observed in the past, but human-to-human transmissions were also recorded, and they formed the basis for the present pandemic [3].
SARS-CoV-2 is a positive-stranded RNA ((+) RNA) virus, one of a broad family of viruses that includes the Zika, chikungunya, and hepatitis C viruses. (+) RNA viruses bundle their genomes as messenger-sense RNA in infectious virions and only replicate these genomes in replication complexes by RNA intermediates [4]. Upon binding to the human host cell receptor (angiotensin-converting enzyme 2), the virus enters the host cell. The virus releases the RNA into the cytoplasm using host cell systems to begin the translation process; ORF1a and ORFab of the RNA are translated into polyprotein Three H bonds and four hydrophobic interactions were observed in the complex stabilization of NuBBE-242-NSP15. PRO206, GLY165, and ARG207 were found in H-bond interactions, and four different aa, ARG199, LYS205, THR167, and VAL166, were found in hydrophobic interactions.

Pharmacokinetic Properties
To explore the behavior of the selected compounds within an organism, their absorption, distribution, metabolism, and excretion (ADME) were predicted by SwissADME [17] and pkCSM tools. Several therapeutic compounds failed to reach clinical trials because of unfavorable ADME parameters. Lipinski's rule of five [18], Veber's rule [19], Egan's rule [20], polar surface area (TPSA), and a number of rotatable bonds were predicted and are explained in Table 3. The binding energies of complexes NuBBE-1970-NSP15 and NuBBE-242-NSP15 were found to be −483.6 and −305.872 kcal/mol, respectively, which were quite high among the selected compounds. The GOLD fitness scores of these two compounds were found to be 72.2588 and 78.5397, respectively. Besides this, further docking analysis by Autodock [16] was performed to check the inhibition of the selected compounds along with different parameters, namely, their binding energy, inhibition constant, intermolecular energy, and electrostatic energy, as shown in Table 2.

Pharmacokinetic Properties
To explore the behavior of the selected compounds within an organism, their absorption, distribution, metabolism, and excretion (ADME) were predicted by SwissADME [17] and pkCSM tools. Several therapeutic compounds failed to reach clinical trials because of unfavorable ADME parameters. Lipinski's rule of five [18], Veber's rule [19], Egan's rule [20], polar surface area (TPSA), and a number of rotatable bonds were predicted and are explained in Table 3.

Molecular-Dynamics Simulation
MD simulation was performed for 60 ns for both complexes. The constancy of the trajectories was examined via root-mean-square deviation (RMSD). Figure 3 reports the RMSD of the backbone atoms for the receptor and ligand separately for selected complexes NuBBE-1970-NSP15 and NuBBE-242-NSP15. For backbone atoms, it was found that RMSD changes in all systems were initiated at the same point, 0.125 nm. The black line represents the RMSD value change in the backbone of NuBBE-1970-NSP15 that indicates few oscillations until~30 ns; after this time, RMSDs were observed to be steady and remained so until the end of the simulation, with an average value of around 0.37 nm.

Molecular-Dynamics Simulation
MD simulation was performed for 60 ns for both complexes. The constancy of the trajectories was examined via root-mean-square deviation (RMSD). Figure 3 reports the RMSD of the backbone atoms for the receptor and ligand separately for selected complexes NuBBE-1970-NSP15 and NuBBE-242-NSP15. For backbone atoms, it was found that RMSD changes in all systems were initiated at the same point, ~0.125 nm. The black line represents the RMSD value change in the backbone of NuBBE-1970-NSP15 that indicates few oscillations until ~30 ns; after this time, RMSDs were observed to be steady and remained so until the end of the simulation, with an average value of around 0.37 nm. RMSD values in the NuBBE-242-NSP15 backbone ( Figure 3, red plot) showed turbulence from ~26 ns that remained until ~57 ns; however, a little oscillation was observed for ~2 ns at the end of the simulation. Considering the difference between RMSDs of 1 nm from ~40 ns to the end of the simulation, the RMSD of the backbone atoms in both complexes appeared to be stable.
The RMSD plot of the ligands within both complexes ( Figure 3, right plot) showed a sudden elevation at a time within the first 10 ns, and the maximal RMSD value was noted as 0.3 nm. However, RMSD values for the ligand in NuBBE-1970-NSP15 (black) formed a plateau after ~10 ns, and the average value was 0.26 nm; the ligand in NuBBE-242-NSP15 was found to be steady after ~18 ns. Consequently, the RMSDs of the ligands of the two above complexes were found to have constant trajectories during the simulation, unlike those of the backbone atoms, indicating that the ligands showed high stability within the close contact residues of the protein receptor.
For the stability of the complex, RMSD was determined, and found to become steady and remain so until the end of the simulation; the average value was around 0.37 nm. Ignoring residues of the Nand C-terminals, the maximal RMSF value for all simulations reported here (~0.35 nm) was detected in neighboring residues of the active site. For NuBBE-1970-NSP15 and NuBBE-242-NSP15, the range in residues 135-175 was employed in the width of the high RMSF peak. The corresponding pattern was found to be in good agreement with the secondary structure and with the B-factors reported for the initial structure ( Figure 4). RMSD values in the NuBBE-242-NSP15 backbone ( Figure 3, red plot) showed turbulence from 26 ns that remained until~57 ns; however, a little oscillation was observed for~2 ns at the end of the simulation. Considering the difference between RMSDs of 1 nm from~40 ns to the end of the simulation, the RMSD of the backbone atoms in both complexes appeared to be stable.
The RMSD plot of the ligands within both complexes ( Figure 3, right plot) showed a sudden elevation at a time within the first 10 ns, and the maximal RMSD value was noted as 0.3 nm. However, RMSD values for the ligand in NuBBE-1970-NSP15 (black) formed a plateau after~10 ns, and the average value was 0.26 nm; the ligand in NuBBE-242-NSP15 was found to be steady after~18 ns. Consequently, the RMSDs of the ligands of the two above complexes were found to have constant trajectories during the simulation, unlike those of the backbone atoms, indicating that the ligands showed high stability within the close contact residues of the protein receptor.
For the stability of the complex, RMSD was determined, and found to become steady and remain so until the end of the simulation; the average value was around 0.37 nm. Ignoring residues of the Nand C-terminals, the maximal RMSF value for all simulations reported here (~0.35 nm) was detected in neighboring residues of the active site. For NuBBE-1970-NSP15 and NuBBE-242-NSP15, the range in residues 135-175 was employed in the width of the high RMSF peak. The corresponding pattern was found to be in good agreement with the secondary structure and with the B-factors reported for the initial structure ( Figure 4). The minimal distances between residues of the protein receptor and both ligands were calculated, and found to have average values of 0.191 nm for NuBBE-1970-NSP15 and 0.26 nm for NuBBE-242-NSP15 ( Figure 5). Fascinatingly, in the NuBBE-242-NSP15 complex from 38 to 43 ns (for 5 ns), the ligand was moved from the site and found more than 1 nm from the binding site. Although it was for a very short while, the ligand came back to site and remained stable till end of the simulation. This observation suggests that complex-1, i.e., NuBBE-1970-NSP15 was stabler than complex-2, i.e., NuBBE-242-NSP15. Besides this, the robustness of ligands was measured by H-bond approximation at the cavity ( Figure 6) using a cut-off of 0.35 nm. On average, two H bonds were formed between the cavity residues of NuBBE-1970-NSP15 and the atoms of the ligand during the simulation. Two H bonds were found in NuBBE-1970-NSP15 for most of the simulation time, around ~50 ns. Between ~10 and 50 ns of the simulation, two H bonds were considered constant; afterward, fair fluctuation occurred. In NuBBE-242-NSP15, an average of three H bonds was obtained, demonstrating more robust interaction between cavity residues and ligand. During the simulation, the number of bonds exponentially improved after ~30 ns and remained consistent until the end of the simulation time. The minimal distances between residues of the protein receptor and both ligands were calculated, and found to have average values of 0.191 nm for NuBBE-1970-NSP15 and 0.26 nm for NuBBE-242-NSP15 ( Figure 5). Fascinatingly, in the NuBBE-242-NSP15 complex from 38 to 43 ns (for 5 ns), the ligand was moved from the site and found more than 1 nm from the binding site. Although it was for a very short while, the ligand came back to site and remained stable till end of the simulation. This observation suggests that complex-1, i.e., NuBBE-1970-NSP15 was stabler than complex-2, i.e., NuBBE-242-NSP15. The minimal distances between residues of the protein receptor and both ligands were calculated, and found to have average values of 0.191 nm for NuBBE-1970-NSP15 and 0.26 nm for NuBBE-242-NSP15 ( Figure 5). Fascinatingly, in the NuBBE-242-NSP15 complex from 38 to 43 ns (for 5 ns), the ligand was moved from the site and found more than 1 nm from the binding site. Although it was for a very short while, the ligand came back to site and remained stable till end of the simulation. This observation suggests that complex-1, i.e., NuBBE-1970-NSP15 was stabler than complex-2, i.e., NuBBE-242-NSP15. Besides this, the robustness of ligands was measured by H-bond approximation at the cavity ( Figure 6) using a cut-off of 0.35 nm. On average, two H bonds were formed between the cavity residues of NuBBE-1970-NSP15 and the atoms of the ligand during the simulation. Two H bonds were found in NuBBE-1970-NSP15 for most of the simulation time, around ~50 ns. Between ~10 and 50 ns of the simulation, two H bonds were considered constant; afterward, fair fluctuation occurred. In NuBBE-242-NSP15, an average of three H bonds was obtained, demonstrating more robust interaction between cavity residues and ligand. During the simulation, the number of bonds exponentially improved after ~30 ns and remained consistent until the end of the simulation time. Besides this, the robustness of ligands was measured by H-bond approximation at the cavity ( Figure 6) using a cut-off of 0.35 nm. On average, two H bonds were formed between the cavity residues of NuBBE-1970-NSP15 and the atoms of the ligand during the simulation. Two H bonds were found in NuBBE-1970-NSP15 for most of the simulation time, around~50 ns. Between~10 and 50 ns of the simulation, two H bonds were considered constant; afterward, fair fluctuation occurred. In NuBBE-242-NSP15, an average of three H bonds was obtained, demonstrating more robust interaction between cavity residues and ligand. During the simulation, the number of bonds exponentially improved after~30 ns and remained consistent until the end of the simulation time. The trajectories of both systems were subjected to cluster analysis (cut-off value, 0.15 nm) for the temporal distribution to secure convergence. Cluster distribution was categorized by a sequence of flat bars, each continuing for a longer time, suggesting a steady evolution towards a stable relaxed state of the complex. Cluster distribution patterns (Figure 7) were determined by using the trajectories of both systems in fast thermal equilibrium. However, infrequently, the motion of the interior atom determined the creation of one less probable cluster. In Figure 7, a total of six slopes were found, and interestingly, the first slope of structure decline in the NuBBE-1970-NSP15 case was perceived from ~2900 to ~1100 structures; this observation possesses strong agreement with the subplot of cluster distribution with simulation time.
In the end, an average of 10 clusters were found, where ~100 structures of similar confirmations were spotted. In the case of NuBBE-242-NSP15, a total of five slopes were marked. The first slope was observed from ~4100 to ~900 structures, where these structures were reserved in fewer than three clusters. At the end of the simulation, approximately 50 structures were observed, which were held in reserve for ~7 clusters.
The leading clusters in Figure 7 were in rapid equilibrium, with almost full-lifetime structures that lasted for 60 ns; this implied that the system was in fast equilibrium, sampling many conformations that did not differ in their stability. For NuBBE-1970-NSP15, the backbone of the protein receptor was noisily executed, whereas the trajectory of NuBBE-242-NSP15 was less The trajectories of both systems were subjected to cluster analysis (cut-off value, 0.15 nm) for the temporal distribution to secure convergence. Cluster distribution was categorized by a sequence of flat bars, each continuing for a longer time, suggesting a steady evolution towards a stable relaxed state of the complex. Cluster distribution patterns (Figure 7) were determined by using the trajectories of both systems in fast thermal equilibrium. The trajectories of both systems were subjected to cluster analysis (cut-off value, 0.15 nm) for the temporal distribution to secure convergence. Cluster distribution was categorized by a sequence of flat bars, each continuing for a longer time, suggesting a steady evolution towards a stable relaxed state of the complex. Cluster distribution patterns (Figure 7) were determined by using the trajectories of both systems in fast thermal equilibrium. However, infrequently, the motion of the interior atom determined the creation of one less probable cluster. In Figure 7, a total of six slopes were found, and interestingly, the first slope of structure decline in the NuBBE-1970-NSP15 case was perceived from ~2900 to ~1100 structures; this observation possesses strong agreement with the subplot of cluster distribution with simulation time. In the end, an average of 10 clusters were found, where ~100 structures of similar confirmations were spotted. In the case of NuBBE-242-NSP15, a total of five slopes were marked. The first slope was observed from ~4100 to ~900 structures, where these structures were reserved in fewer than three clusters. At the end of the simulation, approximately 50 structures were observed, which were held in reserve for ~7 clusters.
The leading clusters in Figure 7 were in rapid equilibrium, with almost full-lifetime structures that lasted for 60 ns; this implied that the system was in fast equilibrium, sampling many conformations that did not differ in their stability. For NuBBE-1970-NSP15, the backbone of the protein receptor was noisily executed, whereas the trajectory of NuBBE-242-NSP15 was less However, infrequently, the motion of the interior atom determined the creation of one less probable cluster. In Figure 7, a total of six slopes were found, and interestingly, the first slope of structure decline in the NuBBE-1970-NSP15 case was perceived from~2900 to~1100 structures; this observation possesses strong agreement with the subplot of cluster distribution with simulation time. In the end, an average of 10 clusters were found, where~100 structures of similar confirmations were spotted. In the case of NuBBE-242-NSP15, a total of five slopes were marked. The first slope was observed from~4100 to~900 structures, where these structures were reserved in fewer than three clusters. At the end of the simulation, approximately 50 structures were observed, which were held in reserve for 7 clusters.
The leading clusters in Figure 7 were in rapid equilibrium, with almost full-lifetime structures that lasted for 60 ns; this implied that the system was in fast equilibrium, sampling many conformations that did not differ in their stability. For NuBBE-1970-NSP15, the backbone of the protein receptor was noisily executed, whereas the trajectory of NuBBE-242-NSP15 was less disorderly, and it was a plane of continuous change to a new conformation, signified by a series of clusters each having a long lifetime.
Principal-component analysis of the relaxation (calculated for backbone atoms) is presented in Figures 8 and 9. The projection of eigenvectors is depicted in Figure 8. In this presentation, the molecule was allowed to move along the largest vector, corresponding with the dominant conformation transitions during the relaxation of the structure. Figure 9 describes the collection of frames of the protein backbone, where the width of the bands was proportional to the amplitudes the atom makes; narrow bands represent sections that hardly moved, while wide bands are the regions most affected by the transitions. As seen in Figure 9, the red circle, which was supposed to move towards the cavity where ligands had interacted, signifies that, in NuBBE-1970-NSP15, the region occupying the red ellipse showed a transition (wide band) towards the ligand within the cavity. On the other hand, the red-ellipse region demonstrated a significant move towards the cavity in NuBBE-242-NSP15, in strong agreement with the H-bond calculation. However, the black ellipses revealed wide-and narrowband regions for NuBBE-1970-NSP15 and NuBBE-242-NSP15, respectively, suggesting that the width of the band illustrated high amplitude, executing significant functionality.
Molecules 2020, 25, x FOR PEER REVIEW 9 of 17 disorderly, and it was a plane of continuous change to a new conformation, signified by a series of clusters each having a long lifetime. Principal-component analysis of the relaxation (calculated for backbone atoms) is presented in Figures 8 and 9. The projection of eigenvectors is depicted in Figure 8. In this presentation, the molecule was allowed to move along the largest vector, corresponding with the dominant conformation transitions during the relaxation of the structure. Figure 9 describes the collection of frames of the protein backbone, where the width of the bands was proportional to the amplitudes the atom makes; narrow bands represent sections that hardly moved, while wide bands are the regions most affected by the transitions. As seen in Figure 9, the red circle, which was supposed to move towards the cavity where ligands had interacted, signifies that, in NuBBE-1970-NSP15, the region occupying the red ellipse showed a transition (wide band) towards the ligand within the cavity. On the other hand, the red-ellipse region demonstrated a significant move towards the cavity in NuBBE-242-NSP15, in strong agreement with the H-bond calculation. However, the black ellipses revealed wide-and narrowband regions for NuBBE-1970-NSP15 and NuBBE-242-NSP15, respectively, suggesting that the width of the band illustrated high amplitude, executing significant functionality.  In the case of NuBBE-1970-NSP15 backbone transitions, atoms attempted to allocate within 0.4 nm of the ligand-binding site, signifying greater functionality and stability with the ligand. Interestingly, one red and one black ellipse of either complex showed stability and interaction with ligands during the simulation. A video of the principal components of both complexes can be found in the Supplementary Material. Considering the area space between the black and red ellipses in both complexes indicated a noteworthy role in finding close contact residues within 0.4 nm of the ligand atoms. Therefore, NuBBE-1970-NSP15 showed greater robustness and functionality with the ligand, although during the initial stage of simulation, more conformational changes were examined until the ligand in the binding site found stability. In the case of NuBBE-1970-NSP15 backbone transitions, atoms attempted to allocate within 0.4 nm of the ligand-binding site, signifying greater functionality and stability with the ligand. Interestingly, one red and one black ellipse of either complex showed stability and interaction with ligands during the simulation. A video of the principal components of both complexes can be found in the supplementary material. Considering the area space between the black and red ellipses in both complexes indicated a noteworthy role in finding close contact residues within 0.4 nm of the ligand atoms. Therefore, NuBBE-1970-NSP15 showed greater robustness and functionality with the ligand, although during the initial stage of simulation, more conformational changes were examined until the ligand in the binding site found stability.

Discussion
With the rising global number of SARS-CoV-2 infection cases, it is imperative to discover potential therapies in order to cure this devastating disease. Due to their rarity, phytocompounds are important, but due to their weak pharmacokinetics, only few reach the production stage [21]. Recently, chloroquine and hydroxychloroquine were found to have promising in vitro and clinical results [22], but these drugs are highly toxic, and disrupt heart and neurological functions [23,24]. To address the current problems, we utilized a novel approach by collecting natural compounds to check their potential inhibitory action against the critical SARS-CoV-2 target.
The screening was performed using the Discovery Studio 2020 (DS 2020) suite with LibDock as a high-throughput screening tool, followed by docking optimization by CDOCKER and final scoring by calculating binding energies. The prepared library was also screened by GOLD docking software. NuBBE-1970 and NuBBE-242, observed as the common top-scoring compounds (from GOLD and DS screening), were selected and further studied in detail for the inhibition of the NSP15 target protein of COVID-19. The different scores for the ligand, protein, and its complex are shown in Table 4.

Discussion
With the rising global number of SARS-CoV-2 infection cases, it is imperative to discover potential therapies in order to cure this devastating disease. Due to their rarity, phytocompounds are important, but due to their weak pharmacokinetics, only few reach the production stage [21]. Recently, chloroquine and hydroxychloroquine were found to have promising in vitro and clinical results [22], but these drugs are highly toxic, and disrupt heart and neurological functions [23,24]. To address the current problems, we utilized a novel approach by collecting natural compounds to check their potential inhibitory action against the critical SARS-CoV-2 target.
The screening was performed using the Discovery Studio 2020 (DS 2020) suite with LibDock as a high-throughput screening tool, followed by docking optimization by CDOCKER and final scoring by calculating binding energies. The prepared library was also screened by GOLD docking software. NuBBE-1970 and NuBBE-242, observed as the common top-scoring compounds (from GOLD and DS screening), were selected and further studied in detail for the inhibition of the NSP15 target protein of COVID-19. The different scores for the ligand, protein, and its complex are shown in Table 4.
Considering binding-site residues within protein cavities, we recruited residues for further calculation; those were performed using these residues against the corresponding ligand. From a cut-off distance of 0.4 nm within the protein, we recruited close contact residues for NuBBE-1970-NSP15 and NuBBE-242-NSP15; these were GLU203, ASN200, LYS205, SER198, and ILE169. Remarkably, in the case of NuBBE-1970-NSP15, some atoms of the ILE169, THR167, ARG91, and LYS90 residues showed greater fluctuations from the ligand atom; consequently, the minimal distance plot demonstrated moderate inflation at the time of approximately 15 to 42 ns, but later (~17 ns) recovered a suitable distance of less than 0.25 nm. This observation promoted significant additional robustness and an equivalent binding phenomenon with ligands in both cases of the complex. However, ligands recovered the distance up to less than 0.25 nm in the last 20 ns and remained consistent until the end of the simulation.
We concluded that the NuBBE-1970-NSP15 complex was stabler and more robust after binding with NuBBE-1970. Moreover, this work was performed for comparative studies of two potential ligands to find the binding affinity with the protein, and our docking results suggested that NuBBE-1970 binds better. To confirm observation, we performed a molecular-dynamics simulation to check the conformational changes that occur after binding with both complexes. After 15 ns of simulation, a very stable trajectories were achieved, which confirmed that there no more conformational changes occurred during the whole length of simulation.

Data Sources and Preparation
The Nuclei of Bioassays, Ecophysiology and Biosynthesis of Natural Products Database (NuBBE) was accessed, and only phytocompounds (n = 1624) were downloaded. The NuBBE database has proven to be a major platform for research in drug design. The content is publicly accessible and comprises validated multidisciplinary information, chemical descriptors, species origins, geographical locations, spectroscopic data, and pharmacological properties [33]. These compounds were downloaded, imported into Discover Studio 2020 in mol2 format, and processed using the ligand-preparation tool.

Target Structure Preparation
The crystal structure of NSP15 endoribonuclease (PDB ID: 6VWW) was retrieved from the protein data bank (https://supplementary.rcsb.org/structure/6VWW) and resolved at 2.20 Å [34]. The structure was downloaded in PDB format as homo 6-mer. Water molecules were deleted, and the 3D structure was prepared in monomer form using the Discover Studio (DS) 2020 protein-preparation wizard and saved in pdb format for further screening purposes.

Virtual Screening and Docking Experiment
From the multistep route of virtual screening of the prepared NuBBE library, we filtered out the 10 best hits via successive screening by LibDock, docking optimization by CDOCKER, and final scoring by calculating binding energies using inbuilt tools in DS 2020. The GOLD program was also used to screen the prepared NuBBE library to explore the full range of ligand conformational flexibility with partial flexibility of the protein [35,36].
Indepth interaction analysis of the top-ranked molecules with the active site of NSP 15 was performed using Auto Dock Tools [37]. The grid was set using x, y, z of 60, 60, 60, and x, y, z centers −69.377, 26.349, and −34.546, respectively. For docking, the final ligand and protein were prepared using Auto Dock tools. The remaining parameters of the program were kept as their defaults considering a movable ligand and rigid protein. The docking result was visualized using DS 2020.

Molecular-Dynamics Simulation
GROMACS 4.6.7 packages [38,39] were used for preparing the system and performing MD simulations using the gromos53a6 force field [40]. The protein solute was solvated by explicit SPC216 water [41] in a dodecahedron box with a margin of 10 Å between solute and box walls. Systems were brought to neutrality by the addition of sodium counter ions. A 1 nm cut-off distance was taken under the particle-mesh Ewald method [42] to calculate long-range electrostatic interactions, and a 1 nm cut-off distance was considered to evaluate van der Waals interactions. The LINCS algorithm of fourth-order expansion was used to constrain bond lengths [43]. The steepest-descent algorithm was applied to optimize for 10,000 steps to remove steric clashes between atoms. The system was equilibrated for 1 ns with position restraints of all heavy atoms. Berendsen weak-coupling schemes were used to maintain the system at 300 K and 1 atom using separate baths for the system. Initial velocities were randomly generated using Maxwell-Boltzmann distribution corresponding to 300 K. Lastly, the production run was performed for 60 ns. Furthermore, xmgrace (http://plasma-gate.weizmann.ac.il) was used for preparing graphs. Ligand-topology preparation was implemented by using the PRODRG server with the option of choosing no chirality, full charge, and no energy minimization [44]. The equilibrated structure obtained after equilibration was considered as the reference structure, and trajectories were fitted to the backbone of this structure. VMD was used to load and visualize all frames of the principal components.

Conclusions
Due to its growing infection and mortality rates, with no approved drugs or vaccines available, the COVID-19 pandemic is a serious public-health threat. The purpose of this research was to computationally explore natural substances that could potentially inhibit SARS-CoV-2 NSP15, a well-known protein for viral replication in host cells. After multiscoring virtual-screening cascade using DS and GOLD tools, followed by molecular-docking and MD simulations of the Natural Compounds Library (NuBBE database) against NSP15, we described the top-ranked compounds, NuBBE-1970 and NuBBE-242, showing the strongest binding and stable complexation with NSP15. Eleven aa of the NSP15 active site interacted with NuBBE-1970, while 13 aa were found to interact with NuBBE-242. Molecular-dynamics analysis showed that NuBBE-1970 was stabler with NSP-15 than NuBBE-242 was. These compounds were found to be nontoxic and to satisfy drug-likeness properties in our in silico analysis. On the basis of our results, we conclude that these natural compounds should be considered for further studies in the search for COVID-19 therapies.

Conflicts of Interest:
The authors declare no conflict of interest.