Repurposing potential of Ayurvedic medicinal plants derived active principles against SARS-CoV-2 associated target proteins revealed by molecular docking, molecular dynamics and MM-PBSA studies.

All the plants and their secondary metabolites used in the present study were obtained from Ayurveda, with historical roots in the Indian subcontinent. The selected secondary metabolites have been experimentally validated and reported as potent antiviral agents against genetically-close human viruses. The plants have also been used as a folk medicine to treat cold, cough, asthma, bronchitis, and severe acute respiratory syndrome in India and across the globe since time immemorial. The present study aimed to assess the repurposing possibility of potent antiviral compounds with SARS-CoV-2 target proteins and also with host-specific receptor and activator protease that facilitates the viral entry into the host body. Molecular docking (MDc) was performed to study molecular affinities of antiviral compounds with aforesaid target proteins. The top-scoring conformations identified through docking analysis were further validated by 100 ns molecular dynamic (MD) simulation run. The stability of the conformation was studied in detail by investigating the binding free energy using MM-PBSA method. Finally, the binding affinities of all the compounds were also compared with a reference ligand, remdesivir, against the target protein RdRp. Additionally, pharmacophore features, 3D structure alignment of potent compounds and Bayesian machine learning model were also used to support the MDc and MD simulation. Overall, the study emphasized that curcumin possesses a strong binding ability with host-specific receptors, furin and ACE2. In contrast, gingerol has shown strong interactions with spike protein, and RdRp and quercetin with main protease (Mpro) of SARS-CoV-2. In fact, all these target proteins play an essential role in mediating viral replication, and therefore, compounds targeting aforesaid target proteins are expected to block the viral replication and transcription. Overall, gingerol, curcumin and quercetin own multitarget binding ability that can be used alone or in combination to enhance therapeutic efficacy against COVID-19. The obtained results encourage further in vitro and in vivo investigations and also support the traditional use of antiviral plants preventively.

Introduction

Beginning on 17 November 2019 (media reports on unpublished Chinese government data), a novel coronavirus from Hubei's capital Wuhan, China, designated as SARS-CoV-2, has caused an international outbreak of serious and infectious respiratory illness termed COVID-19 [1]. Until September, 30th, 2020, the pandemic virus [2] has infected at least 33.8 million people, and killed 1.01 million worldwide. The COVID-19 manifests from mild cough, fever, self-limiting respiratory tract illness to severe progressive pneumonia in both lungs which can lead to multi-organ failure, and ultimately death occurs [3], [4].

The absence of explicit treatment for COVID-19 leads the population over many regions in the world to use medicinal herbs known in ethnopharmacology as antiviral. Teams around the world are working hard to develop active substances against SARS-CoV-2. This paper reports the state of-the-art antiviral choices against COVID-19, with a specific focus on the recent advances of Ayurvedic antiviral medication repurposing procedures and their potential therapeutic suggestion [5]. As per modern Ayurvedic literature, the root of Ayurveda can be traced back to around 6000 BCE when it originated as an oral tradition known as the “Mother of All Healing” [6]. Several clinical drugs with limited or no observable side effects have been developed and practiced from Ayurveda since ancient times to modern practice as ‘tradition to trend’ [7]. Until today, the capability of Ayurvedic prescriptions has not been completely researched and needs further investigation with modern scientific validation approaches for better therapeutic leads.

The selected bioactive compounds from various plants recorded in Ayurveda with potent antiviral activity against different types of viruses [8], [9], [10], [11] are detailed in Table 1. Several researchers also reported the traditional use and research findings of all the plants related to antiviral activities; Allium sativum [12], [13], Allium cepa [14], Zingiber officinale [15], [16], Syzygium aromaticum [17], Mentha piperita [18], [19], Foeniculum vulgare [20], Ocimum sanctum [21], Origanum vulgare [22], [23], Curcuma longa [21], Tinospora cordifolia [24], [25], Cinnamomum cassia [26] and Cordyceps militaris [27].

Table 1

Antiviral Ayurvedic medicinal plants with their bioactive compounds tested against closely related human viruses with the reported mode of action (NCBI-Pubmed literature).

Sl. nos.	Plants	Bioactive principles	Antiviral activity tested on	IC50	Mechanism of action	References
1.	Allium sativum L. (Garlic)	Allicin, diallyl trisulfide, ajoene and alliin	Infectious bronchitis virus (Coronavirus), HIV-1, Influenza B, Human cytomegalovirus (HCMV), Parainfluenza virus type 3, Herpes Simplex type 1 and 2, Vaccinia virus, Dengue virus (DENV-2), Vesicular stomatitis virus and Human rhinovirus type 2	10–34 μM (HIV-1)	TNF-α inhibitor, Immunomodulatory and anti-inflammatory.	[8], [9], [10], [11]
2.	Allium cepa L.(Onion)	Quercetin, zalcitabine, allicin and ribavirin	Poliovirus and Dengue virus type -2 (DENV-2)	35.7 μg/mL (DENV-2)	RNAs (-) cannot be transcribed to generate more copies of positive polarity.	[14], [84], [85]
3.	Zingiber officinale Roscoe(Ginger)	Gingerol and zingerone	H1N1 Flu (Swine flu), Influenza, Anti-avian influenza virus H9N2 and Human respiratory syncytial virus	73.3 mg/mL (HRSV)	Enhance the secretion of IFN-β.	[86], [87], [88]
4.	Syzygium aromaticum L. (Clove)	Eugenol, E-cinnamaldehyde, carvacrol and thymol	Hepatitis C viruses and Herpes simplex virus (HSV-1 and HSV-2)	25.6 μg/mL (HSV-1) and 16.2 μg/mL (HSV-2)	Inhibition of viral replication, damage to viral envelopes of freshly formed virions.	[89], [90], [91]
5.	Mentha piperita L. (Peppermint)	Menthol and rosmarinic acid	Herpes, Influenza A, Newcastle disease virus, herpes simplex virus (HSV-1) and Vaccinia virus	0.25 μg/mL (HSV-1)	Inhibits specifically the viral DNA polymerase during the replication cycle.	[80], [92]
6.	Foeniculum vulgare Mill. (Fennel)	Trans-anethole	Herpes simplex virus (HSV-1) and Parainfluenza type-3 (PI-3)	0.025 μg/mL (HSV-1)	Not reported	[20], [80], [93]
7.	Ocimum sanctum L. (Holy Basil or Tulsi)	Apigenin	Herpes simplex viruses (HSV), adenoviruses (ADV), Hepatitis B virus, Coxsackievirus B1 (CVB1), Enterovirus 71 (EV71) and H1N1Flu virus	39–55 μg/mL (H1N1)	Increase in the levels of IFN-γ, IL-4 and percentages of T-helper and NK-cells.	[94], [95]
8.	Origanum vulgare L. (Oregano)	Carvacrol	Acyclovir-resistant herpes simplex virus type 1 (ACVR-HHV-1), Acyclovir-sensitive HHV-1, Human respiratory syncytial virus (HRSV), Bovine herpesvirus type 2 (BoHV-2), Bovine viral diarrhea virus (BVDV), Herpes simplex virus type 1 (HSV-1), murine norovirus (MNV), Canine distemper virus (CDV), Canine adenovirus (CAV), Canine coronavirus (CCoV), Respiratory syncytial virus (RSV) and Coxsackievirus B3 (CVB3)	23.1 μM (RSV long strain) and 81.7 μM (RSV A strain)	Capsid protein disintegration.	[96], [97], [98]
9.	Curcuma longa L. (Turmeric)	Curcumin	Parainfluenza virus type 3 (PIV-3), Feline infectious peritonitis virus (FIPV), Vesicular stomatitis virus (VSV), Flock house virus (FHV), Kaposi's sarcoma virus, Respiratory syncytial virus (RSV), HIV-1, HIV-2, Influenza viruses PR8, H1N1, H6N1, Herpes simplex virus type 1 (HSV-1), HPV-16, HPV-18, Human T-cell leukemia virus type 1 and Japanese encephalitis virus (JEV)	100 μM (HIV-1) and 250 μM (HIV-2)	Inhibits herpes virus replication, inhibits HIV-1 LTR-directed gene expression, and inhibits transcription of HPV-18.	[99], [100], [101]
10.	Tinospora cordifolia (Willd.) Miers (Giloy or Guduchi)	Tinosporin	H1N1 flu, HIV-1, Herpes simplex virus (HSV-1) and Hepatitis-A virus	40–50 μg/mL (Hepatitis-A Virus)	Increases IgG antibodies, inhibit HIV-1 protease enzyme activity.	[102], [103], [104]
11.	Cinnamomum cassia (L.) J. Presl (Cinnamon bark)	E-cinnamaldehyde	Human respiratory syncytial virus (HRSV), H1N1 avian influenza virus subtype H7N3, Human immunodeficiency virus (HIV-1) and Herpes simplex virus type-1	125 μg/mL (H7N3 influenza A virus)	Cleavage of viral RNA and prevents its translation to the viral protein, inhibiting viral attachment to host cells.	[105], [106], [107]
12.	Cordyceps militaris L. (Caterpillar fungus or Kida jari)	Cordycepin	Human immunodeficiency virus type 1 (HIV-1), Influenza virus, Murine leukemia virus, and Epstein-Barr virus (EBV)	125 μM (EBV in SNU719 cells)	HIV-1 reverse transcriptase inhibitor and selectively inhibits influenza viral genome replication.	[27], [108], [109], [110]

Due to the gravity of the situation and worldwide rapid spread of COVID-19, an urgent and complementary effort from researchers is necessary to find therapeutic agents and new preventive methods. Drugs like hydroxychloroquine, lopinavir, ritonavir, arbidol, remdesivir, and favipiravir are currently undergoing clinical trials to test their effectiveness and safety in the treatment of coronavirus disease 2019, and some promising results have been achieved thus far [28].

During COVID-19, spike protein, main protease, human angiotensin-converting enzyme-2, RdRp and furin play vital roles during SARS-CoV-2 infection [29], replication [30] and multiplication [31] in the host body [32]. Therefore, in the present study, these proteins have been targeted with selected medicinal plants to identify potential compounds through in silico methods. In the first phase, virtual screening of compounds was performed using molecular docking and pharmacophore model generation. In the second phase, molecular dynamics (MD) simulation study was performed to further validate the molecular docking results. The binding free energy for these selected protein-ligand complexes was also calculated through the molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method. Structure alignment based on 3D descriptors of potent compounds and machine learning algorithms were also used to validate intermolecular interactions. The information obtained through computational approaches provides promising results that reveal the repurposing potential of bioactive compounds from Ayurvedic plants against SARS-CoV-2.

Materials and methods

Selection of antiviral plants

The candidate plants with antiviral activities were selected based on related literature in Ayurveda with historical roots in the Indian subcontinent and the available literature in PubMed-NCBI and Google (References against each candidate plants have been shown in Table 1). The priority of selection was based on one or more criteria:

Claimed efficacies in treating various lung diseases (asthma, bronchitis, pneumonia and severe acute respiratory syndrome) in Ayurveda and traditional medicine.

Availability of valid experimental proofs (in vitro or in vivo) on lungs protective potentials.

Published antiviral activities against genetically-close human viruses, such as influenza (A and B), respiratory syncytial virus, coronaviruses, rhinoviruses, parainfluenza viruses, adenoviruses, herpes simplex, measles and chickenpox [33].

A systematic review of scientific studies, following guidelines of the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA statement), was carried out [34], [35] on the selected medicinal plants. Furthermore, the PubMed-NCBI database (http://www.ncbi.nlm.nih.gov/pubmed) was systematically searched for experimental trials and related bioassays.

Acquisition of chemical compound information

As shown in Table 1, the potential compounds, which showed strong antiviral activity against genetically-close human viruses, were identified from the PubMed-NCBI literature. The structures of the natural compounds were downloaded in SDF/Mol2 file format embedded with 3D properties from the PubChem Compound [36] database (http://pubchem.ncbi.nlm.nih.gov/), DrugBank (https://www.drugbank.ca/) [37] and ZINC database (https://zinc.docking.org/) [38]. The molecular arrangement and geometry of all the compounds were fully optimized using the semiempirical quantum chemistry method (PM3) [39]. Finally, the optimized structure of all the ligands was exported in Mol2 format and used for molecular docking study. Fig. 1 shows the chemical structure of potent antiviral compounds identified from Ayurvedic antiviral medicinal plants.

Fig. 1

Chemical structures of potent antiviral compounds identified from Ayurvedic plant sources.

Target selection

The SARS-CoV-2 associated target proteins, namely spike protein receptor-binding domain (ID: 6VW1) [40], main protease (6LU7) [41], angiotensin-converting enzyme-2 (ACE2) (ID: 1R42) [42], human furin (ID: 4RYD) [43] and RdRp (ID: 6M71) were used in the present study. The crystal structures of all the target proteins were downloaded from the RCSB Protein Data Bank (PDB) (Source: http://www.rcsb.org/). All the aforesaid target proteins (viral and host) have been reported to intercede SARS-CoV-2 infection, replication, survival and multiplication in the host body [44], [45].

Molecular docking

The molecular docking between potent antiviral compounds and SARS-CoV-2 associated target proteins were studied using Molegro Virtual Docker (MVD, Molexus IVS, Denmark) 2010.4.0 software for Windows [46], [47]. The molecular docking was performed using a 15 Å grid radius; the number of runs: 10; maximum interactions covered: 1500; maximum population size generated: 50; maximum steps followed: 300; neighbor distance factor: 1.00; the maximum number of poses generated: 5 to cover the ligand-binding site of the target proteins structure [48]. As per MVD algorithm, the docking score is presented as an arbitrary unit; the lower the score is the better the binding affinity [48]. The protein-ligand complex and atomic level chemical interactions were further analyzed and visualized by Chimera software [49] and Discovery Studio [50]. The output data of molecular interaction was shown in radar display and scatter plot using Microsoft Excel-2007 program (Supplementary figures). During graph plot (both radar and scatter plot) X-axis was set as the name of ligands and Y-axis was the MVD docking scores (mean of five poses).

Molecular dynamics simulations

Molecular dynamics (MD) simulations have been carried out by using Groningen Machine for Chemical Simulations (GROMACS 5.1.5) software [51] to study the dynamic behavior of the top-ranked protein-ligand complexes. The topology parameters of protein and ligands were generated by CHARMM27 [52] using an all-atom force field in GROMACS and SwissParam [53] respectively. Simulations were carried out with the TIP3P water model in dodecahedron boxes. The charged system was neutralized with NaCl solution. The processed setup was then energy minimized with the steepest descent algorithm by running 50,000 steps to remove bad contacts. Before the simulation run, each system was equilibrated for NVT (constant number of atoms, volume and temperature) using V-rescale thermostat [54] for 500 ps at 300 K temperature and NPT (constant number of atoms, pressure and temperature) was performed at 1.0 bar by Parrinello-Rahman barostat [55] for 1000 ps. The protein backbone was restrained while the solvent and counter ions were allowed to move during the equilibration phase. The LINCS algorithm [56] was used for all bound constraints. Particle mesh Ewald (PME) method was used for long-range electrostatics. During simulation run position restraints were removed. Finally, systems were subjected to 100 ns simulation run each under periodic boundary conditions at 300 K temperature and 1.0 bar pressure.

Binding free energy calculations

Protein-ligand binding affinity was analyzed using the classical simulation method molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) [57]. For current study, the g_mmpbsa tool [58] of GROMACS was employed to calculate the binding affinity of the simulated protein-ligand complexes by taking 50 snapshots from last 50 ns MD trajectories at regular intervals.

The total binding energy (ΔGbinding) of the complex in solvent can be defined as [59]:

Here, Gcomplex is the total energy of the protein-ligand complex and Greceptor and Gligand are individual energies of respective molecules. Further, free energy for each component is defined as [60]:

In the above equation, x denotes the ligand, protein or protein-ligand complex. EMM and Gsolvation are the average molecular mechanics potential energy in the vacuum and free energy of solvation, respectively. EMM can be calculated by using the following equation [61]:where, Ebonded represents bonded interactions composed of the bond, angle, dihedral and improper interactions. Enon-bonded comprises van der Waals (Evdw) and electrostatic (Eelec) interactions. Also, Gsolvation or free energy of solvation is the combination of Gpolar and Gnon-polar contribution and in MM-PBSA method, implicit solvent model was used to calculate free energy of solvation [62]:

The electrostatic component, Gpolar is calculated by solving the Poisson–Boltzmann (PB) equation [63] and non-electrostatic component, Gnon-polar was calculated with solvent accessible surface area (SASA) model equation [60]:Here, γ is a coefficient related to the surface tension of the solvent; b is a fitting parameter.

The final ΔGbind values for protein-ligand complexes were the average values from the 50–100 ns of MD simulation trajectories.

Pharmacophore modeling

Pharmacophore features of top-ranked compounds (higher affinity with target proteins) were determined using the Ligandscout platform, which also reveals possible structural activity relationship (SAR) [64] and optimal molecular interactions with a specific biological target(s). The fully optimized 3D chemical structure of all the compounds in Mol2 format was loaded into Ligandscout software and key pharmacophore features were identified based on 3D descriptors such as H–bond donor, H–bond acceptor, hydrophobic, aromatic, halogen bond donor, positively and negatively ionizable groups [65].

Validation based on structural superimposition

Molecular shape similarity analysis is based on the similarity of small molecules as the index for determination of minimum common structure from a set of compound libraries. Like pharmacophore theory, the molecular shape similarity index can effectively utilize the overall structural features of integrated compounds [66]. All the selected ligands in Mol2 format were zipped and loaded into PharmaGist web server [67]. The current study emphasizes more to obtain potent compound(s) considering the rule that molecules with similar structures may have similar or the same biological activity [68].

Validation based on machine learning algorithms

Bayesian machine learning models (FDA-approved drug screening) from Assay Central platform (https://assaycentral.github.io/#) was used to identify the applicability score of potent compounds that may work against SARS-CoV-2. Assay Central is a tool used for building machine learning models (Bayesian, Random Forest and Deep Neural Networks, etc.) that in turn can be used to filter and score target-specific lead compounds prior to testing [69]. It is a collection of predictive Bayesian and Random Forest models in a self-contained executable for non-experts to evaluate the likelihood of activity in a target of interest. The output includes a prediction score and applicability score (AS) (percentage of molecular fragments present in the model input). The prediction score (Pm) is the summation of fingerprint i and contributions, Ci. Contribution is based on the number of active compounds with the fingerprint, Ai, out of the total compounds with the fingerprint, Ti, and the total number of active compounds in the dataset, R [69].

.

.

Statistical analysis

The experimental results were expressed as mean ± S.D., n = 3. The data analysis was performed by an Analysis of Variance (ANOVA) followed by Tukey's test considering *P ≤ 0.05 as the statically significant values.

Results and discussion

Thus far, there are no specific therapeutic agents for controlling coronavirus infections. Drug repurposing is the most efficient alternative in a situation like a pandemic to develop a new cure quickly [30]. A lot of safety data, doses, pharmacokinetics, and bioassays results are already available for the molecule, which allows for a substantially rapid drug development process. Therefore, the present study was undertaken to determine potential therapeutic agents against COVID-19 from the available Ayurvedic antiviral drug candidates by performing a target-based molecular interaction study followed by molecular dynamic and MM-PBSA binding energy calculation.

The molecular docking method is useful for modeling the interaction between a small molecule and target proteins at the atomic level, so that the behavior of small molecules at the binding site of target proteins and fundamental biochemical processes can be elucidated [46], [47]. The docking process includes two basic steps: determination of the conformation of the ligand, its position and orientation (pose) and assessment of its binding affinity [70]. These two parameters have been successfully determined in the present study during sampling and scoring functions and are discussed in detail in forthcoming sections.

Docking analysis of angiotensin-converting enzyme-2 (ACE2)

Wan et al. [71] demonstrated that residue 394 (glutamine) in the SARS-CoV-2 receptor-binding domain (RBD) can be recognized by critical lysine 31 on the human ACE2 receptor [72]. Moreover, SARS-CoV-2 binds to ACE2 ten times more strongly than other coronaviruses, making SARS-CoV-2 more infectious than others [73]. So the spike protein of SARS-CoV-2 was predicted to have a strong binding affinity to human ACE2. Thus, target angiotensin-converting enzyme-2 (ACE2) receptor with small molecules may be a viable solution that may restrict the coronavirus from entering into the cells [30]. In this context, molecular docking was performed and the result showed that curcumin has strong binding affinity (−88.9) for ACE2 followed by rosmarinic acid (−88.1) and gingerol (−74.3) (Table S1 and Fig. S1). To have a better understanding, the comparative docking scores relative to the center point of all the twenty antiviral compounds with target proteins are represented in radar display (Fig. S1a). Alternatively, a self-explanatory graphical representation in the form of a scatter plot has been shown, where the individual mean docking scores (poses) are presented in a matrix with top ranked ligands (yellow oval area) (Fig. S1b). In another study, it has been reported that inhibition of myocardial fibrosis by curcumin is associated with modulating expression of angiotensin II (Ang II) receptors and angiotensin-converting enzyme 2 (ACE2) in male Sprague Dawley rats [74].

Docking analysis of furin

For achieving a successful fusion of SARS-CoV-2 with host target cells, the spike protein needs to be cleaved by proteases; this step is critical, as it allows the fusion sequences to be exposed and is generally, mediated by the host furin (α protein convertase) [32]. RRAR, a special furin-like cleavage site (FCS) in SARS-CoV-2 spike protein (S) that is absent from other lineage B βCoVs, such as SARS-CoV, has been reported to be responsible for its high infectivity and transmissibility [75]. Furin is found to be expressed in significant concentrations in the lungs [76]. Thus, this enzyme can be used by viruses in the respiratory tract to transform and activate their own surface glycoproteins [75], [77] that ultimately facilitate interaction with the cell surface receptor ACE2 [78]. This makes their role in viral protein processing noteworthy [75]. Furins are also known to regulate avian influenza A virus infection in which HA glycoprotein cleavage is needed to enter into the host cell [77]. Therefore, furin becomes an additional drug target in terms of inhibiting viral fusion with host cells [32]. Molecular docking study revealed that rosmarinic acid (−89.17) bound tightly to furin followed by curcumin (−87.36) and quercetin (−79.60) (Table S1 and Fig. S2). Several studies indicated that rosmarinic acid has anti-inflammatory and antiviral activity in mice infected with the Japanese encephalitis virus [79], [80]. In another study, it has been mentioned that the rosmarinic acid derivative has shown anti- HIV-1 activity by inhibiting HIV-1 integrase [81].

Docking analysis of SARS-CoV-2 spike protein

The large transmembrane spike glycoprotein (type I) of SARS-CoV-2 accounts for notable features; it is heavily-glycosylated, which is thought to mediate viral entry into the host body [29]. Therefore, its inhibition may be associated with decreasing the viral multiplication. The spike glycoprotein has two functional domains, designated as S1 and S2, both of which are necessary for a coronavirus to enter into a cell [29] successfully. The molecular docking study revealed that curcumin has a strong binding affinity (−109.3) with SARS-CoV-2 RBD domain of spike protein (Table S1) followed by gingerol (−89.59) and ribavirin (−85.92) (Fig. S3). Curcumin formed chemical interactions at RBD domain-human ACE2 interface with Tyr505, Ala387, Asp38, Gln493, Glu 35, His34, Glu 37 and Arg393 in the active site. It is worth mentioning that all these amino acids are localized in the interface region of spike glycoprotein and host receptors, ultimately facilitates receptor-mediated endocytosis [82] during primary infection. In a recent study, it has been reported that pan-coronavirus fusion inhibitor (EK1C4), targeting spike protein successfully restricted (IC50 range: 1.3–15.8 nM) viral entry into the host body against different types of coronaviruses such as SARS-CoV, MERS-CoV, SARS-CoV-2, HCoV-OC43 and SARSr-CoVs [83].

Docking analysis of SARS-CoV-2 main protease

Among coronaviruses, an alluring drug target is the ~ 306 amino acids long main protease (Mpro, 3CLpro), to forestall the spread of disease by restraining the cleavage of the viral polyprotein. Mpro is essential for processing the polyprotein that led to the proteolytic activation of the viral functional proteins [31]. An in silico analysis showed that the curcumin (−115.69) binds strongly with Mpro followed by rosmarinic acid (−105.43) and quercetin (−96.04) (Table S1 and Fig. S4). The Mpro active site amino acids such as Thr190, Pro168, Met165, Glu166 and Cys145 are predicted to play a major role during chemical interactions with curcumin.

Docking analysis of SARS-CoV-2 RNA dependent RNA polymerase (RdRp)

The 3D crystal structure of a SARS-CoV-2 RNA Dependent RNA Polymerase (RdRp) has recently been solved by X-ray crystallography. The single-chain of core RdRp relies on virus-encoded Non-structural protein cofactors such as nsp7 and two units of nsp8 for its optimum function. The crystal structure of RdRp was downloaded from PDB (ID 6M71) with a reported resolution of 2.90 Å. RNA dependent RNA polymerase (RdRp) is involved in the replication and transcription of the SARS-CoV-2 genome. It is the cleavage product of the polyproteins 1a and 1ab from ORF1a and ORF1ab [111]. The minimum complex necessary for its proper functioning is completed by attachment of three additional protein peptides (nsp7-nsp8, and one additional nsp8) to the core polymerase which is chain A and contains 851 amino acids residues. The relative binding free energy between remdesivir and ATP was calculated to be −2.80 ± 0.84 kcal/mol, where remdesivir bound much stronger to SARS-CoV-2 RdRp than the natural substrate ATP [112]. The ∼ 100-fold improvement in the Kd from remdesivir over ATP indicates an effective replacement of ATP in blocking of the RdRp preinsertion site. Key residues Asp618, Ser549, and Arg555 are found to be the contributors to the binding affinity of remdesivir [113]. These findings suggest that remdesivir can potentially act as a SARS-CoV-2 RNA-chain terminator, effectively stopping its RNA replication, and hence remdesivir was used in the present study as a reference ligand to compare with other phytochemicals affinity with target proteins. An in silico result showed that the gingerol (−112) binds strongly with RdRp followed by curcumin (−97) and quercetin (−88) (Tables S1, S3 and Fig. S5). The RdRp active site amino acids such as Asn691, Asp623 and Arg624 are predicted to play a major role during chemical interactions with gingerol. Interestingly, the docking score of gingerol (−112) is found to be comparable with reference ligand, remdesivir (−134).

Molecular dynamic simulations

Among the different compounds screened through molecular docking studies, curcumin, gingerol, rosmarinic acid, ribavirin and quercetin showed strong intermolecular interaction with the target proteins with high binding scores (−74 to −115) (Table S1). Therefore, these compounds were selected as the potential inhibitors and used for further analysis through MD simulation. A 100 ns simulation was performed for each system to get deep insights into the stability and more reliable binding mode of the ligand-protein complex. The overall stability of each complex was determined through the backbone root mean square deviation (RMSD) of proteins. Further, these compounds were filtered by the MM-PBSA method, which is a reliable method for binding free energy calculations. In last, the binding mode for the selected potential inhibitors from MD simulations and binding free energy analysis was done.

MD simulation analysis for ACE2

The target protein angiotensin-converting enzyme-2 (ACE2), which was found to have desirable interactions with curcumin, rosmarinic acid and gingerol during molecular docking study were subjected to 100 ns simulation run with respective ligands. After the MD simulation run, the stability of each simulated complex was determined by analyzing the RMSD of protein backbone atoms. The RMSD values for three protein-ligand systems curcumin, rosmarinic acid and gingerol were 0.17, 0.20 and 0.22 nm respectively ( Fig. 2A). The analysis revealed that all three systems obtained a steady-state after convergence and RMSD values were observed within the acceptable value of < 0.3 nm. Further, binding free energy calculations were performed through the MM-PBSA method using the last 50 ns stable RMSD trajectories; a total of 50 snapshots were obtained at regular intervals. The different types of energies which contributed to binding such as van der Waals, electrostatic, polar solvation and SASA energy were calculated; the average values of these energies are depicted in Table S2. The overall average values of the binding energy for curcumin, gingerol and rosmarinic acid were −69.51, −41.62 and −29.18 kJ/mol (Table S2 and Fig. 3A). Our analysis revealed that curcumin displayed stable behavior with ACE2 protein throughout the simulation (100 ns) and displayed the lowest binding free energy and therefore could be considered for further study.

Fig. 2

Analysis of molecular dynamics simulation results. A) Root mean square deviation (RMSD) of backbone atoms of ACE2 in complex with curcumin, gingerol and rosmarinic acid, B) RMSD of backbone atoms of furin in complex with curcumin and quercetin, C) RMSD of backbone atoms of spike protein RBD-ACE2 in complex with gingerol, D) RMSD of backbone atoms of Mpro in complex with quercetin and curcumin.

Fig. 3

Binding free energy analysis from MM-PBSA method over the last 50 ns trajectories. A) The ACE2 complexes with curcumin (orange), gingerol (purple), and rosmarinic acid (green) are represented in different color schemes, B) Furin protein complexes with curcumin (orange) and quercetin (red) are represented in different color schemes, C) Spike RBD-ACE2 in complex with gingerol represented here in purple color, D) The Mpro complexes with quercetin (red) and curcumin (orange) are represented in different color schemes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

MD simulation analysis for furin

The third target protein furin was found to have better interaction with rosmarinic acid, curcumin and quercetin during molecular docking analysis. To further validate these results, all three protein-ligand complexes were simulated for 100 ns independently. The MD simulation revealed that furin protein backbone atoms showed RMSD for bound ligands 0.16 (rosmarinic acid), 0.16 (curcumin) and 0.14 (quercetin) nm respectively. The RMSD values were found within the acceptable value of < 0.3 nm. In last, binding free energy calculation using the MM-PBSA tool was performed similarly as stated for other proteins. The different types of energies which contribute to binding such as van der Waals, electrostatic, polar solvation and SASA energy were calculated and the average values of these are depicted in (Table S2). The analysis revealed that curcumin and quercetin displayed the average binding free energy −89.72 and −34.12 kJ/mol (Fig. 3B and Table S2). Unfortunately, rosmarinic acid displayed a positive binding energy value (1.97 kJ/mol), which confirmed its weak affinity towards furin protein and hence was removed from further analysis. Detailed analyses of curcumin and quercetin have been shown in (Figs. 2B and 3B).

MD simulation analysis for spike protein

The second target protein in the present study, large transmembrane spike glycoprotein (type I) of SARS-CoV-2 in complex with humans ACE2 protein, was found to have better interaction with curcumin, gingerol and ribavirin during molecular docking studies. Each protein-ligand complex was further subjected to 100 ns MD simulation independently. The system's stability was then analyzed by RMSD values of protein backbone atoms. The RMSD values of protein-ligand complexes were 3.61, 3.88 and 2.61 nm respectively. Higher values of RMSD were observed but found to be stable. Subsequently, binding free energy calculations were performed through MM-PBSA method. For calculation last 50 ns trajectories of stable RMSD values were selected and a total of 50 frames were generated at an equal interval. The different types of energies which contribute to binding such as van der Waals, electrostatic, polar solvation and SASA energy were calculated, the average values of these energies are depicted in Table S2. Binding free energy calculations revealed that gingerol had an average binding free energy value of −50.59 kJ/mol, which is better than other tested ligands. Rosmarinic acid displayed a positive value of average binding free energy, which is unfavorable for binding. Unfortunately, the curcumin bound complex was not able to give binding free energy values due to instability of the complex and hence discontinued from further analysis. Based on dynamic characteristics and binding free energy score, gingerol displayed better results than other tested ligands. Therefore, only gingerol was selected for detailed analysis (Figs. 2C and 3C). In a recent study, a reference drug remdesivir was studied against spike RBD-ACE2 protein, and recorded average binding free energy −71.88 kJ/mol from the MM-GBSA method [114]. Our hit compound displayed comparatively high binding energy but still under acceptable limit [115] and hence can be considered for further analysis.

MD simulation analysis for Mpro

The fourth targeted protein in this study is Mpro, which was found to have better interaction with curcumin, rosmarinic acid and quercetin during molecular docking analysis. To further validate these results, a 100 ns simulation run was performed and the stability of the protein-ligand complexes were analyzed by protein backbone atoms RMSD. The protein backbone RMSD was found to be 0.34, 0.27 and 0.22 nm for curcumin, rosmarinic acid and quercetin. Curcumin displayed RMSD values slightly greater than 0.3 nm during the simulation run, but was observed to become stable during the later stage of simulation (Fig. 2D). Binding free energy through MM-PBSA was calculated following a similar protocol as mentioned before for Mpro bound ligands. The different types of energies which contribute to binding such as van der Waals, electrostatic, polar solvation and SASA energy were calculated, the average values of these energies are depicted in Table S2. The average binding free energy for potential inhibitors was −80.05 kJ/mol, −40.41 kJ/mol for quercetin and curcumin respectively (Fig. 3D). Unfortunately, rosmarinic acid displayed a positive average binding free energy value (1.37 kJ/mol), which is not favorable for binding and hence, it has been removed from further analysis. In a recent study, the binding free energy for COVID-19 main protease (PBD: 6LU7) bound inhibitor N3 was found to be −51.07 kJ/mol [115]. In contrast, our hit compound quercetin showed very low binding free energy (−80.05 kJ/mol). This observation strongly suggests that quercetin may bind with Mpro more efficiently than reference inhibitor N3 and therefore considered for further analysis.

MD simulation analysis for RdRp

RdRp, which was found to have better interaction with gingerol, curcumin and quercetin during molecular docking analysis. To further validate these results, a 100 ns simulation run was performed and the stability of the protein-ligand complexes were analyzed by protein backbone atoms RMSD. The protein backbone RMSD was found to be 0.29, 0.24, 0.31 and 0.19 nm for remdesivir, gingerol, curcumin and quercetin ( Fig. 4A). Binding free energy through MM-PBSA was calculated following a similar protocol as mentioned above for other receptors. The different types of energies which contribute to binding such as van der Waals, electrostatic, polar solvation and SASA energy were calculated, the average values of these energies are depicted in Table S2. The average binding free energy of potential hit compounds was −123.22 kJ/mol, −114.4 kJ/mol, −105.5 kJ/mol and −96.3 kJ/mol for remdesivir, gingerol, curcumin and quercetin respectively (Fig. 4B and Table S2). In a present study, the binding free energy for SARS-CoV-2 RNA dependent RNA polymerase (RdRp) with reference inhibitor, remdesivir was found to be −123.2 kJ/mol. Interestingly, our hit compound gingerol showed very close binding free energy with reference inhibitor, remdesivir. This observation strongly suggests that gingerol may bind with RdRp efficiently and therefore considered for further analysis.

Fig. 4

Analysis of molecular dynamics simulation results. A) Root mean square deviation (RMSD) of backbone atoms of RdRp in complex with remdesivir, gingerol, curcumin, and quercetin B) Binding free energy analysis from MM-PBSA method over the last 50 ns trajectories of RdRp. The complexes are represented in different color schemes with remdesivir (olive), gingerol (purple), curcumin (orange), and quercetin (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Binding mode analysis

Molecular dynamics simulation is an attractive approach to study molecular interactions and binding mode analysis [116]. In the present study, binding mode analysis was performed by calculating an average structure from the last 10 ns MD trajectories for all four target proteins complexed with respective ligands. Only those ligands were subjected to binding mode analysis, which displayed an acceptable amount of average binding free energy (Table S3). Interestingly, when ligands were superimposed, it was found that only the best-bound inhibitor from MM-PBSA was able to maintain its interaction inside the active site of the protein, and the rest of the inhibitors lost their interactions or showed binding at different sites. Therefore, the molecules displaying the lowest binding free energy were considered and further analyzed for intermolecular interactions inside the respective target protein’s active site.

Binding mode for ACE2 protein

To assess the binding mode of selected compounds from detailed MD simulations and binding free energy calculations, the average structure was calculated from the last 10 ns MD trajectories and compounds were superimposed. Unfortunately, gingerol was not able to show interactions with desirable active site residues and rosmarinic acid found bound outside the active site which is not important for inhibition of protein. Hence, these two compounds were excluded from further analysis ( Fig. 5A). The selected compound, curcumin was found to make three hydrogen bonds with ACE2 active site residues Ser43, Asn103 and Ser106 (Fig. 5B and Table S3). Our investigation also revealed that curcumin makes two types of non-polar interactions with active site residues. The first one is hydrophobic interactions involving Ser44, Ala46, Ser47, Gly66, Trp69, Ser70, Lys74, Ser77, Glu110. The second type of interaction is pi-alkyl, which was observed with Met62 and Leu73 of ACE2 (Fig. S6A and Table S3). A similar pattern of hydrophobic interactions with other inhibitors of ACE2 has been reported recently [117]. Our analysis showed that curcumin shows better binding affinity and molecular interactions such as hydrogen bonds, hydrophobic interactions, and van der Waals interactions with ACE2 active site residues throughout the simulation. Hence, it can be used for further studies against ACE2 protein.

Fig. 5

Binding mode analysis for human ACE2 and furin proteins. A) ACE2 (cyan color) shown in solid ribbon representation in complex with curcumin (orange color), B) The close depiction of active site of ACE2 with curcumin, here active site residues are shown in cyan sticks C) Furin protein (light pink) shown in solid ribbon representation in complex with curcumin (orange color), D) The close depiction of furin bound curcumin complex, background is shown with protein in line ribbon (pink color). Active site interacting residues are shown in pink sticks. Hydrogen bonds are indicated with green dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Binding mode for furin

Compounds such as curcumin, quercetin and rosmarinic acid were found to have good interaction during molecular docking but during binding free energy rosmarinic acid displayed unfavorable binding free energy values and therefore, excluded from further analysis. Consequently, last 10 ns MD trajectories data was used to calculate average structure for both the ligands, i.e., curcumin and quercetin. During superimposition of structures, it was observed that quercetin moved outside the active site, so forth only curcumin was considered for binding mode analysis (Fig. 5C). During analysis, it was observed that curcumin formed two hydrogen bonds with active site residues Ala292 and Asn295 of furin protein (Fig. 5D and Table S3). Similar residues were also targeted in other studies providing a strong potentiality to curcumin against furin protein[118]. Apart from hydrogen bonds, various hydrophobic interactions were also observed. Curcumin forms van der Waals interaction with residues His194, Ser253, Gly255, Pro256, Trp291, Ser293, Cys305, Thr309, Thr367, Ser368 and pi-pi interaction with Trp254 (Fig. S6B and Table S3).

Binding mode for spike protein

A similar type of analysis was performed for this Spike-ACE2-RBD of SARS-CoV-2. The results obtained after MD simulations and binding free energy calculation revealed that out of three tested ligands, gingerol performs better than other compounds (Table S3). Hence, we further selected gingerol for binding mode analysis. The average structure for the gingerol bound protein was calculated from the last 10 ns of MD simulation trajectories and further analyzed ( Fig. 6A). The analysis revealed that gingerol forms two hydrogen bonds with spike protein RBD domain key residues Ser494 and Gln493 (Fig. 6B and Table S3). A recent study reported Gln493 and Ser494 of spike protein RBD domain interacted with hotspot residues of ACE2 protein, i.e., Lys31 and Lys353 [115], [117]. Similar interactions of gingerol with these residues provide a support to current findings. Gingerol also makes hydrophobic interactions with both the interacting proteins. Hydrophobic interactions with spike RBD were observed with Lys403, Tyr453, Tyr495 and Tyr505; similarly, hydrophobic interactions for ACE2 protein were of two types van der Waals interaction with Glu37, Asp38, Phe356, Ala387, Arg393 and pi-pi interactions with His34, Lys353, Ala386 (Fig. S6C and Table S3). The interaction pattern of gingerol with key residues strongly suggests its inhibitory effect against the Spike-ACE2-RBD complex.

Fig. 6

Binding mode analysis for SARS-CoV-2 spike and Mpro protein. A) ACE2 protein is shown in cyan and RBD domain of SARS-CoV-2 is shown in green color solid ribbon representation. B) The close depiction of protein complex displays gingerol (purple color) in stick representation inside the spike protein RBD active site. Protein in the background is shown as ribbon. Active site interacting residues of ACE2 protein are shown in cyan sticks whereas, RBD domain active site residues are shown in green sticks. C) Mpro (light purple color) shown in solid ribbon representation bound with quercetin, D) The close depiction of protein displays quercetin (red color) in stick representation inside the Mpro active site. Protein in the background is shown as ribbon (light purple color). Active site interacting residues are shown in purple sticks. Hydrogen bonds are indicated with green dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Binding mode for Mpro

Binding free energy calculation revealed that quercetin and curcumin showed the highest binding affinity against Mpro. So quercetin and curcumin binding mode was assessed using average structures calculated from the last 10 ns MD trajectories. Superimposition of Mpro complexes with quercetin and curcumin displayed that only quercetin was found bound inside the active site of Mpro (Fig. 6C). Unfortunately, curcumin was bound outside the active site of protein hence excluded from further analysis. Quercetin forms a hydrogen bond interaction with His164, which is one of the key residues reported for the inhibition of Mpro (Fig. 6D). Quercetin also forms van der Waals interactions with the active site residues Cys44, Tyr54, Cys85, Cys145, Gly174, Thr175, Phe181, Phe185, Val186, Asp187, Gln189, Gln192 and pi-pi interaction His41, Met49, Met165, Arg188 (Fig. S6D and Table S3). A similar type of binding pattern for Mpro protein inhibitors was recently noticed [119], [120]. Hence, our findings confirmed that quercetin binds tightly to Mpro and can inhibit its activity through several molecular interactions.

Binding mode for RdRp

Binding free energy data revealed that gingerol, curcumin and quercetin (Table S2) showed the strong binding affinity against RdRp. Therefore average structure was calculated from the last 10 ns MD trajectories for all mentioned ligands along with remdesivir. It was observed that gingerol formed similar binding mode with RdRp as remdesivir does. Gingerol forms hydrogen bond interactions with Thr591, Lys593 and Asp865 which are one of the key residues reported for the inhibition of RdRp ( Fig. 7, Table S3). Gingerol also forms van der Waals interactions with the active site residues Val588, Gly590, Ser592, Trg598, Met601, Ser759, Cys813, Ser814, Gln815, Arg836 and pi-pi interaction Ile589, Ala688, Leu758 (Fig. S6F and Table S3). A reference inhibitor, remdesivir showed H-bond interactions with key residues Lys593, Cys813, and Asp865 with RdRp (Figs. 7 and S6E). A similar type of binding pattern for RdRp protein inhibitors was recently noticed [112], [113]. Our findings confirmed that gingerol binds tightly to RdRp and can inhibit its activity through several molecular interactions.

Fig. 7

Binding mode analysis for SARS-CoV-2 RdRp protein. A) RdRp protein is shown in peach color solid ribbon representation. B) The close depiction of protein complex displays remdesivir (olive color) in stick representation inside the RdRp active site C) RdRp protein is shown in peach color solid ribbon representation. D) The close depiction of protein complex displays gingerol (purple color) in stick representation inside the RdRp active site. Protein in the background is shown as ribbon. Active site interacting residues of RdRp protein are shown in peach sticks whereas. Hydrogen bonds are indicated with green dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

A pharmacophore model of a bioactive compound is a collection of steric and electronic descriptors that is needed to ensure the stability of the ligand-receptor complex and associated biological response [121]. Pharmacophore features also led to the development of more stable and targeted drugs with no or negligible side effects in the host, as they are supposed to be more accurately targeted [122]. The top-ranked ligands (based on MD result), which showed high binding affinity with all the receptors were further subjected to predict pharmacophore features. The 2D and 3D pharmacophore features of top-ranked compounds were determined by LigandScout 4.1 and shown in Fig. 8. The inspected pharmacophoric features of all the compounds include hydrogen bond donors and acceptors as directed vectors, positive and negative ionizable regions, and lipophilic areas illustrated by spheres. All these above-reported pharmacophore features are essential for binding with receptors. These features would enable repurposing COVID-19 drugs candidates at the same time open up new avenues to design new and more potent SARS-CoV-2 inhibitor(s) as prospective antiviral agents.

Fig. 8

Pharmacophore features (2D and 3D) of top-ranked compounds (based on MD simulation) predicted using Ligandscout software. The pharmacophore color code is yellow for hydrophobic regions, red for hydrogen acceptors and green for hydrogen donors. 2D Pharmacophore features represent HBA: hydrogen bond acceptor, HBD: hydrogen bond donor, AR: aryl, NI: negative ionizable and H: hydrophobic center. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Structural superimposition plays a major role in determining the structural template with a common structure. In the present study, the PharmaGist web server has been used to align all the aforesaid antiviral compounds in 3D orientation, followed by the determination of the pivot molecule(s) [123]. The output file consists of a set of aligned structures with scores based on the common descriptors computed by multiple flexible alignments of the input ligands [67]. After superimposition, the obtained result showed that gingerol, quercetin and curcumin represent pivot molecules and share most of the pharmacophore features of submitted ligands ( Fig. 9) and hence strong intermolecular interactions were observed during docking and MD simulation.

Fig. 9

PharmaGist based ligands superimposition and determination of pivot molecule (key structure). Curcumin, gingerol and quercetin showed structural similarity with high align score and hence may be able to bind SARS-CoV-2 associated target proteins more efficiently than others.

Assay Central (AC) is a stand-alone and user friendly set of predictive Bayesian and Random Forest models for non-experts to determine the likelihood of action of bioactive compounds against the target of interest. After developing the model, each molecule in the selected 'project' receives a relative score known as applicability number that indicates fraction of structural features shared with the training set [69]. Additionally, the prediction score for each compound is also generated that ranges from 0 to 1; a higher score indicates that the compound is more likely to be active at the modeled target. The current threshold value used to consider a compound as active is 0.5 or higher. Furthermore, the applicability score refers to the percentage of fragments of the expected molecule present in the model. The lower applicability score (≤ 0.3) implies that the prediction is not effective because the molecules comprising the model do not have adequate structural similarity to the predicted compound. Such machine learning models could help during prioritization of compounds for testing in vitro and in vivo. The attribute of the ROC or the receiver operator can be interpreted as a percent probability of correctly predicting compound activity against target. Additional metrics included the F1 Score, Cohen's Kappa and Matthew's Correlation Constant (MCC) ranging from 0 to 1 and the value closer to 1 are considered as the most viable model [69]. The Assay Central Bayesian model [69] for coronavirus disease (COVID-19) associated ACE2 receptor has a 5-fold cross-validation receiver operating characteristic (ROC) of 0.84, precision 0.72, recall 1, specificity 0.55, F1-score 0.84, Cohen’s kappa (CK) 0.57, Matthews correlation coefficient (MCC) 0.63 and Random Forest square of the correlation (R2) 0.73 ( Fig. 10A). For ACE2 protein, prediction score (Pm) and applicability score (AS) of potent compound, curcumin was 0.80 and 0.72 respectively. SARS spike protein has a 5-fold cross-validation receiver operating characteristic (ROC) of 0.97, precision 0.68, recall 0.93, specificity 0.92, F1-score 0.79, Cohen’s kappa (CK) 0.74, Matthews correlation coefficient (MCC) 0.76 and Random Forest square of the correlation (R2) 0.62 (Fig. 10B). For spike protein, prediction score (Pm) and applicability score (AS) of potent compound, gingerol was 0.84 and 0.78 respectively. The SARS main protease model has a 5-fold cross-validation ROC of 0.92, precision 0.85, recall 0.98, specificity 0.78, F1-score 0.90, CK 0.77, MCC 0.79 and Random Forest square of the correlation (R2) 0.93 (Fig. 10C). For SARS-CoV-2 main protease prediction score (Pm) and applicability score (AS) of potent compound, quercetin was found to be 0.70 and 0.77 respectively. Moreover, the Assay Central platform lacks data for furin and RdRp therefore; it was not validated by the same method. Hence, based on the prediction (Pm) and applicability score (AS) of Assay Central Bayesian model (Fig. 10A–C) it is revealed that curcumin, gingerol and quercetin has strong affinity for ACE2, spike protein and main protease respectively, interestingly; the result further corroborate the molecular dynamic study. The results were further corroborated by comparing reported IC50 with interaction scores. Interestingly, the compound with low IC50 (reported) has shown strong chemical interactions with target proteins and similar was the case with curcumin, gingerol and quercetin. Thus, based on the machine learning algorithm [69] it can be suggested that aforesaid compounds (curcumin, gingerol and quercetin) possess strong binding affinity against COVID-19 associated target proteins.

Fig. 10

A Assay Central Visualization. A) Machine learning model performance summary of a Bayesian (left) and Random Forest (right) algorithms. Bayesian machine learning models using ECFP6 fingerprint was used for scoring and selecting compounds having potent binding affinity with ACE-2 target protein. B Assay Central Visualization. A) Machine learning model performance summary of a Bayesian (left) and Random Forest (right) algorithms. Bayesian machine learning models using ECFP6 fingerprint was used for scoring and selecting compounds having potent binding affinity with SARS-CoV-2 spike protein. C Assay Central Visualization. A) Machine learning model performance summary of a Bayesian (left) and Random Forest (right) algorithms. Bayesian machine learning models using ECFP6 fingerprint was used for scoring and selecting compounds having potent binding affinity with SARS-CoV-2 main protease.

Though synthetic antiviral drugs are significant in treatment of viral infections, there are innumerable side effects with alterations due to progressive viral complications. Hence, medicinal plants are considered excellent alternative sources of antiviral compounds [124]. Curcuma longa and Zingiber officinale are reported to have strong antioxidant, and antiviral properties and various reports suggested their use in inhibiting SARS-CoV-2 viral replication and pathogenesis [124]. Various biological effects of curcumin have been documented including anti-tumor, anti-hypertensive and antiviral activities [125], [126]. Curcumin mediated antiviral activity was observed against various viruses, including hepatitis viruses, Zika virus, Chikungunya virus, human papillomavirus, human immunodeficiency virus, herpes simplex virus-2 and respiratory influenza virus [127]. Previously, curcumin has been shown to hinder the replication of SARS-coronavirus as well [128] and identified as a constructive inhibitor of protease. There are ample experimental research findings that support the efficiency of curcumin, gingerol and quercetin against the human viruses including SARS-CoV-2 pathogenesis [84], [85], [129]. Gingerol and curcumin extract demonstrated 45% and 40% of protease inhibition at 0.2 mg/mL concentration in an in vitro protease inhibitory assay and the result was further corroborated by pharmacokinetic properties [129]. SARS-CoV-2 genomic resemblance with SARS-coronavirus (> 80%) and Middle East Respiratory Syndrome-coronavirus (50%) collectively suggests the possibility of curcumin efficacy against COVID-19 [130]. Experimental evidences of curcumin in therapeutic and prophylactic respiratory infection management as well as pathophysiology have been identified and support the implementation of curcumin in COVID-19 [131]. Recent experimental evidence from computational approach predicted the direct binding of curcumin with receptors mediating entry of the virus into host cells as well modulation of associated inflammatory events [132], [133], [134]. The reported ability of curcumin to modulate cellular and physiological consequences [127], [128], [135] associated with SARS-CoV-2 infection along with recent in-silico predictions [132], [133], [134] further support its potential. Synergistic antiviral effects of quercetin and ascorbic acid have been effective in preventing viral entry and replication [136]. Experimental and computational evidence augment the efficacy of curcumin, gingerol and quercetin as anti-protease activity [86], [87], [88], [129]. The present study demonstrated the repurposing possibilities of curcumin, gingerol and quercetin based on their binding ability with multiple COVID-19 targets. Furthermore, their inhibitory activity using molecular dynamic simulations, pharmacophore features and Bayesian machine learning models supported their ability to function as SARS-CoV-2 inhibitors.

Since the emergence of COVID-19, SARS-CoV-2 virus has spread enormously across the globe causing a pandemic. Vaccination or the development of herd immunity seems to be the only way to break the chain of virus spread; however, the appearance of similar mutated strain(s) and failure of initially recommended drugs raised concerns in the near future. Therefore, effective therapies are desperately required to alleviate the burden of this pandemic and reduce mortality rates. Various repurposable drugs were proposed to effectively combat COVID-19. In March 2020, emergency use authorization of hydroxychloroquine was approved by the US food and drug administration (FDA), due to its therapeutic immunomodulatory effects [137]. In a study published in May 2020, remdesivir, an inhibitor of viral RNA dependent RNA polymerase (RdRp) that was used in the treatment of ebola in 2004, was widely considered against COVID-19 [138]. A study by [139] showed the significance inhibition of SARS CoV-2 replication by remdesivir and lopinavir in Vero E6 cells with an EC50 of 23.15 and 26.63 μM respectively. Though, lopinavir-ritonavir, reduced the duration of ICU stay in severely affected patients, they also showed severe gastrointestinal side effects [140]. Early interferon therapy was found more promising by reducing viral counts and lessening symptoms but did not reduce fatality [141]. However, the solidarity trial of treatment evaluation published by WHO in October declared that the four drugs remdesivir, hydroxychloroquine, lopinavir/ritonavir and interferon had no significant effect on covid-19 (www.who.int). It is imperative to understand the mandate of extensive research that focuses on repurposing of drugs and the current study has evaluated a wide range of Ayurvedic components with good efficacy to halt viral replication, attachment and entry. Thus, the repurposing of existing drugs is the fastest and cost effective approach to mitigate the COVID-19 burden for the near future. It remains important to identify medications that can be safely used during the early stage of infection to avoid progression to serious illnesses and deaths.

Conclusion

The present study aimed to identify molecules from medicinal plants reported in ancient literature, Ayurveda that may inhibit SARS-CoV-2 by acting on the major target proteins such as spike protein, main protease, angiotensin-converting enzyme-2, furin and RdRp. The obtained results by molecular docking, molecular dynamic simulation, binding free energy calculations, pharmacophore model and structure-based alignment methods showed that curcumin, gingerol and quercetin served as pivot molecules and may work against COVID-19. Furthermore, the binding affinities of all the compounds were also compared with a reference ligand, remdesivir against RdRp target. In addition, we validated molecular docking and molecular dynamic findings using a Bayesian machine learning models using ECFP6 fingerprints and interestingly, prediction score (Pm) and applicability score (AS) of potent compound, gingerol against spike RBD was found to be 0.84 and 0.78 respectively. Furthermore, prediction score (Pm) and applicability score (AS) of curcumin against ACE2 protein was found to be 0.80 and 0.72. Conversely, the prediction score (Pm) and applicability score (AS) of potent compound quercetin against Mpro was turned out to be 0.70 and 0.77. Thus, the results of machine-learning algorithms further corroborate the findings of molecular dynamic simulation and open up new avenues for preclinical investigations of aforesaid compounds in suitable in vivo and in vitro model(s). This study also supports the traditional use of antiviral plants as preventive measures. Hence, it is suggested that the world community should undertake this study further to repurpose clinical studies to test efficacy and safety in the treatment of COVID-19, which is an urgent need of the hour. However, further research is needed to materialize their dosage and safety.

Author contributions

A.K.V. designed the study, performed research, analyzed the data and wrote the paper. V.K. performed molecular dynamic simulation and binding free energy calculations from MM-PBSA and also contributed to writing their parts. S.S., A.S. and I.C. are involved in receiver operating characteristic curve (ROC) validation and analysis, machine learning model validation, data analysis and writing their parts. Professor B.C.G. & K.W.L. contributed in data interpretation and analyzed results. All the authors read, agreed and approved the final manuscript for submission.

Declaration of Competing Interest

There is no conflict of interest for the publication of this article.