Modelling studies reveal the importance of the C-terminal inter motif loop of NSP1 as a promising target site for drug discovery and screening of potential phytochemicals to combat SARS-CoV-2.

COVID-19 pandemic causative SARS-CoV-2 coronavirus is still rapid in progression and transmission even after a year. Understanding the viral transmission and impeding the replication process within human cells are considered as the vital point to control and overcome COVID-19 infection. Non-structural Protein 1, one among the proteins initially produced upon viral entry into human cells, instantly binds with the human ribosome and inhibit the host translation process by preventing the mRNA attachment. However, the formation of NSP1 bound Ribosome complex does not affect the viral replication process. NSP1 plays an indispensable role in modulating the host gene expression and completely steals the host cellular machinery. The full-length structure of NSP1 is essential for the activity in the host cell and importantly the loop connecting N and C-terminal domains are reported to play a role in ribosome binding. Due to the unavailability of the experimentally determined full-length structure of NSP1, we have modelled the complete structure using comparative modelling and the stability and conformational behaviour of the modelled structure was evaluated through molecular dynamics simulation. Interestingly, the present study reveals the significance of the inter motif loop to serves as a potential binding site for drug discovery experiments. Further, we have screened the phytochemicals from medicinal plant sources since they were used for several hundred years that minimizes the traditional drug development time. Among the 5638 phytochemicals screened against the functionally associated binding site of NSP1, the best five phytochemicals shown high docking score of -9.63 to -8.75 kcal/mol were further evaluated through molecular dynamics simulations to understand the binding affinity and stability of the complex. Prime MM-GBSA analysis gave the relative binding free energies for the top five compounds (dihydromyricetin, 10-demethylcephaeline, dihydroquercetin, pseudolycorine and tricetin) in the range of -45.17 kcal/mol to -37.23 kcal/mol, indicating its binding efficacy in the predicted binding site of NSP1. The density functional theory calculations were performed for the selected five phytochemicals to determine the complex stability and chemical reactivity. Thus, the identified phytochemicals could further be used as effective anti-viral agents to overcome COVID-19 and as well as several other viral infections.

Introduction

Animal species are predominantly susceptible to Coronaviruses (CoVs); particularly causes severe respiratory diseases in humans [1]. The first epidemic report of CoVs was registered in 2003 for Severe Acute Respiratory Syndrome (SARS) caused by a family member of CoV, Severe Acute Respiratory Syndrome CoV (SARS-CoV). Similarly, another CoV that had a major impact on the human health of the Middle Eastern region during 2012 is caused by Middle East Respiratory Syndrome CoV (MERS-CoV), which is referred as Middle East respiratory syndrome and initially reported in Saudi Arabia. A distinct variant of the human coronavirus that emerged in Wuhan in late 2019 is later named Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), which is being the causative agent for the present pandemic coronavirus disease 2019 (COVID-19). As with the other two CoVs (SARS-CoV and MERS-CoV), SARS-CoV-2 also has the possibility of zoonotic transmission to humans [2]. These CoVs belongs to the Coronaviridae family are classified under the order Nidovirales and further the sub-family Coronavirinae is broadly divided into four genera such as Alphacoronavirus (α-CoV), Betacoronavirus (β-CoV), Gammacoronavirus (γ-CoV) and Deltacoronavirus (δ-CoV) [3].

Among these four genera, α-CoVs and β-CoVs largely affect humans and also found commonly in mammals, whereas γ-CoV and δ-CoV are identified for their presence in birds. The common reservoir of CoVs are bats, which has a significant role in the evolution and distribution of CoVs within mammals. SARS-CoV-2 belongs to a β-coronavirus sub-genus and disseminated rapidly through people to people contact as well as airway transmission all over the world in early 2020. In March 2020, the World Health Organisation (WHO) declared the COVID-19 as pandemic and still the progressive infection is being reported worldwide. As of December 23, 2020, WHO reports, so far 79,931,215 cases have been clinically confirmed for COVID-19, which includes 1,765,265 deaths worldwide. The leader board of WHO statistics depicts the United States of America (18,972,813) and India (10,244,303) has more confirmed cases than other counterparts. In general, CoVs contain long positive-sense single-stranded RNA as genetic material. Comparatively, CoVs have a very long genome of about 27–32 kb in size than other RNA viruses. The 29,903 bp full-length genome of SARS-CoV-2 isolated from Wuhan-Hu-1 strain was deposited in NCBI by Wu et al. in 2020 [4]. The 29.9 kb genetic material codes for the structural and non-structural proteins (NSPs), which are essential for viral replication, transcription, structural assembly and pathogenesis.

All CoVs have a very long ORF1 in the 5′ end, which occupies almost two-thirds of the entire viral genome, composed of an overlapping two open reading frames (ORF1a and ORF1b). The translated two long precursor polyproteins 1a (pp1a) and 1 ab (pp1ab) of the viral genome in the host cell are further proteolytically cleaved to form 16 mature non-structural proteins in accordance with the position in the precursor polypeptide [5]. The majority of the NSPs are playing a crucial role in the replication and transcription of the viral genome inside the host cell. Alike other coronaviruses, SARS-CoV-2 also impede the host translation process to synthesis viral proteins instantly upon infection.

Despite a few NSPs, the majority of them share the conserved structural fold as well as sequence within the CoV family. NSP1 is one of the proteins which has sequence divergence and length difference within the family. Even-though, the amino acid sequence and length of NSP1 significantly differs within the CoVs family, the molecular function of the protein is evolutionarily conserved. Extensive studies are carried out on the NSP1 of SARS-CoV among the other family members [2,6]. The virulent nature of NSP1 and its implication in disease progression of SARS-CoV, MERS-CoV and SARS-Cov-2 mediated infections have been widely reported by several groups [[7], [8], [9], [10]]. Inducting the host mRNA degradation and inhibition of host translation are the two vital process employed by NSP1. Though the NSP1 has reported no endonucleolytic activity, it plays a role in hampering the expression of genes by captivating unknow endonuclease for host mRNA degradation. Moreover, the specificity in the degradation of host mRNA in contrast to viral RNA is still unclear. NSP1 tend to form a complex with the 40s ribosome subunit by blocking the mRNA binding site to inhibit host translation machinery [[11], [12], [13]]. Apart, NSP1 has been reported to inhibit IFN expression and also blocks several anti-viral signaling pathways, which aid the virus to evade host innate defence mechanisms. Several host proteins are found to be interacting with SARS-CoV-2 NSP1, among which four proteins namely (POLA1, POLA2, PRIM1 and PRIM2) forms the α-complex of DNA polymerase, which may lead to the antagonistic activity on the host immune system [14]. A recent study on SARS-CoV-2 NSP1 has shown to conquer the host anti-viral action by antagonizing interferon induction. Invading immune response and prevention of viral mRNA degradation in host cells is expected to be carried by N-terminal, and attachment of C-terminal domain in the mRNA binding site of 40s ribosome leads to the inhibition of host translation process [11,15,16].

Computer assisted drug designing is a powerful, fast and handy tool to explore the chemical features in the identification of potent lead molecules during crisis and pandemics. Numerous synthetic chemical molecules were designed to antagonize the key viral targets, however the ADME and toxicity profiles of the synthetic molecules in the in vitro condition are uncertain. To overcome the ADME and toxicity barriers, molecules derived from Indian herbal plants resources are explored to a larger extent as an alternative for the synthetic approach [[17], [18], [19]]. Inhibition of NSP1 has reported increasing the host innate immunity to counteract the SARS-CoV-2 infection [8,10,11,16,20,21]. The inter motif loop which connects the N-terminal and C-terminal regions are highly dynamic and facilitates the binding of C-terminal helices to the 40S ribosome unit. Due to the higher flexibility of the C-terminal inter motif loop, so far the full-length structure of NSP1 was not determined experimentally. Based on this hypothesis, we have constructed the first full-length three-dimensional structure of NSP1 with the promising drug target site inter motif loop through homology modelling and the conformational behaviour of the modelled protein are accessed through large scale molecular dynamics (MD) simulations. Subsequently, the phytochemicals were screened to identify the potent molecules using molecular docking. Further, MD simulation studies were employed to decode the protein-ligand stability in a biological mimic environment. The selected phytochemicals were analysed for the reactive groups to explore the avenues for lead enrichment to have better efficacy and high inhibition activity against NSP1.

Materials and methods

Sequence retrieval and physico-chemical properties

The full-length amino acid sequence of SARS-CoV-2 NSP1 corresponding to the first reported sequence of Wuhan-Hu-1 strain was retrieved from the NCBI database using the accession number: YP_009725297. Each one representative protein sequences of SARS-CoV-2 NSP1 corresponds to the strains isolated from Australia, Brazil, Canada, Czech Republic, France, India, Italy, Japan, Philippines, Russia, South Africa, Spain, Sri Lanka, Taiwan, Thailand, United Kingdom and the United States of America were analysed to compare the compositional difference with the reference Wuhan-Hu-1 strain isolated from China. The physico-chemical properties of the NSP1 protein of the SARS-CoV-2 Wuhan-Hu-1 strain were computed using PDB Goodies and Expasy-Protparam Tool [22,23].

Homology modelling and validation

BLAST search using the FASTA format sequence of NSP1 (YP_009725297) in NCBI-BlastP tool [24] was employed to find the structural homologs in Protein Bata Bank (PDB). The partial structures of NSP1 in the PDB repository were selected for the construction of a full-length three-dimensional structure of SARS-CoV-2 NSP1. The target and template sequence alignment were performed in Clustal omega [25]. Based on the aligned regions, the full-length three-dimensional model was constructed using the Prime tool in the Schrödinger suite [26,27]. The stereochemical [28] quality of the modelled structure was validated using the Procheck server [29].

Molecular dynamics simulation

The full-length model of SARS-CoV-2 NSP1 was subjected to MD simulation using GROMACS v4.5.5 (Groingen Machine for Chemical Simulations) package with GROMOS96 43a1 force field. In order to mimic the biological environment, Single Point Charge 216 (SPC216) water model was used to equilibrate the system within the orthorhombic shaped boundaries with the dimension of 10 × 10 × 10nm3. The methodology with the parameters described by Amala et al., 2019 [30] was employed for the apo and complex simulations of NSP1. The final production MD run of NSP1 apo structure was simulated for 100ns. The parameters and ligand topologies for the ligands were generated using PRODRG server. The final production MD run of ligand-bound NSP1 complexes was simulated for 50ns. All trajectories in both apo and complex simulations were recorded at 2ps interval. The trajectory files were analysed for the assessment of Root Mean Square Deviation (RMSD), residual fluctuations, compactness of the protein and H-bond contacts. All the trajectory values were visualized and plotted in a 2D graph using OriginPro Tool.

Receptor preparation and binding site prediction

A validated model of SARS-CoV-2 NSP1 was prepared using the Protein Preparation Wizard available in the Schrödinger modelling suite [31,32]. The addition of hydrogen atoms, assigning proper bond orders and formal charges, neutralization of side chains that are not involved in salt bridge formations and assessment of protonation states were achieved in preparation. The structure was optimized and minimized using the Optimized Potential Liquid Simulation 2005 (OPLS-2005) force field. Identification of the proper binding site of the target is essential and mandate in structure-based drug designing. SiteMap module embedded in the Schrödinger modelling suite was used to predict the possible binding pockets in NSP1 [33,34]. The SiteMap tool analyse the protein surface initially and ranks the regions which can favour the binding of ligands on the receptor. Then, computes hydrophobic and hydrophilic cores based on the amino acids on the identified protein surfaces and accordingly generate the contour maps. Further, it identifies the metal-binding sites, hydrogen bond acceptor and donor sites within the predicted hydrophilic regions. The SiteMap ranks the predictions based on the Sitescore, which is calculated from various properties such as site size and receptor degree of enclosure. A grid covering the entire binding site of the receptor was generated using the Grid Generation panel in Schrödinger. The 20 Å grid in X, Y and Z coordinates allows the small molecules to move flexibly within the defined binding site [35].

Phytochemicals collection and ligand preparation

Phytochemicals are an excellent source of alternatives to synthetic molecules and are already being in use for a few thousand years in the treatment of simple to more complex ailments. The therapeutic efficiency of Indian and Chinese natural medicines is well documented in the literatures of their respective languages and because of its very long traditional lineage, even today it is being effectively practised in several parts of the world [[36], [37], [38], [39], [40], [41], [42], [43]]. On this background, phytochemicals from Indian medicinal plants, which are reported by various research groups in literatures are included in this study. The phytochemicals reported in the literature which are available in the Pubchem database [44] were retrieved using the accession number as well as by using the phytochemicals name. The phytochemicals which are not available in Pubchem were drawn using Marvinsketch software. Apart from the literature, selected phytochemicals were retrieved from the first Indian medicinal plant database IMPPAT [45]. All the phytochemical molecule collection available for computational screening in two-dimensional form were prepared using the LipPrep module of the Schrödinger suite. The molecules were corrected for bond orders, bond lengths and a maximum of 32 conformers per molecules were generated. The conformers of the phytochemicals were optimized and energy minimized further in the LigPrep tool [46,47].

Structure-based virtual screening

Virtual screening workflow in Schrödinger suite has Qikprop and Glide (Grid-based ligand docking with energetic) modules for ADME filtration and molecular docking process. ADME (Absorption, Distribution, Metabolism and Excretion) properties of the phytochemicals were predicted in the QikProp tool [48,49]. Initially, QikProp neutralizes the phytochemicals for normal mode analysis and then computes 44 descriptive properties pertaining to QPPMDCK, QplogHERG and QPPCaco to evaluate the ADME property of the molecules. QikProp performs calculation in normal and fast modes, the normal mode analysis in QikProp allows us to calculate all the descriptive properties for effective ADME prediction compared to the fast mode. The primary ADME acceptability of the molecules was relying on Lipinski’s rule of five, an initial filtration method used to rule out the unlikely molecules from the phytochemicals collection. Software filters the molecules for screening only when it satisfying the following criteria: (a). Not more than 10 hydrogen bond acceptors, (b). Not more than 5 hydrogen bond donors, (c). Less-than 500Da molecular weight and (d). QPlogPo/w value less-than 5. Schrödinger accomplishes virtual screening in a three-stage process: (i). High-throughput virtual screening (HTVS), (ii). Standard Precision (SP) and (iii). Extra Precision (XP). HTVS protocol executes a rapid docking procedure to assess the molecules binding energy and interaction within the defined binding pocket. The top-scoring hits from the HTVS process were further subjected to SP mode followed by the final XP mode of re-docking procedures [50,51]. The final hits were ranked based on the GlideScore and the top scored molecules were selected for further complex MDs, binding free energy and DFT studies.

Prime MM-GBSA

Ligand binding and strain energies of the set of receptor-bound phytochemicals were analysed in Prime Molecular Mechanics-Generalized Born Model and Solvent Accessibility (MM-GBSA) [26,52]. Prime MM-GBSA employs OPLS-AA, SGB solvation model for polar solvation, Molecular Mechanics Energies and nonpolar solvation which includes different nonpolar SASA and van der Waals interactions. The binding free energies upon binding of phytochemicals in the receptor was calculated by ΔG bind = G complex − (G protein + G ligand)

Gcomplex = complex energy; G protein = receptor energy; G ligand = unbound ligand energy.

The Macromodel script embedded in Prime MM-GBSA calculated the energies of the atoms interacting in protein and phytochemicals.

Density functional theory

Schrödinger-Jaguar tool was used for the density functional theory (DFT) calculations. The same parameters described in Samurkas et al., 2020 [53] was replicated in geometrical optimization of the phytochemicals during DFT studies, which uses Becke’s 3 parameter exchange potential, Lee-Yang-Parr correlation function (B3LYP) and basis set 6-31G∗∗. Aqueous condition Poisson Boltzmann Solvers (PBS) was used to calculate the energies of the phytochemicals. Based on the energy values, the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) were calculated along with the difference in energy gaps.

Results and discussion

Sequence retrieval and comparative analysis

SARS-CoV-2 NSP1 sequences correspond to the respective strains isolated from various countries across the continents were retrieved using accession numbers are shown in Fig. 1 . Multiple sequence alignment of the 17 sequences with the reference sequence corresponds to the China Wuhan strain is shown in Fig. 1. Irrespective of its country representation all the residues are 100% identical to the reference NSP1 sequence except a single amino acid change (D75E) that retains the negative charge in the Taiwan NSP1 sequence. The conserved pattern observed in the multiple sequence alignment of NSP1 isolated across the world illustrates the stringent nature during the evolution. Since no significant mutations have been observed among the protein sequences across countries, NSP1 can be potentially targeted for the discovery of drugs that can hamper the rapid progression of COVID-19. Also, major challenges in the vaccine development for the virus is to tackle the rapid mutations on the target, which can be solved by targeting NSP1 for developing the vaccine [54,55]. The protein sequences of NSP1 from the alpha and beta family of coronaviruses were shown to have divergence among the families (α and β), however, partially conserved nature was reported within the family. The NSP1 sequence of SARS-CoV-2 has high sequence similarity with the SARS-CoV (84.4%), RaTG13 (96.7%) and pangolin-CoV (95.6%). However, the comparative analysis of SARS-CoV-2 NSP1 with other members of the β-CoV sub-family has very low sequence identities such as MERS-CoV: 13.9%, HCoV-HKU1: 11.1% and MHV-A59: 15.0% [8]. Notably, the absence of NSP1 in the gamma and delta family of coronaviruses shows its significant role in alpha and beta families. The high conservation in protein sequences and a very low frequency of amino acid mutations even after the worldwide transmission indicates the importance of NSP1 in SARS-CoV-2 survival and disease progression. Targeting SARS-CoV-2 NSP1 protein would be the ideal strategy to develop potent drug molecules [55], which has a very minimal chance for reduced therapeutic activity due to mutations.

Fig. 1

Multiple sequence alignment of the selected NSP1 with the reference strain isolated from Wuhan, China showing the conservation of the protein. The conserved residues were marked with ∗. The uniport accession number of the NSP1 sequences along with the geographical location of isolates are labelled.

Physiochemical properties and in silico characterization

ORF1a polypeptide of SARS-CoV-2 is fragmented into several non-structural polypeptides. The first 180 amino acids of the long polypeptide form the leader protein named Non-structural protein 1. The physiochemical properties of the NSP1 are listed in Table 1 . The predicted molecular weight of the protein is 19,775.31 Da and has the calculated isoelectric point of 5.36, which states the protein is acidic in nature. The 10-arginine and 9-lysine residues together contribute to the positive charge of the protein. Similarly, 9-aspartic acid and 18-glutamic acid residues collectively provide negative charges to NSP1. The instability index value of NSP1 is 28.83 and based on the instability index the NSP1 is categorized as a stable protein. The extinction coefficient and aliphatic index of the target protein is 12,950 and 89.72, respectively. Amino acids such as tyrosine, tryptophan, phenylalanine predominantly contribute to the extinction coefficient and to a lesser extent cystines are also accounted for. The GRAVY value of NSP1 is −0.378 and in accordance with the GRAVY value, the protein is hydrophilic in nature. NSP1 is having 22-Glycines which accounts for 12.2% of the entire composition; 21 residues of each histidine (11.7%) and valine (11.7%) present in NSP1. In the least level, each one cystine and tryptophan are in the NSP1 accounted to 0.6%. Based on the protein sequence, the NSP1 is predicted to be localized in Cytoplasm and SOUSI predicted the soluble nature of the protein.

Table 1

Physico-chemical properties of NSP1 computed based on the sequence information.

Length	180
Molecular Weight	19775.31 Da
Theoretical pI	5.36
Positively charged residues (Arg + Lys)	19
Negatively charged residues (Asp + Glu)	27
Extinction coefficients	12,950
Estimated half-life	30 h (mammalian reticulocytes, in vitro); >20 h (yeast, in vivo); >10 h (Escherichia coli, in vivo)
Instability Index	28.83 (stable)
Aliphatic index	89.72

The SARS-CoV-2 NSP1 sequence searched against the structural repository PDB using the BlastP tool resulted in identifying two partial protein structures of SARS-CoV-2 NSP1 (PDB ID: 7K7P and 6ZLW) [16,56]. The sequence of protein structure (PDB ID: 7K7P) has shown the query coverage of 65% with the identity of 100% for 118 residues from 10th position to 127th position in the N-terminal region. Similarly, the PDB ID: 6ZLW has query coverage of 18% with the identity of 100% for 33 residues from 148th position to 180th position in the C-terminal region. Based on the identity, query coverage and score, both the templates are selected for modelling the full-length three-dimensional structure of NSP1. The sequence-structure alignment of NSP1 with the identified templated structures is shown in Fig. 2 .

Fig. 2

Structure-based sequence alignment of template structures (PDB ID: 7K7P and 6ZLW) with target protein NSP1 sequence. The template region (PDB ID: 7K7P) matches with the target protein sequence are highlighted in yellow. The cyan colour depicts the region aligned between the template (PDB ID: 6ZLW) and the target NSP1 sequence. The residues that correspond to non-template regions are in pink colour. The corresponding secondary structural elements of the template and target aligned regions were pictorially represented. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

The alignment of both the template structures are accurate, the full-length structure of SARS-CoV-2 NSP1 was modelled in Schrödinger-Prime module and subsequently, the long loop (128–147) lies between the N and C-terminal domains were refined in the Prime loop refinement module. The overall structure of SARS-CoV-2 NSP1 is dissected into three core regions namely N-terminal (1–127), C-terminal (148–180) and a inter motif loop (128–147) which connects the N-and C-terminal regions (Fig. 3 ).

Fig. 3

a). Modelled structure of full-length NSP1 showing the long loop (dark blue) that connects the C-terminal helical motif with the N-terminal region. b). Superposition of modelled protein structure (orange) with the templates PDB ID: 7K7P (cyan) and 6ZLW (pale green). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

The modelled protein contains three helices (α1, α2 and α3), out of which α1 (36–48) present in the N-terminal and remaining two helices (α2: 153–159 and α3: 166–178) in the C-terminal region. The structure of NSP1 composed of seven strands namely β1 (12–20), β2 (51–55), β3 (67–74), β4 (83–92), β5 (94–97), β6 (102–110) and β7 (116–125) and all of the strands are located in the N-terminal region. The evaluation of modelled structure in the Ramachandran plot generated using Procheck (Supplementary Fig. 1) shown 89.3% of residues in favoured regions, 10.1% of residues in additionally allowed regions, 0.7% of residues in the generously allowed region and none of the residues were mapped in the disallowed region [57]. The evaluation of protein structure in the Ramachandran plot without errors shows that the good quality of the NSP1 modelled structure for further analysis.

Molecular dynamics simulation of full length NSP1

The homology modelling method generates a static model based on the template information, MDs was performed to evaluate the conformational changes and stability of NSP1 structure in the dynamic mode. RMSD, Root Mean Square Fluctuation (RMSF), Rg and solvent accessible surface area of the protein during the simulation period were calculated from the MD trajectories. RMSD, RMSF and Rg serve as the principal factors in measuring the conformational and structural stability of the protein model. An average displacement of backbone Cα atoms from the initial structure is referred as RMSD and the average RMSD calculated from whole trajectories compared with the initial structure is plotted in Fig. 4 a. Increased RMSD pattern was observed up to 5ns until the model structure accommodates well in the environment. Later in the simulating environment, protein has minimal movement throughout the timeframe with an average RMSD of 0.76 nm and the SD of 0.11 nm. The graphical inspection of the values in the RMSD plot shows the modelled protein structure is quite stable throughout the 100ns simulation.

Fig. 4

The conformational stability analysis of NSP1 through 100ns molecular dynamics simulation. a). RMSD - average displacement of Cα backbone atoms plotted against the time scale. b). RMSF - fluctuation of residues over the time of simulation. c). Rg - compactness of the protein detected during the 100ns simulation.

The residual movement of NSP1 in a dynamic environment was monitored by analysing the RMSF and the values are plotted in the graph as shown in Fig. 4b. The amino acids in the helical regions of the protein shown minimal fluctuations, whereas the residues in loop regions shown higher fluctuations throughout the simulation. The residues fluctuated on an average scale of 0.23 nm during the entire simulation period with a SD of 0.11 nm. Most significant fluctuations on the higher scale of 0.6 nm were observed for residues in the loop which connects the N and C-terminal regions. The C-terminal helices have also fluctuated to the scale of 0.5 nm during the simulation, which shows the flexible property of the C-terminal region to bind into the ribosome unit. Several disordered regions were predicted in the C-terminal regions [58] and the disordered linker regions connecting the N and C-terminal domains has a high propensity to form various conformations of the protein due to its dynamic nature [59]. The flexibility of the linker loop facilitates the swing of the C-terminal domain, which pivotally binds in the mRNA binding channel of the ribosome [59]. Rg is computed to evaluate the compactness of the protein-based on mass-weighted RMSD of atoms from the central core geometry. The protein Rg was estimated to be 1.47 nm (SD: 0.02 nm) on average and the lesser Rg depicts the higher protein compactness (Fig. 4c). The NSP1 structure remained stable at 300K temperature during the simulation process was evident from the less Rg value. The model structure is also validated by computing the energy profiles. The simulation result shows that the model structure is reliable and quite stable in terms of conformation stability, which can be used for structure-guided screening of inhibitors.

Receptor preparation and binding site analysis

Receptor preparation is a mandated step before the structure-based virtual screening process since the results are completely relying on the input receptor structure. The modelled protein was prepared by correcting the bond orders and adjusting bond distances to minimize the structural deformities and also to avoid steric clashes. The modelled structure was optimized initially using the OPLS 2005 force field to restrain the geometrical orientations. The rotating hydroxy and thiol hydrogen atoms were optimized, protonation and tautomerization states were generated for the aspartic acid, glutamic acid, and histidine residues and chi flips were performed for asparagine, glutamine and histidine residues. Finally, the optimized protein was minimized to attain the stable biologically relevant conformation. Prediction of the proper binding site is essential for the identification of selective and potent inhibitors in structure-guided drug designing. SiteMap predicted the probable residues which could potentially involve in the binding of small molecules. Top scored NSP1 binding site has the SiteScore, drugability score and volume of 0.944, 0.911 and 495.978, respectively. The residues predicted to form the binding site are Lys72, Glu96, Arg99, Ser100, Gly101, Glu102, Ala131, Gly133, His134, Ser135, Tyr136, Gly137, Ala138, Asp139, Leu140, Lys141, Phe142, Asp147, Glu148, Leu149, Gly150, Thr151, Asp152, Tyr154, Glu155 and Glu159. The binding cavity lies between the N- and C-terminal regions of NSP1 and the residues in the long loop severs as the binding site residues. The hydrogen bonds and van der Waals interactions between the receptor and ligand were formed by the polar residues in the binding sites, whereas hydrophobic residues in the receptor involve in the formation of pi-pi stacking and pi-cation interactions with the small molecule. Hence, the predicted functionally involving amino acids are considered to generate the Gird panel of 20 × 20 × 20 Å for screening inhibitors.

Virtual screening of molecules

Computation screening of molecules from the large collection of ligands against the drug targets are widely used to achieve the identification of hit molecules in a short time and currently, it fastens the drug discovery process by 50%. A total of 5638 phytochemicals from various Indian medicinal plants were prepared in Ligprep module. Conformational sampling of phytochemicals is essential to predict the proper binding and accurate orientation of molecules within the receptor binding site. Also, various conformers of the phytochemicals were generated during the preparation process, which increases the possibility of phytochemicals binding in NSP1 protein. The molecules which passed the stringent thumb rule of Lipinski’s were subjected to predict the presence of reactive groups. The filtering and ligand preparation process has generated 11,609 conformations for 5638 phytochemicals and all these 11,609 conformations were subjected to virtual screening process. Initially, 11,609 conformations were randomly docked into the predicted binding site of NSP1 through HTVS mode. The top 10% of hits (1,160) were obtained from the HTVS mode was selected for the second round of refined docking process referred as SP mode. The 10% of hits (116 molecules) obtained from the SP was finally subjected to XP mode. The XP mode docking pipeline predicts the accurate binding affinity and interaction between the phytochemicals and NSP1. Based on the GlideXP score, XP mode returned 30% of hits (34 molecules) at the end of the virtual screening pipeline. Among these 34 molecules, top-scoring five phytochemicals molecules which are showing higher GlideXP score compared to other phytochemicals are represented in Table 2 and its 2D interactions profiles are shown in Fig. 5 . The top-scoring phytochemicals are dihydromyricetin (CID161557), 10-demethylcephaeline (CID185699), dihydroquercetin (CID10185), pseudolycorine (CID443689) and tricetin (CID5281701). The binding mode and interactions of the top-scored phytochemicals were discussed in binding mode analysis.

Table 2

Top scored phytochemicals identified from virtual screening and their interactions (Hydrogen bonds and salt bridges) with protein atoms.

ID	Glide Score (kcal/mol)	Glide Energy (kcal/mol)	Ligand and Protein Interactions	Bond length (Å) and its binding strength	ΔG Bind (kcal/mol) [MM-GBSA]
CID161557	−9.68	−51.03	C–O⋯HN (LYS72)OH⋯O <svg xmlns="http://www.w3.org/2000/svg" version="1.0" width="20.666667pt" height="16.000000pt" viewBox="0 0 20.666667 16.000000" preserveAspectRatio="xMidYMid meet"><metadata> Created by potrace 1.16, written by Peter Selinger 2001-2019 </metadata><g transform="translate(1.000000,15.000000) scale(0.019444,-0.019444)" fill="currentColor" stroke="none"><path d="M0 440 l0 -40 480 0 480 0 0 40 0 40 -480 0 -480 0 0 -40z M0 280 l0 -40 480 0 480 0 0 40 0 40 -480 0 -480 0 0 -40z"/></g></svg> C (GLN96)OH⋯OC (GLY137)C=O⋯HN (ASP139)OH⋯OC (GLU148)OH⋯OC (GLU155)	2.31 (Strong)1.73 (Strong)1.98 (Strong)1.95 (Strong)2.56 (Moderate)1.77 (Strong)	−42.01
CID185699	−9.18	−46.67	OH⋯OC(LYS141)NH⋯O–C (ASP152)HN⋯O–C (ASP152)	1.81 (Strong)2.03 (Strong)2.95 (Salt bridge)	−44.92
CID10185	−9.11	−49.73	C–O⋯HN (LYS72)OH⋯OC (GLN96)OH⋯OC (GLY137)OH⋯OC (ASP139)NH⋯OC (ASP139)OH⋯OC (GLU155)OH⋯OC (GLU155)	2.35 (Strong)1.74 (Strong)2.04 (Strong)2.34 (Strong)1.95 (Strong)2.64 (Moderate)1.79 (Strong)	−37.23
CID443689	−8.71	−44.31	HO⋯HN (ARG99)OH⋯OC (GLY137)OH⋯OC (GLY137)C–O⋯HN (ASP139)C–O⋯HN (LYS141)HN⋯C–O (ASP152)	2.20 (Strong)2.05 (Strong)1.71 (Strong)2.38 (Strong)2.18 (Strong)4.33 (Salt bridge)	−45.17
CID5281701	−8.69	−46.86	C–O⋯HN (LYS72)OH⋯OC (GLN96)OH⋯OC (GLY137)CO⋯HN (ASP139)OH⋯OC (GLU148)OH⋯OC (GLU155)OH⋯OC (GLU155)	2.22 (Strong)1.72 (Strong)2.23 (Strong)2.02 (Strong)2.50 (Strong)2.31 (Strong)1.82 (Strong)	−42.89

Fig. 5

a) Superposition of top five ranked phytochemicals in the binding pocket of NSP1, which encapsulated by the long inter motif loop. 2D profiles of five top-scored phytochemicals, b) dihydromyricetin, c) 10-demethylcephaeline, d) dihydroquercetin, e) pseudolycorine and f) tricetin, showing interaction pattern with NSP1 residues. The inter motif loop residues Gly137, Asp139, Lys141 and Glu148 were predominantly found in the ligand interactions. The H-bonds formed between the phytochemical and NSP1 were shown in the arrow, whereas the line shows between the atoms of phytochemicals and amino acids represent the salt bridges.

Binding mode analysis

The strength of the molecule binding in the receptor is determined by the hydrogen bond distance [60,61]. Hydrogen bond formed between the acceptor and donor atoms at a distance less than or equivalent to 2.5 Å are stronger bonds, which has covalent characteristics. Whereas the hydrogen bonds formed in the distance between 2.5 Å to 3.2 Å are moderate and have electrostatic characters [62]. The dihydromyricetin (CID 161557) has a higher GlideXP score of −9.68 kcal/mol, which has strong interaction with the receptor protein NSP1 through six hydrogen bonds interactions (Lys72, Gln96, Gly137, Asp139, Glu148 and Glu155). The H-bond was established between the hydrogen atom in the amine group of Lys72 with the oxygen atom in the carboxyl group of dihydromyricetin at a distance of 2.31 Å. The H-bond with a distance of 1.73 Å has formed between the carboxyl oxygen in Gln96 and hydrogen atom in dihydromyricetin hydroxyl group. Similarly, the carboxyl group oxygen atom of Gly137 formed H-bond at a distance of 1.98 Å with the hydrogen atom in hydroxyl of dihydromyricetin. Asp139 has interacted with the oxygen atom of the carboxyl group in dihydromyricetin through the atom hydrogen in its amino group at a distance of 1.95 Å. A bond between the Glu148 carboxyl oxygen atom and dihydromyricetin hydroxyl hydrogen atom was formed at the distance of 2.56 Å. Glu155 has established a H-bond (distance: 1.77 Å) between the carboxyl oxygen atom and hydroxyl group of dihydromyricetin. Moreover, the dihydromirycetin has been already in clinical trial Phase 2 for its activity against glycemic control and secretion of insulin in Type 2 diabetes mellitus (Clinical Trial Identifier: NCT03606694). Also, dihydromirycetin exhibits potential antimicrobial properties along with antioxidant activity were experimentally reported earlier [9,63].

The 10-demethylcephaeline (CID185699) possess the GlideXP score of −9.18 kcal/mol and the phytochemical tightly trapped within the binding pocket of NSP1 by two H-bonds (Lys141 and Asp152) and a salt bridge (Asp152). Hydrogen in the hydroxyl group of 10-demethylcephaeline established H-bond at the distance of 1.81 Å with the oxygen atom in Lys141 carboxyl group. The hydrogen atom in the amine group of 10-demethylcephaeline formed H-bond (distance: 2.03 Å) interaction with carboxy oxygen in Asp152. Additionally, Asp152 has formed a salt bridge (C–O⋯HN) with 10-demethylcephaeline at a distance of 2.95 Å. The biological activity of 10-demethylcephaeline isolated from Alangium longiflorum has shown promising anticancer activity against lung carcinoma and breast adenocarcinoma [64].

The (2r)-2-(3,4-dihydroxyphenyl)-3,5,7-trihydroxy-2,3-dihydro-4h-chromen-4-one (CID10185-also referred as dihydroquercetin) has the GlideXP score of −9.11 kcal/mol. Dihydroquercetin bound within the pocket of NSP1 protein by establishing seven H-bonds (Lys72, Gln96, Gly137, Asp139, Glu155). Lys72 amine group forms a hydrogen bond in the distance of 2.35 Å with carboxyl group oxygen of dihydroquercetin. A hydrogen atom in the hydroxyl group of dihydroquercetin forms an H-bond with a carboxyl oxygen atom in Gln96 (distance of 1.74 Å). Similarly, hydroxyl group hydrogen atom in dihydroquercetin formed H-bonds with carboxyl oxygen atoms in Gly137 (distance: 2.04 Å) and Asp139 (distance: 2.34 Å). Asp139 carboxyl group oxygen atom has also formed H-bond in the distance of 1.95 Å (NH⋯OC) with amine group in dihydroquercetin. Glu155 has formed two H-bonds with dihydroquercetin in the distance of 2.64 Å (OH⋯OC) and 1.79 Å (OH⋯OC), respectively. dihydroquercetin has been shown to have higher activity against several diseases such as inflammation, hepatitis, cardiovascular disease, cancer, oxidative cellular injury and liver disease [63,65]. Also, recently a review on the use of antioxidants in the treatment of COVID-19 has briefed the potentials of dihydroquercetin and its significant impact on the SARS-CoV treatment [66].

Pseudolycorine (CID443689) exhibited the GlideXP score of −8.71 kcal/mol and showed five H-bond interactions with Arg99, Gly137, Asp139 and Lys141. Arg99 amine group hydrogen interacts with the hydroxyl group oxygen of pseudolycorine. The distance between the hydrogen of the amino group and the oxygen of the hydroxyl group is 2.20 Å. Hydrogen atoms in hydroxyl groups of pseudolycorine have formed two individual H-bonds with an oxygen atom in the carboxyl group of Gly137 in the distance of 2.05 Å and 1.71 Å, respectively. Amine hydrogen in Asp139 interacts with a hydroxyl group in phytochemical through H-bond in the distance of 2.38 Å (C–O⋯HN). Similarly, hydrogen in the amino group of Lys141 favoured the H-bond formation with carboxyl oxygen at a distance of 2.18 Å (C–O⋯HN). In addition to the five H-bonds, pseudolycorine has established one salt bridge with Asp152 at a distance of 4.33 Å (HN⋯C–O). The therapeutic potential of pseudolycorine in cancer, skin diseases and bactericides was experimentally proven by various research groups. The toxic properties of pseudolycorine were also been evaluated experimentally [[67], [68], [69]].

Tricetin (CID5281701) was strongly connected with the receptor NSP1 protein through seven hydrogen bonds (Lys72, Gln96, Gly137, Asp139, Glu148 and Glu155) and the binding pose has resulted in the GlideXP score of −8.69 kcal/mol. Amine group hydrogen in Lys72 has formed H-bond at the distance of 2.22 Å with carboxyl oxygen in tricetin. Gln96 carboxyl group interacts with a hydroxyl group in tricetin, where the oxygen atom of amino acids interacts with a hydrogen atom at the distance of 1.72 Å. The hydrogen of the hydroxyl group in tricetin established H-bond with 2.23 Å distanced carboxyl group in Gly137 (OH⋯OC). 2.02 Å H-bond is formed between Asp139 and oxygen of tricetin (CO⋯HN). Glu148 carboxyl group interacts with the hydroxyl group of tricetin and forms H-bond (OH⋯OC) at a distance of 2.50 Å. Two hydroxyl groups hydrogen atoms 31 and 32 of tricetin interacts with the carboxyl oxygen atom of Glu155 in the distance of 1.82 Å and 2.31 Å, respectively. The pharmacological properties of the phytochemical tricetin were reported to have potential effects in anti-inflammation, anti-diabetic and also regulate the neoplasm metastasis [[70], [71], [72]]. The dietary supplemented tricetin has shown to suppress the metastasis from non-small cell lung to bone [73]. The experimental evidence of all these five phytochemicals have been reported against various diseases and the computational ADME values support it. Hence, these five phytochemicals can be directly used as anti-viral molecules.

The docking protocol was validated by computing the relative binding free energies of the docked protein-ligand complex [74]. Prime MM-GBSA calculates various energy components of the protein, ligand and protein-ligand complex. The energy values include minimized energy, surface area energy and solvation energy [75]. The phytochemicals docked in NSP1 was analysed to estimate the relative binding affinities of the respective molecules in the receptor-binding site. Interaction of the ligand with the receptor or formation of a bond between the atoms resulting in the formation of binding energy (ΔGbind), indirectly shows the stability of the protein-ligand complex. Prime MM-GBSA analysis calculated the relative binding free energy of the phytochemicals bound NSP1 complexes in the range between −45.17 kcal/mol to −37.23 kcal/mol (Table 2). The molecule CID443689 has shown a higher binding free energy value of −44.23 kcal/mol, whereas CID185699 molecule has a relatively similar binding energy of −44.0 kcal/mol. Among the five selected phytochemicals, CID10185 has shown the least binding free energy of −38.94 kcal/mol. Even-though the binding free energy of CID16157 is comparatively lesser than CID443689 and CID185699, it has shown higher docking energy and score, which reflects the strong binding compared to other molecules. The results of the binding free energy calculation showed the coulomb (ΔGCoulomb) and van der Waals (ΔGvdW) energies have significantly contributed to the ligand binding. Three snapshots from the trajectories of the simulation were collected and prime MM-GBSA calculations were performed. The complex (1–5) showed the average free energy values of −50.49 kcal/mol, −51.23 kcal/mol, −38.58 kcal/mol, −41.99 kcal/mol and −42.71 kcal/mol, respectively. Our results show that there are no significant changes in the binding free energy between the initial docking and simulated complex, which depicts the strong binding of ligand molecules in the binding pocket of NSP1.

DFT analysis

All five selected phytochemicals were optimized using the B3LYP/6–31∗∗G basic set embedded in PBF solvation was analysed based on the frontier molecular orbital distribution, which predicts the likely regions responsible for the active interaction with receptor protein NSP1. Molecular features of the molecules computed in terms of HOMO, LUMO and molecular electrostatic contours were significantly correlate with the biological activity [76]. HOMO and LUMO play an indispensable role in charge transfers during the chemical reactions. The predicted energy values of the HOMO, LUMO and energy gap are shown in Fig. 6 . The HOMO (EHOMO) values of the selected phytochemicals were ranges from −0.2215eV to −0.2178eV. HOMO values of the phytochemicals directly relate to the propensity of strong binding in the receptor protein and higher HOMO values depict the ionization potential. The LUMO energies (ELUMO) computed for the phytochemicals are in the range from −0.0649eV to −0.0154eV. The low LUMO energies indicate the electron affinity of the phytochemicals. The HOMO and LUMO energy contours likely display the electrophilic and nucleophilic attack sites of the molecules. The minimal differences between the HOMO and LUMO energies depict the chemical reactivity, polarizability and stability of the molecule.

Fig. 6

3D molecular orbital contour map of the phytochemicals. The Highest Occupied Molecular Orbital (HOMO) of the five phytochemicals represents the most active and least active regions in (a), (c), (e), (g) and (i). Lowest Unoccupied Molecular Orbital (LUMO) representing the 3D contour maps of most active and least active regions in (b), (d), (f), (h) and (j). The differences between the LUMO and HOMO energies are shown as Energy Gap.

The gap (ELUMO- EHOMO) calculated from the differences in LUMO and HOMO energies of the selected phytochemicals are in the range of 0.2025eV–0.157eV. The distribution plot shows that the HOMO and LUMO energies are occupied in the distal ends of the phytochemicals. The better inhibitors are readily donating the electrons to unoccupied electron regions in the residues and also willing to receive the free-electron from receptor residues. A minimal energy gap could pave the ways for the excitement of electron from HOMO to LUMO and vice-versa, which depicts the strong binding affinity between the NSP1 and phytochemicals.

NSP1-ligand complex dynamics simulation

The inhibition activity of the identified phytochemicals completely relies on the stable binding within the receptor cavity. The NSP1-phytochemical complexes were examined in dynamic simulation to explore the binding nature and stability. The trajectories of the 50ns simulations of the complexes were analysed and the corresponding RMSD and RMSF were plotted is shown in Fig. 7 a. The backbone RMSD of five complexes (Complex 1: NSP1+CID16157; Complex 2: NSP1+ CID185699; Complex 3: NSP1+ CID10185; Complex 4: NSP1+CID443689; Complex 5: NSP1+CID5281701) were analysed for comparative analysis. The backbone RMSD of the protein-ligand complexes shows initial raise during the simulation process and after 5ns, the complexes attained stable equilibrium. Among the five complexes, Complex 1, Complex 2, Complex 3 and Complex 5 were having an average RMS deviation of 0.6 nm from the initial docked complex structure, whereas complex 4 has shown comparatively lesser RMS deviation from the initial conformation to the average scale of 0.3 nm. Average RMS deviations of the Complex 1, Complex 2, Complex 3, Complex 4 and Complex 5 are 0.59 nm (SD - 0.04 nm), 0.51 nm (SD - 0.03 nm), 0.61 nm (SD - 0.07 nm), 0.2 nm (SD - 0.05 nm) and 0.59 nm (SD - 0.05 nm), respectively. NSP1 structure has attained significant stability upon binding of the ligand, which is evident from the comparison of the RMSD plot (apo and complex). The apo form of the protein has shown the RMS deviation on an average scale of 0.8 nm, whereas the phytochemical bound NSP1 have shown a lesser deviation.

Fig. 7

Stability analysis of the top-scored phytochemicals bound NSP1 complexes. (a) Comparative analysis of five phytochemical bound to NSP1 protein based on the backbone RMS deviations. (b) Residual fluctuations of NSP1 upon binding of selected five phytochemicals. (c) The compactness of the phytochemical bound NSP1 was mapped in the Rg graph.

Comparison of backbone RMSD of protein-ligand complex and RMSD of ligand molecule indicates that the complex showed an increase in RMSD up to 0.61 nm, whereas ligand RMSD were decreased and reached a maximum up to 0.2 nm. The decrease in RMSD indicates a stronger binding of the ligands into the binding site of the NSP1. Average RMSD of the ligands CID16151, CID185699, CID10185, CID443689 and CD5281701 are 0.05 nm (SD-0.01 nm), 0.14 nm (SD-0.03 nm), 0.13 nm (SD-0.03 nm), 0.03 nm (SD-0.00 nm) and 0.07 nm (SD-0.02 nm), respectively. Interestingly, the ligand CID443689 maintained stable RMSD throughout the MD from the initial simulation, whereas all other ligands showed initial drift in their RMSD values and then reached stability in the remaining simulation. Also, the ligands CID5281701 and CID185699 showed maximum deviation throughout the simulation. The ligand RMSD deviations during the simulation period are shown in Supplementary Fig. 3. Our results show that the phytochemicals strongly bound in the NSP1 protein and have stable interactions, which restricted the movement of ligands during the simulations.

The residue level fluctuations were analysed to decode the local changes in the protein structure. RMSF values of the protein complexes were shown in Fig. 7b. The average residual fluctuations of the Complex 1 was 0.17 nm with a SD of 0.11 nm. Complex 2 has fluctuations of residues to 0.20 nm in average and the SD of the Complex 2 is 0.11 nm. Similarly, complexes 3, 4 and 5 have shown average fluctuations of 0.18 nm, 0.23 nm and 0.18 nm with the corresponding SD of 0.15 nm, 0.18 nm and 0.09 nm, respectively. The free loop residues in the N-terminal region has shown higher fluctuations up to 0.7 nm, which had less movement in apo simulation. The long loop regions and their connecting secondary structures on either side of the protein have shown higher movement in the apo simulation has limited fluctuations in the phytochemical bound complex due to the strong binding of phytochemicals. The residual fluctuations of C-terminal helices in apo-protein up to 0.7 nm was highly stabilized upon tight binding of ligand and the movement of residues limited to 0.3 nm on average.

The Rg depicts the core compactness and its stability of protein structure upon ligand binding. The Rg graph of the protein-ligand complexes are shown in Fig. 7c. The average Rg of the complexes (1–5) are 1.49 nm, 1.49 nm, 1.48 nm, 2.51 nm and 1.49 nm. Similarly, the SD values of the complexes (1–5) are 0.03 nm, 0.02 nm, 0.02 nm, 0.03 nm and 0.01 nm. The minimal SD shows the protein-ligand complexes are highly stable during the simulation process and a comparison with apo-protein shows that the ligand binding on the flexible loop region stabilized the core structure of the protein.

The H-bonds formed between the NSP1 and phytochemicals were analysed, which are mainly responsible for the stability of the protein-ligand complexes and moreover these H-bonds contribute substantially in specificity and ADME. The H-bond interaction profiles of the complexes during the 50ns simulation are shown in Supplementary Fig. 3. The complex 1 has shown to the maximum of eight H-bonds during the simulation process. Five hydrogen bond interactions were observed throughout the simulation process of Complex 2, Complex 3, Complex 4 and Complex 5. Furthermore, a sixth hydrogen bond was observed at several instances in Complex 3 and Complex 5 during the simulation process. The strong H-bond interactions between the NSP1 and phytochemicals were the reason for the stable complex and the binding of lead molecules relatively stabilise the protein structure as well.

Conclusions

Direct therapeutics targeting the coronavirus to control COVID-19 is still yet to be discovered and moreover, the immediate requirement of drugs is not achievable through a regular drug discovery pipeline. Computational methods are widely used to speed up the drug discovery process to some extent; however, the drug testing and trial phases are time-consuming. Repurposing of existing drugs or phytochemicals being practised for several hundred years in Indian and Chinese medicines were the only available options to minimize the time and immediate usage. The clinical isolates from various geographical locations showed no significant mutation in NSP1 during viral transmissions serves as the prominent drug target to develop anti-viral molecules. The flexible inter motif loop connecting N and C terminal regions were found vital for the NSP1 biological function and play a significant role in ribosome binding. The present study decoded the potential inter motif loop as a promising drug binding site of NSP1 and screened for the potential inhibitors from the phytochemicals available in medicinal plants. The selected five phytochemicals are having better binding affinity and higher stability along with good pharmacokinetic properties. Moreover, the top-scoring phytochemical dihydromyricetin is already in Phase 2 Trials for diabetes, similarly, the other four phytochemicals have been experimentally evaluated by different research groups against various cancers. Since the experimental evaluation of these phytochemicals are done earlier, it can surpass the preliminary testing process and facilitate directed testing of anti-viral activity. Hence, these phytochemicals could be used as potent molecules against COVID-19 and other viral infections.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.