High throughput virtual screening reveals SARS-CoV-2 multi-target binding natural compounds to lead instant therapy for COVID-19 treatment.

The present-day world is severely suffering from the recently emerged SARS-CoV-2. The lack of prescribed drugs for the deadly virus has stressed the likely need to identify novel inhibitors to alleviate and stop the pandemic. In the present high throughput virtual screening study, we used in silico techniques like receptor-ligand docking, Molecular dynamic (MD), and ADME properties to screen natural compounds. It has been documented that many natural compounds display antiviral activities, including anti-SARS-CoV effect. The present study deals with compounds of Natural Product Activity and Species Source (NPASS) database with known biological activity that probably impedes the activity of six essential enzymes of the virus. Promising drug-like compounds were identified, demonstrating better docking score and binding energy for each druggable targets. After an extensive screening analysis, three novel multi-target natural compounds were predicted to subdue the activity of three/more major drug targets simultaneously. Concerning the utility of natural compounds in the formulation of many therapies, we propose these compounds as excellent lead candidates for the development of therapeutic drugs against SARS-CoV-2.

Introduction

<span class="Species">Coronaviruses (<span class="Species">CoV) have a wide range of hosts, including several species of animals and birds [1]. They belong to a group of enveloped viruses that have positive-sense RNA as genetic material [2]. In 1965, Tyrrell and Bynoe isolated the first Human coronaviruses (HCoV) from the nasal washings of a male child suffering from symptoms of a common cold [3]. HCoV came into the spotlight during the outbreak of Severe Acute Respiratory Syndrome (SARS) in China during the fall of 2002 [[4], [5], [6]]. The SARS-CoV appeared to have originated from the bat and civet reservoir [7,8], claimed 8096 probable cases and 774 deaths in >30 countries [9,10]. Ten years after of SARS episode, in 2012, there was another outbreak of SARS-like disease caused by a different verity of HCoV, officially named as Middle-East Respiratory Syndrome Coronavirus (MERS-CoV) [11]. In 2012, MERS-CoV claimed 2279 laboratory-confirmed cases in 27 different countries [12]. Genome sequencing analysis has revealed MERS-CoV was also of bat origin [13,14]. In Dec 2019, complicated cases of pneumonia were reported in Wuhan, China [15]. World Health Organization (WHO) announced the name of the causative virus as 2019 novel coronavirus (2019-nCoV) on January 12, 2020. Later, the Coronavirus Study Group (CSG) of the International Committee on Taxonomy of Viruses named 2019-nCoV as Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) [16]. As of March 24, 2020, SARS-CoV-2 has cumulatively infected over 381,653 people worldwide and killed over 16,558 individuals since it is first detected. Recent evidence reveals that SARS-CoV-2 was a chimeric virus of a bat coronavirus and an unknown coronavirus [17]. The current epidemic and faster transmission to almost 190 countries and territories are making the situation alarming worldwide. The threat of the currently emerging SARS-CoV-2 necessitates therapeutic strategies targeting proteins involved in the replication and transcription processes of the viral RNA genome.

Since the SARS outbreak of 2002, a tremendous amount of basic and clinical research on <span class="Species">coronaviruses has been performed, which led to the identification of many potential drug targets [[18], [19], [20]]. Out of the proteases involved in polyprotein proteolysis, <span class="Gene">3C-like protease (3CLpro) (EC 3.4.22.69) has been well studied [21]. Antiviral activities of a series of 3CLpro inhibitors have been reported for Feline Coronaviruses and Feline Caliciviruses [22]. Furthermore, recently, the X-ray crystal structure of 3CLpro of SARS-CoV-2 and its complex with an α-ketoamide inhibitor have been reported, one of the α-ketoamide inhibitors has shown pronounced inhibition of SARS-CoV replication in human Calu-3 lung cells [23].

<span class="Species">Coronavirus <span class="Gene">helicase catalyses the processive separation of double-stranded DNA and RNA in a 5′-to-3′ direction [24,25]. Moreover, helicase (EC 3.6.4.12) inhibitors such as bananin and its derivatives have shown promising anti-cornonaviridae activity [26]. The non-structural protein 15 have endoribonuclease (EC 3.1.26.3) activity that cleaves preferentially at uridine residues [27,28]. A non-structural protein (nsp) 14 has also been shown to have exoribonuclease (ExoN) (EC 3.1.13.1) activity and involved in CoV replicative machinery. The active site mutations of the protein revealed severe defects in CoV RNA synthesis [29]. In the SARS-CoV replication study, it has been shown that AdoMet dependent methyltransferase inhibitors inhibit SARS-CoV replication [30] by targeting N-methyltransferases (EC 2.1.1.35) of the virus [31]. Due to strictly virus-specific activity, the RNA-dependent RNA polymerase (RdRp) (EC 2.7.7.6) has been considered as a promising wide-spectrum drug target for antiviral drug development. Several studies discovered that inhibitors against RdRp could effectively intervene in the coronavirus lifecycle [32,33]. RNA viruses are unique pathogens in terms of exceptionally genetically variable. They accumulate mutations in their genome by employing viral reverse transcriptase that lacks a proofreading mechanism. This unique feature of RNA viruses makes it challenging to design active therapeutic agents. Most of the antiviral drugs in use are designed to specifically target single viral enzyme, which is essential for viral replication or invasion [34,35]. However, the high rate of mutations in the viral drug targets has been accounted for the reduced susceptibility of currently available antiviral drugs [36]. A combination of drugs that have different molecular targets can be a better choice; however, some times, the combination therapy is not safe due to unwanted drug interactions [37]. Besides, drugs designed for multiple protein targets are extensively used for the treatment of both infectious and inherited diseases [[38], [39], [40], [41], [42]]. For the immediate drug requirement for SARS-CoV-2, biologically active drug-like molecules that target multiple viral enzymes of the viral replication cycle are highly wanted.

Natural compounds of plant-based origin have been studied as an exciting class of pharmacologically active molecules, some of them have an ancient history of antiviral activity [[43], [44], [45], [46], [47], [48], [49], [50]]. Extensive studies have been performed over the past decades to identify anti-<span class="Species">CoV agents using natural products [51]. The saikosaponins <span class="Gene">(A, B2, C, and D), for example, exerts antiviral activity by interfering viral attachment and penetration, against HCoV-22E9 [52]. Besides, several natural products such as Lycoris radiata, Artemisia annua, Pyrrosia lingua, and Lindera aggregate have been reported to display significant anti–SARS-CoV properties [53]. Moreover, inhibitors from natural origin have been identified against the SARS-CoV enzymes, such as helicase and 3CLpro and viral RdRp [[54], [55], [56], [57]]. NPASS database is freely accessible (http://bidd2.nus.edu.sg/NPASS/) that provides the literature-reported experimentally-determined activity (e.g., IC50, Ki, EC50, GI50, and MIC) values of the natural products against macromolecule or cell targets along with the taxonomy of the species sources of 35,032 unique natural products [58]. In the heart of the current Corona Virus Disease 2019 (COVID-19) outbreak, these NPASS compounds may be used for capable therapy as they can remarkably reduce the time taken to design a therapeutic regimen.

Structure-based drug design by virtual screening and molecular docking studies has become a valuable primary step in the identification of novel lead molecules for the treatment of diseases [59,60], and proven to be a very efficient tool for antiviral [[61], [62], [63], [64]] and antibacterial [65,66] and antiprotozoal [67,68] drug dis<span class="Species">covery. Therefore, a virtual screening experiment was conducted to determine the interaction of natural ligands of the NPASS database within the binding pocket of putative drug targets of the virus that was calculated in terms of docking scores and MM-GBSA values. Our high throughput virtual screening revealed 21 natural compounds having higher docking scores and MM-GBSA values for six potent therapeutic targets of <span class="Species">SARS-CoV-2 over the known chemical inhibitors. Remarkably, we suggested three natural compounds that able to bind the catalytic site of three/more crucial viral enzymes with a relatively high affinity, which ultimately can be used for the development of instant drugs for the emerging COVID-19.

Material and methods

Selection of different drug targets of SARS-CoV-2 and its sequence similarity with SARS coronavirus

For developing the structure of <span class="Species">SARS-CoV-2 functional enzy<span class="Chemical">mes, the amino acid sequences of SARS-CoV-2 (accession NC_045512.1) were downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/) in the FASTA format. All of the six proteins namely helicase (accession YP_009725308.1), endoribonuclease (accession YP_009725310.1), exoribonuclease (accession YP_009725309.1), RNA dependent RNA polymerase (RdRp) (accession YP_009725307.1), methyltransferase (accession YP_009725311.1) and 3C-like proteinase (accession YP_009725301.1) belong to the replication complex of the deadly virus SARS-CoV-2. The amino acid sequences obtained from NCBI were aligned with SARS coronavirus using the BLASTp server (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins) [69].

Homology modeling of the selected drug targets, refinement, and validation of structure

Since the crystal structures of the selected drug targets were not available on the protein data bank (<span class="Gene">PDB), the 3D structures were modeled using SWISS-MODEL (https://swissmodel.expasy.org/). For this purpose, the amino acid sequences of each target were submitted in the SWISS-MODEL server to develop the tertiary structure [70]. Here, we had selected the template similar to the query sequence <span class="Species">coverage and identity of the sequence-based upon their Global Model Quality Estimate (GMQE) [71] and Quaternary Structure Quality Estimate (QSQE) scores. The homology modeling technique, which we use to predict the tertiary structure of the protein, is the widely used method. However, accurate estimation of the three-dimensional position of individual atoms in a protein sequence is tough even with the best-matched template and target sequence alignment [[72], [73], [74], [75], [76]]. The homology model generally deviates from their native structure concerning their atomic coordinates [77]. Therefore, the refinement of the homology model is a very crucial step to identify the accurate near-native structure [78]. It is known that the geometrical/predicted structure of the target sequence affects the function of the protein, which also includes pharmacophore drug designing.

Here, we have used the 3D-refine server (http://sysbio.rnet.missouri.edu/3Drefine/) for the refinement of the modeled structures of each target protein of <span class="Species">SARS-CoV-2. This refinement server works on the two-step protocol, which reliably brings the predicted homology model closer to its native structure [[79], [80], [81], [82], [83], [84]]. Where the first step is the optimization of the <span class="Chemical">hydrogen bond network, and second is the minimization of atomic-level energy of optimized homology models using the knowledge-based force fields [85]. This server requires a homology model in PDB format as an input query. The server provides the five refined models as an output with best on the top in PDB format. There are several parameters for the selection of best-refined models, which include the Molprobity score, which is the measure of the local quality of the structure [86], the lower score shows the excellent model. Another given parameters as output scores are GDT-TS score, GDT-HA score, and RMSD, which tells about the positioning of Calpha atoms. GDT-TS and GDT-HA score of the refined model, which ranges from 0 to 1, where the lower score shows a good quality model [87]. The .pdb files of refined structures of respective target proteins were uploaded on the RAMPAGE server (http://mordred.bioc.cam.ac.uk/~rapper/rampage.php), which had generated the Ramachandran plot. This plot is generated to obtain the value of the total count of amino acids present in the allowed, disallowed, and favorable regions [88].

Preparation of protein and ligands using Schrödinger suite

Refined homology models of the <span class="Disease">COVID-19 drug targets (proteins) were used for the further downstream steps. These protein structures were defined as a receptor and opened via protein preparation wizard in maestro v11.9 [89]. The standard refined model of protein does not fit for the docking with ligand molecules vis-a-vis calculation of the binding affinity. Therefore, models were preprocessed and optimized by filling the missing side chain using prepwizard module of the Schrodinger suite 2019-1. Heavy <span class="Chemical">water molecules were removed with less than 3H-bonds followed by energy minimization, which convert the heavy atoms to RMSD 0.30 Å.

Natural compounds from the NPASS database [58] were used as ligands to target the <span class="Disease">COVID-19 replication complex proteins. The compounds having druggable properties were identified (3963 natural compounds), and 2D structures in .sdf format were downloaded from the NPASS database. For using these compounds as ligands for docking purposes, there is a need to remove the <span class="Chemical">salt, addition of hydrogen atoms, and deprotonation. LigPrep v4.7 module of the Schrodinger suite was used to prepare the ligands in the three-dimensional structure. Charge neutralization to attain the biological relevant pH (pH 7), desalting was done using Epik v4.7, which is based upon Hammett and Taft methodologies. Torsional bond of ligands was released, and a maximum one stereoisomer per ligand was generated, which resulted in the final 4570 ligands upon LigPrep of natural compounds from the NPASS database.

In this study, one reference/control ligand for each target protein was included. Those ligands were chosen based upon the previous studies. Out of six, four inhibitors namely MLS001181552 (inhibitor of <span class="Gene">helicase) [90], <span class="Chemical">remdesivir (inhibitor ofRdRp) [91], sinefungin (inhibitor of methyltransferase) & TG-0205221 (3C-like proteinase) [92] were well-known for SARS coronavirus protein targets while rest two namely RO-7 (inhibitor of endoribonuclease) [93] and 2-morpholine-4-ylethanesulfonate (inhibitor of exonuclease) [94] inhibitors were known for Influenza virus protein and member of DEddH family of exonucleases (Fig. 1 ).

Fig. 1

Showing the chemical formula of reference/control ligand for each target protein.

Generation of grid

The generation of the grid is a very crucial step for the binding of a ligand to the receptor. Here, a 3-dimensional boundary for the ligand binding was generated in the receptor using Glide, version 8.2 [95] of Maestro, Schrodinger. The receptor grid for each target protein was generated by indicating the active sites amino acid residues, which were searched from previously reported studies appropriate to each target. Amino acid residues involved as the active site with respect to the target, along with their references, are available in Supplementary Table 1. The size of the receptor grid was set at default, i.e., 20 Å. Active site amino acids are crucial because they may affect the entry of the ligand, followed by binding to the target protein.

Molecular interaction/docking studies of six target proteins of COVID-19 against natural compounds

After grid generation, the natural compounds/ligands were docked to the target protein/receptor using the protocol of Glide grid version v8.2 [95]. The software internally generates different conformations, which passes across several filters viz. euler angles, grid-based force field evaluation, and Monte Carlo energy minimization. Lastly, the evaluation of conformers takes place based upon the docking score, and one of the best conformations per ligand is generated as an output. The docking was performed in three steps, which include high throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP). The first two steps of docking, HTVS, and SP utilize the self-same scoring function, whereas XP reduces the intermediate conformations and thoroughness of the torsional refinement and sampling [96]. A total of 4570 natural compounds were docked against each target protein with their respective positive control inhibitor (mentioned earlier). Subsequently, top scorers were set forth for SP docking, and the output of SP docking was put forward in XP docking with positive control. The docking score is a binding energy/affinity given in the kcal/mol. All the statistical data of top scorers for each target protein wrt control obtained is presented in result tables.

Absorption, Distribution, Metabolism, Excretion (ADME) properties evaluation/calculation

It has been shown that the natural compounds do not follow Lipinski's rule as these compounds tend to keep low hydrophobicity as well as the potential of donating the intermolecular H-bond [97]. We have used QikPropv5.9 of the Schrodinger suite, which analyzes the druggable properties of the ligands. This module takes ligprep file as an input and predicts the ADME properties in the form of several physicochemical properties and principal descriptor values, which were found in the acceptable range.

Generation of molecular mechanics, the generalized born solvent accessibility (MM-GBSA) score and ROC plot analysis

Another application of Schrodinger suite software viz., Prime, was used to estimate the binding energy and represented as MM-GBSA score of the selected docked complex of the natural compound and target protein obtained from extra precision (XP) step of docking. This post-docking application provides the correct ranking of ligands binding based upon the MM-GBSA energy score [98].

Further, there are several other metrics, which are currently in use to evaluate the performance of virtual screening docking protocol, for instance, enrichment factor (EF), area under the ROC curve, and Boltzmann-enhanced discrimination of receiver operating characteristic (BEDROC) etc. [99]. In this study, receiver operating characteristic (ROC) plots were also generated to assess the performance of the SP and XP docking methodology, which differentiates the active compounds with false positives means decoys [100]. These statistical curves represented the probability and were plotted to correctly classify the ligands into actives and decoys that finally attested our docking performance. It gives the value between 0 and 1, where 0 represents the worst performance of docking, and 1 represents the best performance [99,101]. To plot these curves, sensitivity (on Y-axis) was plotted against 1-specificity (on X-axis). Here, sensitivity signifies the true positive, whereas 1-specificity signifies the decoys/false positive. All-natural compounds/ligands performed best on SP and XP docking were subjected to the docking protocol validation. Separate ROC curves to validate for both XP and SP docking protocols were generated. In addition to this, we have also calculated the Enrichment factor (an enrichment calculator from Schrodinger suite was used), Area Under Accumulation (AUC), and BEDROC values were also predicted.

Molecular dynamics

For comprehending the properties of the structure and their microscopic interaction, molecular dynamics (MD) simulation was performed. For this, systems were built for best three multi-target natural compounds with one of their potent targets viz. <span class="Gene">helicase and <span class="Gene">NPC270578, exoribonuclease, and NPC214620 & methyltransferase and NPC52382 complex using Desmond v5.6, system builder panel (Schrodinger 2019-1). The optimized potentials were used for the liquid simulations at the OPLS3e force field, and these solvated systems were opened in the molecular dynamics panel [102]. Firstly, the solvated system builder was used to selecting a single point charge as a water model, and ligand-protein complexes were neutralized by adding sodium or chloride ions. The NPC214620-exonuclease complex model system consisted of 80,840 atoms and 24,147 water molecules and was neutralized by adding 7 Cl− ions. NPC52382-methyltransferase model system contained 34,012 atoms and 9781water molecules and was neutralized by adding 2 Cl− ions. Whereas, the NPC270578-helicase complex model system consisted of 79,549 atoms and 23,354 water molecules and was neutralized by adding 9Cl− ions. System equilibration was done at a constant temperature, i.e., 300.0 K and standard atmospheric pressure, i.e., 1 bar. Upon the equilibration step, the total time for each simulation was appointed ten nanoseconds with trajectories analysis at regular intervals for the protein-ligand complex. Ultimately, the RMSD was analyzed to identify the stability of the target protein in its natural motion.

Results

Developing the structure of SARS-CoV-2 functional enzymes

The amino acid lengths of the selected protein targets for <span class="Species">SARS-CoV-2 were found to be 601 for <span class="Gene">helicase, 346 for endoribonuclease, 527 for exoribonuclease, 932 for RNA-dependent RNA polymerase, 298 for methyltransferase and 306 for 3C-like proteinase. It was observed that these target proteins from SARS coronavirus were also carrying almost similar amino acid length, i.e., 603, 346, 527, 955, 344, and 308, respectively. In this piece of study, we wanted to investigate the identity, similarity, and percentage of conserved sequences between COVID-19 and SARS coronavirus.

Models for each enzyme such as <span class="Gene">helicase, endoribonuclease, <span class="Gene">exoribonuclease, RNA dependent RNA polymerase, methyltransferase, and 3-C like proteinase were generated, using template 6jyt.2.A, 2h85.1.A, 5nfy.1.A, 6nur.1.A, 2xyr.1.A and 2z9j.1, respectively. The model quality of the generated 3D models was also estimated by the QMEAN score [103] function of SWISS-MODEL, which relies on modeling error quantification and expected model accuracy. The GMQE and QSQE scores generally lie between 0 and 1, where higher the score, the higher the reliability of the modeled structure means expected accuracy of interchain interactions for a given template, whereas, in case of QMEAN, a score lower than 4.0 gives reliability [104]. After the model generation and refinement, each enzyme's structure was validated through the generation of the Ramachandran plot (Table 1 ). Ramachandran plot confirmed that all target models after refinement would be used for molecular docking with ligand except 3C-like proteinase. In case of 3C-like proteinase the model obatined from Swiss-modeling was used for molecular docking. The structure generation confirms the amino acid mapping in the allowed region, which shows the stable nature of existing proteins.

Table 1

Details of scores obtained for target proteins via RAMPAGE and Swiss-model servers.

Proteins	Template used for homology modeling	Sequence identity (%) by Blast (protein seq. with template seq.)	Ramachandran score (% of residue)	GMQE score	RMSD (Å)	Ramachandran score after refinement (% of residue)
Helicase	6jyt.2. A	97.95	84.34	0.98	0.304	88.6
Endoribonuclease	2h85.1. A	88.0	97.96	0.98	0.249	98.5
Exoribonuclease	5nfy.1. A	94.07	89.44	0.98	0.258	93.1
RNA-dependent RNA polymerase	6nur.1. A	96.24	97.5	0.83	0.228	98.1
Methyltransferase	2xyr.1. A	92.57	96.18	0.97	0.252	96.2
3C-like proteinase	2z9j.1.A	96.0	96.51	0.99	0.269	96.4

Details of scores obtained for target proteins via <span class="Gene">RAMPAGE and Swiss-model servers.

Computational Virtual Screening of Natural Ligands against potential therapeutic targets of SARS-CoV-2

To screen effective inhibitors for <span class="Species">SARS-CoV-2, we used compounds of NPASS database against six potential therapeutic targets of the virus, such as <span class="Gene">helicase, endoribonuclease, exoribonuclease, RNA-dependent RNA polymerase, N-methyltransferase, and 3C-like protease in our molecular docking experiment. HTVS, SP, and XP step-wise screening protocol were followed to find potent compounds having higher docking scores and binding energy as follows.

Molecular docking with helicase

In a <span class="Species">SARS-CoV replicon assay, it has been shown that the specific inhibition of <span class="Species">SARS-CoV helicase by SSYA10-001 blocked viral replication [90]. Therefore, we used this compound as a control to find out more effective inhibitors for COVID-19 during our high throughput virtual screening experiment of natural compounds. Interestingly, we found five compounds such as NPC270578, NPC52382, NPC473043, NPC175107, and NPC22192 displayed docking score and MM-GBSA value 2.5 and 1.5 times greater than the control, respectively (Table 2 ). Notably, NPC52382 has been shown to have biological activity against the malaria parasite in the NPASS database.

Table 2

Docking and MMGBSA score of the top targeted compounds.

Molecular docking with endoribonuclease

The antiviral property of endoribonuclease inhibitors has been validated for the <span class="Species">influenza virus in multiple in vitro assay [93]. The essential function of endonuclease in the initiation of viral transcription supports its potential as a promising target for the development of antiviral agents. We find five compounds such as NPC169474, NPC297657, NPC19721, NPC279121, and NPC10737, which have more than two-fold docking scores and better MM-GBSA value than the known viral endoribonuclease inhibitor (Table 2). Subsequent biological activity analysis of these compounds in the NPASS database reveals that NPC279121 has antiviral activity against <span class="Species">Hepatitis C virus by inhibition of RdRp [105], Influenza A virus by inhibition of neuraminidase [106], HIV-1 by inhibition of integrase [107] and Simian virus 40 by unknown mechanism whereas NPC10737 has antiviral property (IC50 = 34 μg/ml) against HIV1 through suppression of transcription from 5′-long terminal repeat including activation via NF-kB [108].

Molecular docking with exoribonuclease

Inhibitors for exonuclease of the Lassa <span class="Disease">fever virus have been identified [94]. Our molecular docking studies identified five compounds, such as NPC137813, NPC191146, NPC3825, <span class="Gene">NPC270578, and NPC52382, bearing 1.8 to 2.7fold higher docking score and 2 to 2.5fold higher MM-GBSA value as compared to the general lead exonuclease inhibitor MES (2-(N-morpholino) ethanesulfonic acid) (Table 2). Surprisingly, we did not get any hints of the antiviral studies of these compounds in the NPASS database.

Molecular docking with RdRP

It is essential to find out novel compounds against <span class="Disease">COVID-19 based on broad-spectrum antiviral compounds. Several studies revealed the potent antiviral activities of <span class="Chemical">remdesivir through the inhibition of RdRp [[109], [110], [111]]. Further, remdesivir has also shown antiviral activities against SARS-CoV and MERS-CoV [111]. To study whether the NPASS database has some compounds active against RdRP of SARS-CoV-2, we conducted molecular docking studies. We found a compound, NPC161224, having a 1.35 fold higher docking score but with a lower MM-GBSA value than remdesivir (Table 2). According to the NPASS database, this compound has also shown RdRP inhibition property against HIV1.

Molecular Docking with 3C-like proteinase (3CLpro)

<span class="Gene">3CLpro is an attractive target for antiviral therapeutics because of its essential role in the processing of <span class="Species">CoV polyprotein [112,113]. Numerous 3CLpro inhibitors have been shown to block CoV replication in cell culture [114,115]. TG-0205221, a known inhibitor of 3CLpro of SARS coronavirus, could significantly reduce the viral titer of SARS CoV in cell culture [92]. The binding mode information of TG-0205221 is employed here towards the discovery of new 3CLpro inhibitors from the NPASS database. From our intensive docking studies, we report five molecules such as NPC19709, NPC61506, NPC107109, NPC130230, and NPC175552 having three-fold higher docking score than TG-0205221. However, these molecules have a relatively lower MM-GBSA value than the control (Table 2). Notable, the antiviral activity (IC50 = 14.48 μg/ml) against Influenza A virus of NPC130230 by inhibition of neuraminidase [116] has been described in the NPASS database.

Molecular docking with methyltransferase

Methylation of the viral mRNA cap structure at the N7 position of the <span class="Chemical">guanine is indispensable for the synthesis of viral proteins [117]. N7- methyltransferase has been shown as a potential antiviral target for <span class="Species">SARS-CoV in a yeast-based screening assay [118]. In our docking experiment, three natural compounds such as NPC226294, NPC270578, and NPC52382 has shown better docking scores and MM-GBSA value than, Sinefungin, the known inhibitor of N7- methyltransferase (Table 2). Besides, we get one compound NPC70622 that has antiviral activities such as IC50 = 300 μg/ml against the Hepatitis B virus by inhibition of HBV-DNA production [119] and EC50 = 0.2 μM against Hepatitis C virus by inhibition of RdRp [105] as described in the NPASS database.

Identification of potential multi-target inhibitors of SARS-CoV-2

Taking the docking scores and MM-GBSA values displayed by natural compounds for three/more targets into consideration, we propose three compounds as multi-target inhibitors against SARS-CoV-2. Compound NPC214620 secures better docking scores than respn>ective control molecules against five different drug targets. The docking scores and MM-GBSA values of the molecule in the binding pocket of respn>ective targets are shown in Fig. 2A. As compared to controls, the second compound NPC52382 showed higher docking scores and MM-GBSA values for four therapeutic targets, as shown in Fig. 2B. The third multi-target compound that exhibits higher docking scores and MM-GBSA values towards helicase, endoribonuclease, and N-methyltransferase is NPC270578, shown in Fig. 2C.

Fig. 2

(A). Showing the contacts between multi-target NPASS compound, NPC214620 with (1) Endoribonuclease (2) Exoribonuclease (3) RNA dependent RNA polymerase (RdRp) (4) Methyltransferase and (5) 3C-like proteinase. The ligand 2D representation in the side panel is showing the rotatable bonds of the ligand in different colors. (B). Showing the contacts between multi-target NPASS compound, NPC52382 with (1) Helicase (2) Endoribonuclease (3) Exoribonuclease (4) Methyltransferase. The ligand 2D representation in the side panel is showing the rotatable bonds of the ligand in different colors. (C). Showing the contacts between multi-target NPASS compound, NPC270578 with (1) Helicase (2) Exoribonuclease (3) Methyltransferase. The ligand 2D representation in the side panel is showing the rotatable bonds of the ligand in different colors.

(A). Showing the contacts between multi-target NPASS compound, NPC214620 with (1) Endoribonuclease (2) Exoribonuclease (3) RNA dependent RNA polymerase (<span class="Gene">RdRp) (4) Methyltransferase and (5) <span class="Gene">3C-like proteinase. The ligand 2D representation in the side panel is showing the rotatable bonds of the ligand in different colors. (B). Showing the contacts between multi-target NPASS compound, NPC52382 with (1) Helicase (2) Endoribonuclease (3) Exoribonuclease (4) Methyltransferase. The ligand 2D representation in the side panel is showing the rotatable bonds of the ligand in different colors. (C). Showing the contacts between multi-target NPASS compound, NPC270578 with (1) Helicase (2) Exoribonuclease (3) Methyltransferase. The ligand 2D representation in the side panel is showing the rotatable bonds of the ligand in different colors.

QikProp analysis of ADME and principal descriptor values

Lipinski's Rule of Five is not a strict criterion for natural compounds [120]. Therefore, some crucial ADME properties for the lead compounds were evaluated using Schrӧdinger QikProp module. Various basic physiochemical properties such as PlogPo/w (Predicted octanol/<span class="Chemical">water partition coefficient), PlogS (Predicted aqueous solubility), PlogHERG (Predicted IC50 value for the blockage of HERG K+ channels), PPCaco (Predicted apparent <span class="CellLine">Caco-2 cell permeability for the gut blood barrier), PlogBB (predicted brain/blood partition coefficient), PPMDCK (Predicted apparent MDCK cell permeability in nm/s), and PlogKhsa (prediction of binding to human serum albumin) of these compounds were predicted. The values of those compounds found in the acceptable range were shown in Table 3 . The principal descriptor values including MW (molecular weight of molecule), SASA(Total solvent accessible surface area), FOSA (Hydrophobic component of the SASA), FISA (Hydrophilic component of the SASA), PISA (π (carbon and attached hydrogen) component of the SASA), PHOA (Percent human oral absorption), DonorHB (Estimated number of hydrogen bonds that would be donated by the solute to water molecules in an aqueous solution), and Accept HB (Estimated number of hydrogen bonds that would be accepted by the solute from water molecules in an aqueous solution) were evaluated for the drug-like behavior of the lead compounds. The values of the compounds fallen in the suitable range were displayed in Table 4 .

Table 3

ADME Properties of the selected NPASS compound.

NPASS compound ID	QPlogPo/w (−2.5 to 6.5)	QPlogS (−6.5 to 0.5)	QPlogHERG (above −5.0)	QPPCaco (<25%-poor, >500-great)	QPlogBB (−3 to 1.2)	QPPMDCK (<25%-poor, >500-great)	QPlogKp (−8.0 to −0.1)	QPlogKhsa (−1.5 to 1.5)
NPC214620	−1.19	−0.09	−2.24	180.67	−0.96	77.83	−4.61	−0.91
NPC52382	−0.78	−1.97	−3.55	384.25	−0.91	175.94	−3.97	−0.88
NPC270578	−0.78	−1.97	−3.55	384.25	−0.91	175.94	−3.97	−0.88
NPC473043	−1.02	−1.92	−3.24	113.29	−1.32	46.99	−5.01	−0.78
NPC175107	−2.75	−1.41	−3.34	65.68	−1.56	26.07	−5.47	−1.29
NPC22192	1.97	−5.89	−5.45	246.55	−1.52	108.91	−4.45	0.10
NPC169474	0.52	−3.52	−4.92	870.98	−0.68	426.10	−3.48	−0.56
NPC297657	0.38	−2.61	−3.92	131.17	−1.50	55.06	−4.79	−0.43
NPC19721	0.02	−1.53	−3.38	54.84	−1.76	21.45	−5.43	−0.61
NPC279121	0.55	−1.89	−2.97	100.45	−1.31	41.27	−5.11	−0.33
NPC10737	0.42	−3.37	−4.68	516.13	−0.95	242.04	−3.82	77.97
NPC137813	−1.52	−2.57	−5.05	713.22	−0.74	352.70	−3.72	−1.48
NPC191146	3.08	−4.94	−4.32	1406.61	−0.43	715.33	−2.98	0.36
NPC3825	0.98	−3.29	−4.12	316.75	−1.05	142.79	−4.23	−0.32
NPC161224	−8.84	2.00	−4.83	5.58	−1.58	4.22	−8.63	−3.37
NPC226294	−1.15	−1.88	−3.67	178.04	−1.31	76.60	−4.53	−0.95
NPC19709	−0.891	−1.973	−3.159	259.635	−0.931	115.172	−4.401	−0.766
NPC61506	1.18	−2.972	−4.431	85.894	−1.889	34.841	−5.047	−0.226
NPC107109	0.436	−0.906	−2.725	607.8894	−0.55	288.864	−3.683	−0.769
NPC130230	0.50	−1.06	−2.60	154.93	−1.13	65.92	−4.65	−0.58
NPC175552	0.08	−1.76	−3.51	203.24	−1.22	88.39	−4.42	−0.60

Table 4

Principal drug characteristic predictions for NPASS compounds.

NPASS compound ID	MW (130–725)	SASA (300−1000)	FOSA (0–750)	FISA (7–330)	PISA (0–450)	Volume (500–2000)	PHOA (<25%-poor, >80%-high)
NPC214620	170.12	273.73	90.35	183.38	0.00	398.07	60.39
NPC52382	450.40	492.93	344.11	148.82	0.00	842.51	42.72
NPC270578	405.40	492.93	344.11	148.82	0.00	842.51	42.72
NPC473043	448.38	468.84	264.08	204.76	0.00	807.40	31.81
NPC175107	492.39	505.02	275.29	229.29	0.00	893.00	17.48
NPC22192	440.48	710.51	541.37	169.15	0.00	1175.41	68.36
NPC169474	320.34	544.04	432.70	111.35	0.00	838.95	69.64
NPC297657	326.35	486.40	288.36	198.05	0.00	793.81	67.08
NPC19721	320.26	413.49	175.51	237.99	0.00	659.32	45.03
NPC279121	286.24	284.77	174.51	210.26	0.00	617.49	65.98
NPC10737	358.39	573.79	438.48	135.31	0.00	934.07	77.97
NPC137813	458.46	633.15	513.79	119.36	0.00	1037.56	56.35
NPC191146	396.44	589.20	499.81	89.39	0.00	1013.24	88.35
NPC3825	354.36	525.85	368.18	157.67	0.00	872.99	64.50
NPC161224	491.18	529.71	191.73	279.05	0.00	884.39	0.00
NPC226294	448.38	497.96	313.91	184.05	0.00	844.39	34.60
NPC19709	432.383	461.385	294.609	166.776	0.00	794.934	51.982
NPC61506	304.256	475.401	257.966	217.435	0.00	724.859	55.51
NPC107109	300.267	375.147	247.332	127.815	0.00	609.057	79.324
NPC130230	286.24	372.91	182.49	190.42	0.00	611.19	69.09
NPC175552	290.27	415.65	237.66	177.99	0.00	655.12	55.75

ROC curve plot and enrichment study

In machine learning and data mining research, Receiver Operating Characteristic (ROC) curves are usually used to study the performance of scoring function [121]. The ROC curve value detects predictive ability, which varies in the range 0.5 for random prediction to 1.0 for perfect prediction [122,123]. For analyzing the performance of our scoring functions, SP and XP data were evaluated by (ROC) curve, and the ROC curve was generated for the virtual screening experiment by employing the pose viewer (pv.- maegz) file of docking output (16 actives and 1000 decoys) of SP and XP results. SP and XP docking reveal ROC values 0.92 and 0.97, respectively (Fig. 3 ), indicating a high predictive performance of the used methods. In the later stage, the overall performance of SP and XP docking methods was measured by the values of the area under the curve(AUC). The discrimination of active and inactive compounds are measured as 100% with an AUC value equal to 1 and random with a value of 0.5 [121]. The value of the AUC was calculated as 0.91 and 0.96 for SP and XP docking, respectively. The BEDROC values at α = 160.9 for SP and XP docking were calculated as 0.658 and 0.866, respectively, indicating the better performance of XP against SP docking, as shown in Table 5 .

Fig. 3

SP and XP docking validation and representation of both actives and decoys compounds by plotting the ROC curve.

Table 5

SP and XP docking validation on the basis of their active counts and ROC value.

Docking mode	Percentage of results	1%	2%	5%	10%	20%	Enrichment		Matric value
SP docking	Active counts	6	6	11	13	14	ROC		0.92
	% of Actives	37.5	37.5	68.8	81.2	87.5	Area under accumulation		0.91
							BEDROC	160.9	0.658
								20.0	0.646
								8.0	0.752
XP docking	Active counts	9	11	13	13	16	ROC		0.97
	% of Actives	56.2	68.8	81.2	81.2	100	Area under accumulation		0.96
							BEDROC	160.9	0.866
								20.0	0.804
								8.0	0.863

MD simulation of multi-target inhibitors of SARS-CoV-2

To understand the conformational changes of the protein and ligand during the course on interactions, we have performed ten nanoseconds of MD simulation using Desmond v 5.7 for each multi-target inhibitor with one of their potent targets. Root mean square deviation (RMSD), Root mean square fluctuations (RMSF) have been analyzed throughout the simulation trajectory. The RMSD computed for NPC214620 and exonuclease complex most of the time found below 3 Å, indicating minimum conformational changes (Fig. 4A) of the complex during the simulation. The RMSD of NPC52382 protein-ligand complex presented below 1.8 Å remained within the acceptable range (0–3 Å) (Fig. 4C) indicates little induction of conformational changes of the protein during the simulation. The RMSD of NPC270578 protein-ligand complex above 3 Å (Fig. 4E) suggests the protein went substantial conformational changes during ligand interaction. To monitor the dynamic behavior of amino acid residues during protein-ligand interaction, we have analyzed the RMSF fluctuations of the simulation studies. Along with the typical N- and C-terminal fluctuations, the internal fluctuation detected for the three complexes were observed in the unstructured regions of the protein (Fig. 4B, D, and F). Next, the graphical analysis of the total number of protein-ligand interactions such as hydrogen bond, ionic interaction, hydrophobic interaction, and water bridges during the simulation for the three compounds with their respective target proteins are presented in Fig. 5 .

Fig. 4

Molecular dynamics simulation. (A) and (B) showing the RMSD and RMSF plot of Exoribonuclease and NPC214620 complex (C) and (D) representing the RMSD and RMSF of Methyltransferase and NPC52382 complex whereas (E) and (F) showing the RMSD and RMSF of Helicase and NPC270578 complex, respectively. The blue and red lines in the background of the RMSF plot showing the beta and alpha helices, respectively.

Fig. 5

Graphical representation of protein-ligand contacts. (A) Showing the interaction between the helicase and NPC270578 (B), showing the interaction between Exoribonuclease and NPC214620, whereas (C) showing the interaction between Methyltransferase and NPC52382 compound.

Molecular dynamics simulation. (A) and (B) showing the RMSD and RMSF plot of Exoribonuclease and NPC214620 complex (C) and (D) representing the RMSD and RMSF of Methyltransferase and <span class="Chemical">NPC52382 complex whereas (E) and (F) showing the RMSD and RMSF of <span class="Gene">Helicase and NPC270578 complex, respectively. The blue and red lines in the background of the RMSF plot showing the beta and alpha helices, respectively.

Graphical representation of protein-ligand contacts. (A) Showing the interaction between the <span class="Gene">helicase and <span class="Gene">NPC270578 (B), showing the interaction between Exoribonuclease and NPC214620, whereas (C) showing the interaction between Methyltransferase and NPC52382 compound.

Discussion

The highly contagious nature of SARS-CoV-2 enables the virus to spread 190 countries and territories of the globe within a short period of two and a half months. In respn>onse to the global pandemic situation, it is urgent and reasonable to employ sophisticated in silico approaches for facilitating the drug discovery for SARS-CoV-2. The highly error-prone replication [124] cycle and exceptionally shorter generation times have become a pressing issue for designing new effectual compounds for the highly evolving RNA virus. In a computational study, it has been shown that herbal leads displayed better binding potential towards probable drug targets of the Ebola virus than the known chemical analogs [125]. In the current effort, we have searched for effective inhibitors against six potent therapeutic targets of SARS-CoV-2 throughout the virtual screening on the NPASS database. Amongst 35,032 screened compounds, twenty-one compounds showed maximum docking scores compared to their respective known inhibitors. The detailed information regarding the natural source and biological activity is described in Supplementary Table 2. Notably, the antiviral properties of compounds NPC279121, NPC10737, NPC161224, NPC130230, and NPC70622 have been documented for many viruses, including Hepatitis B virus, Hepatitis C virus, Simian virus 40, Influenza A virus and HIV1 in the NPASS database. The compounds NPC214620, NPC52382, and NPC270578 are targeting five, four, and three-drug targets, respectively. These molecules are the significant marks in the multi-target-based drug design for the rapidly spreading virus. NPC214620 (PubChem CID: 11874) is 2,3,4-trihydroxybenzoic acid, identified in the course of biomarker analysis of flavonoid intake during metabolic profiling of human urine [126]. Trihydroxybenzene motif has been applied to design anti-inflammatory drugs targeting selectin molecules [127]. NPC52382 (PubChem CID: 442437) is known as neoastilbin and has 14 natural sources, including fungi, sea slug, and plants (NPASS database). The compound has biological activities with five proteins (mitochondrial Glutaminase kidney isoform, Protein disulfide-isomerase, Transcriptional regulator ERG, DNA dC->dU-editing enzyme APOBEC-3G, and Peptidyl-prolyl cis-trans NIMA-interacting1) and with two organisms (malaria parasite and human being) as per the NPASS database. According to the NPASS database, NPC270578 (PubChem CID: 119258) is known as astilbin and can be obtained from 93 natural sources and have biological activity against Xanthine dehydrogenase. Neoastilbin and astilbin are enantiomers and display the same docking score and MMGBSA value for three protein targets, including helicase, exoribonuclease, and methyltransferase as shown in Fig. 2B and C. However, unlike astilbin, neoastilbin shows better docking score and MMGBSA value than the known inhibitor (RO-7) for endoribonuclease of SARS-CoV-19 (Fig. 2C). Similar observation was indicated for the two thiocarbazate enantiomers of cathepsin L. Molecular docking studies revealed the docking score of S-enantioner (9.03 kcal/mol) more than the R-enantiomer (7.02 kcal/mol), the difference was also validated in their biological activity; S-enantiomer, with an IC50 of 56 nM and R-enantiomer, with an IC50 of 33 μM [128]. The importance of enantiomerism in the interaction of compounds with their binding sites, evaluated by docking scores, was also mentioned in other works [129,130]. On the basis of the earlier reports, it can be accepted that endoribonuclease has shown a better binding affinity for one enantiomer than the other in the present molecular docking study.

Reports from in silico work have aided in ranking the molecule/compound, thus plummeting the chances related to the poor selection of lead molecules. Therefore, such studies are worthwhile for research groups working on wet-lab experiments intending to identify novel antiviral lead compounds. Recent advances in computational approaches are playing a substantial role in the screening of drug and their design. Here, we will discuss a story of successful drug identified via structure-based virtual screening of compounds, i.e., a researcher's team successfully dis<span class="Species">covered <span class="Gene">neuraminidase inhibitor against influenza virus [131]. This is a potent, highly selective, novel, and orally active compound against influenza A and influenza B named BCX-1812, RWJ-270201 (cyclopentane peramivir). Further, the group conducted the in vitro studies and experiments on mice, rats, and ferrets and found higher potency than known drug named zanamivir and oseltamivir [131].

After <span class="Species">human clinical trials, the drug, premavir found safe and effective [132], and the team got US patent [133]. This <span class="Chemical">cyclopentane drug premavir, was approved by US-FDA to be administered as an intravenous under an Emergency Use Authorization during H1N1 pandemic in 2009 [134,135] and later on IV peramivir was approved by FDA under the brand name Rapivab TM [136]. Several other major successful drugs are existing in the market, which was discovered by in silico structure-based work includes drug (retroviral protease inhibitor) against HIV-1 [137]. Potential drugs, namely ratirexed [138], an inhibitor of thymidylate synthetase against HIV-1, and amprenavir target the protease of HIV-1 [137,139]. Other popular examples are drugs; namely, Isoniazid targets the InhA enzyme works against tuberculosis, which is also a major global health problem discovered by virtual screening [140] and Norfloxacin targets topoisomerase II & IV enzyme and treats the Urinary Tract infection [141]. All these successful drugs pave us the way to work upon the virtual screening based strategy to combat the ongoing pandemic. By considering this work now, researchers will be able to plan the experimental work on lead compounds identified in this study. However, the only concern is that there are several failure stories also reported in the documents; for instance, an antidepressant (RPX00023) was declared to be an agonist of 5-HT1A receptor while it was found to be an inhibitor [142]. Therefore, in vitro activity experiments with individual protein targets followed by SARS-CoV-2 replication assays of the compound identified from this study may support the development of novel and highly potent therapeutics for the treatment of rapidly spreading COVID-19.

Authors statement

The authors appreciate and extend thanks for the kind consideration and thorough review of this manuscript that provided us the opportunity to revise and improve it. Keeping the reviewers useful and validly raised concerns and questions in mind we revised and improved the manuscript. We are extremely thankful to the editor for the selection of five reviewers and humbly thankful to all reviewers for giving your valuable suggestion which will not only be useful for the improvement of this manuscript but also useful in the fight against minacious <span class="Disease">COVID-19 present situation. We have addressed the questions raised or corrections suggested by reviewers point by point.