Repurposing of the herbal formulations: molecular docking and molecular dynamics simulation studies to validate the efficacy of phytocompounds against SARS-CoV-2 proteins.

Herbal formulations mentioned in traditional medicinal texts were investigated for in silico effect against SARS-COV-2 proteins involved in various functions of a virus such as attachment, entry, replication, transcription, etc. To repurpose and validate polyherbal formulations, molecular docking was performed to study the interactions of more than 150 compounds from various formulations against the SARS-CoV-2 proteins. Molecular dynamics (MD) simulation was performed to evaluate the interaction of top scored ligands with the various receptor proteins. The docking results showed that Liquiritic acid, Liquorice acid, Terchebulin, Glabrolide, Casuarinin, Corilagin, Chebulagic acid, Neochebulinic acid, Daturataturin A, and Taraxerol were effective against SARS-COV-2 proteins with higher binding affinities with different proteins. Results of MD simulations validated the stability of ligands from potent formulations with various receptors of SARS-CoV-2. Binding free energy analysis suggested the favourable interactions of phytocompounds with the recpetors. Besides, in silico comparison of the various formulations determined that Pathyadi kwath, Sanjeevani vati, Yashtimadhu, Tribhuvan Keeratiras, and Septillin were more effective than Samshamni vati, AYUSH-64, and Trikatu. Polyherbal formulations having anti-COVID-19 potential can be used for the treatment with adequate monitoring. New formulations may also be developed for systematic trials based on ranking from these studies.Communicated by Ramaswamy H. Sarma.

Introduction

Coronavirus disease (COVID-19) is a pandemic disease caused by the novel Coronavirus and has been declared as a Global Public Health Emergency by World Health Organization (WHO). As of January 17, 2021, the disease has spread worldwide with 93,217,287 cumulative cases and 2,014,957 deaths, and in India, 10,557,985 cumulative confirmed cases and 1,52,274 cumulative deaths are reported (World Health Organization, 2020). This pandemic situation becomes more challenging to handle in low and middle-income countries because a large proportion of individuals may be at increased risk of infection and face difficulties in receiving quality health services (Dong et al., 2020; Gupta et al., 2020). The governments have prompted various strategies for the prevention and management of disease, but it has become one of the forefronts of health challenges due to the lack of the approved targeted therapy (Guan et al., 2020). Modern molecules such as Remdesivir, Lopinavir/Ritonavir, Favipiravir/Umiferovir, and anti-malarial drugs such as chloroquine, hydroxychloroquine have been reported to have anti-COVID-19 potential, however, reports about their ineffectiveness against COVID-19 is also published (Geleris et al., 2020; Marzolini et al., 2020; Wang et al., 2020). Researchers are actively involved in searching for novel drug molecules alternatives to existing drugs, and pluralistic knowledge systems available worldwide may offer a bright ray of hope (Rastogi et al., 2020). Many of the traditional systems such as traditional Chinese medicine, Kampo, traditional Korean medicine has been practiced in some areas of the world for the treatment of SARS-CoV diseases (Lee et al., 2020; Komuro, 2017; Yang et al., 2020).

In India, The Ministry of Ayurveda, Yoga, and Naturopathy, Unani, Siddha and Homoeopathy (AYUSH), a federal government organization that is actively involved in research and development on Indian traditional medicine (http://ayush.gov.in/) presented a plausible plan of action for Ayurvedic intervention which is not limited to prophylaxis alone but also includes the therapeutic and integrative model of care. They prescribed the use of various formulations such as Chyavanprasha, Brahma Rasayana, Amrit Bhallataka, Swarna prashan, Sanjeevani vati, Chitrakatdi vati, Triphala, Trikatu, Gojihvadi kwath, Kantakari Avaleha, Chitrakadi vati, Vyaghri haritaki, Dashamul kwath, Sitopaladi, Talishadi, Yashtimadhu, Laghu Vasant Malati rasa, Tribhuvan Keerti rasa, Brihata Vata Chintamni rasa, Mrityunjaya rasa, Siddha Makardhvaja, and plants such as Tinospora cordifolia, Zingiber officinale, Curcuma longa, Ocimum sanctum, Glycyrrhiza glabra, Adhatoda vasica, Andrographis paniculata, Swertia chirata, Moringa oleifera, for the treatment of patients with mild COVID-19 to severe COVID-19 symptoms (Rastogi et al., 2020; http://ayush.gov.in/).

Based on these recommendations, we undertook to investigate interactions of few polyherbal formulations with different proteins of SARS-CoV-2, which are involved in initiating, penetration, biosynthesis, maturation, and release of SARS-CoV in host cells. We have selected more than 150 compounds reported from medicinal plants used in these formulations and studied their interactions against target proteins of SARS-CoV-2 using molecular docking. This study identifies potential compounds, plants, and ayurvedic polyherbal formulations that have shown better scores compared to drugs such as remdesivir, hydroxychloroquine, favipiravir. A schematic diagram of the workflow is given in Figure 1.

Figure 1.

Schematic diagram of the workflow.

Materials and methods

Protein preparation

The crystal structures of target proteins were downloaded in pdb (gz) format from the RCSB protein data bank (PDB) (https://www.rcsb.org/). A list of the selected protein targets is given in the supplementary material (Table S1). Protein structures were prepared by removing water molecules and heteroatoms followed by the conversion in .pdb format using BIOVIA Discovery Studio Visualizer v20.1.0.19295 (Dassault Systèmes, San Diego, CA, USA). The protein structures were prepared for molecular docking by adding polar hydrogen atoms and Kollman charges on them. This procedure was followed by the conversion of the pdb structure of protein into pdbqt and grid preparation using AutoDock Tools version 1.5.7 (ADT; Scripps Research Institute, La Jolla, San Diego, USA).

Table 1.

Docking scores of drugs used for the treatment of COVID–19 against SARS-CoV-2 proteins.

SARS-CoV-2 protein targets
	AVG	6LU7	6M03	6Y84	6VSB	6LXT	6LVN	6M3M	6VYO	6VWW	6M17	6VW1	6W01	6NUR	6W4B	6W02	6M71
Remdesivir	–7.3	–7.5	–4.3	–9.0	–6.4	–7.4	–4.8	–8.1	–6.7	–8.3	–8.1	–8.8	–5.8	–8.1	–6.3	–8.4	–8.1
Hydroxychloroquine	–5.8	–5.7	–3.6	–7.3	–5.0	–6.1	–4.0	–7.2	–5.7	–6.7	–6.1	–6.1	–6.6	–6.1	–5.4	–5.4	–6.2
Chloroquine	–5.7	–4.5	–5.5	–7.0	–5.5	–5.7	–3.8	–7.0	–5.3	–6.0	–6.8	–6.2	–6.2	–5.9	–5.4	–4.2	–5.4
Favipiraviar	–5.4	–5.7	–5.3	–5.6	–5.5	–5.7	–3.4	–5.3	–5.6	–5.8	–6.4	–5.3	–5.8	–5.8	–4.6	–4.7	–5.2
Umiferovir	–7.4	–6.4	–7.8	–7.7	–7.6	–7.7	–5.5	–8.2	–6.4	–8.8	–6.5	–8.6	–8.8	–8.3	–7.0	–6.2	–6.9
Lopinavir/Ritonavir	–6.5	–5.2	–4.1	–9.4	–5.4	–6.1	–3.8	–8.1	–5.8	–6.7	–7.7	–8.4	–7.2	–7.3	–6.9	–6.0	–6.6
Azithromycin	–9.4	–5.4	–6.3	–12.2	–8.5	–10.2	–6.5	–10.4	–9.3	–11.9	–11.2	–10.6	–11.0	–10.2	–9.1	–8.9	–9.1
Ivermactin	–8.8	–6.9	–8.4	–9.5	–10.5	–9.2	–5.9	–9.6	–7.7	–9.7	–10.3	–10.9	–10.6	–9.5	–7.6	–6.4	–8.3
Aspirin	–5.6	–5.3	–5.5	–5.8	–6.1	–4.9	–3.8	–5.9	–5.6	–6.6	–5.3	–6.1	–6.7	–5.9	–5.1	–4.7	–5.5
Darunavir	–10.5	–9.2	–9.0	12.1	–10.2	–10.7	–7.3	–12.6	–9.4	–12.7	–11.0	–11.5	–12.1	–11.2	–9.5	–8.7	–10.2
Galidesivir	–6.3	–4.5	–4.8	–7.5	–5.6	–6.3	–4.3	–7.4	–6.2	–7.5	–7.0	–7.8	–7.5	–7.5	–5.2	–5.3	–6.3
Rilpivirine	–1.7	–4.7	–1.0	–2.1	–1.6	–1.4	–1.0	–1.6	–1.3	–1.4	–1.5	–1.3	–1.5	–1.4	–1.2	–1.7	–2.5
Velapatasvir	–3.0	–1.8	–2.6	–3.5	–2.8	–3.2	–2.0	–2.9	–3.0	–3.3	–3.5	–3.6	–3.5	–3.4	–2.5	–2.8	–4.2
Danoprevir	–8.8	–8.1	–5.5	–11.0	–8.3	–8.5	–5.8	–10.1	–8.0	–10.4	–10.2	–10.6	–10.0	–10.0	–7.7	–7.2	–8.8

Ligand preparation

Chemical structures of more than 150 compounds from various polyherbal formulations and 14 FDA-approved drugs were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) in the Spatial Data File (.sdf) format). Chemical structures of ligands in the .sdf format were converted to the .pdb format using BIOVIA Discovery Studio Visualizer. Ligand structures were prepared by adding nonpolar hydrogens, Gasteiger changes, and rotatable bonds and converted in .pdbqt format using AutoDock Tools. A list of formulations is given in supplementary material (Table S2).

Table 2.

Docking scores of top 25 compounds from herbal formulations against SARS-CoV-2 proteins.

Sr. no.	Compounds	Plant name	Common Name	AVG	6LU7	6M03	6Y84	6VSB	6LXT	6LVN	6M3M	6VYO	6VWW	6M17	6VW1	6W01	6NUR	6W4B	6W02	6M71
1.	Liquiritic acid	Glycyrriza glabra	Yashtimadhu	–13.6	–12.3	–12.9	–14.8	–13.2	–13.0	–8.9	–15.1	–14.9	–16.0	–14.4	–12.6	–12.4	–15.0	–14.2	–14.3	–14.1
2.	Terchebulin	Terminalia chebula	Harde	–13.0	–9.8	–13.4	–15.0	–14.1	–12.3	–7.7	–14.7	–14.0	–13.8	–15.5	–14.2	–11.8	–12.6	–12.8	–13.4	–12.7
3.	Liquorice acid	Glycyrriza glabra	Yashtimadhu	–12.7	–11.3	–11.6	–14.7	–12.5	–12.5	–8.4	–14.6	–13.2	–13.5	–15.0	–14.3	–11.8	–13.2	–12.4	–12.7	–12.1
4.	Glabrolide	Glycyrriza glabra	Yashtimadhu	–11.6	–11.8	–13.0	–12.0	–12.9	–11.0	–8.8	–11.0	–10.0	–12.2	–11.4	–10.6	–13.1	–13.2	–11.8	–11.0	–12.4
5.	Casuarinin	Terminalia chebula	Harde	–11.6	–11.0	–11.8	–13.0	–12.0	–11.5	–7.1	–13.2	–12.0	–12.8	–13.0	–13.2	–13.6	–9.9	–9.6	–9.8	–12.3
6.	Corilagin	Phyllanthus emblica	Amla	–11.0	–10.0	–12.0	–11.9	–11.6	–9.2	–7.4	–12.0	–11.2	–13.6	–11.8	–12.0	–11.7	–10.4	–11.5	–9.9	–9.8
7.	Neochebulinic acid	Terminalia chebula	Harde	–10.9	–9.0	–9.7	–12.2	–11.8	–10.4	–7.1	–13.6	–9.8	–12.1	–11.3	–11.0	–10.4	–13.3	–9.6	–11.6	–11.3
8.	Chebulagic acid	Terminalia chebula	Harde	–10.8	–9.8	–10.4	–13.7	–10.2	–10.3	–7.4	–12.2	–10.6	–13.9	–10.5	–10.8	–10.6	–11.1	–10.6	–10.9	–10.4
9.	Daturataturin A	Dhatura metel	Dhaturo	–10.7	–9.3	–8.9	–12.0	–11.4	–10.6	–7.4	–12.2	–9.9	–12.0	–11.5	–11.1	–10.9	–12.6	–9.5	–11.0	–11.3
10.	Taraxasterol	Clitoria ternatea	Gokarni	–10.6	–12.1	–9.8	–10.2	–9.3	–9.7	–8.9	–9.8	–10.3	–10.4	–10.6	–11.4	–11.3	–11.2	–12.3	–10.9	–11.3
11.	Hispaglabridin A	Glycyrriza glabra	Yashtimadhu	–10.5	–8.6	–9.4	–11.2	–9.8	–9.9	–6.5	–11.2	–10.2	–10.9	–11.3	–11.0	–11.6	–12.0	–12.3	–11.1	–11.7
12.	Chebulinic acid	Terminalia chebula	Harde	–10.5	–9.2	–9.9	–12.7	–12.0	–8.6	–6.1	–12.0	–10.4	–11.6	–12.6	–11.4	–9.4	–9.8	–10.0	–11.2	–11.3
13.	Ilekudinoside	Datura metel	Ganthalo Daturo	–10.5	−4.6	–12.0	–8.8	–10.4	–8.9	–7.2	–13.6	–10.4	–11.3	–12.9	–10.5	–12.3	–13.1	–9.1	–10.6	–11.5
14.	Atidin	Aconitum ferox	Vachnag	–10.4	–10.2	–10.7	–11.3	–10.5	–9.3	–6.5	–9.5	–10.6	–12.3	–11.2	–11.4	–11.0	–10.8	–10.7	–10.6	–9.7
15.	Nimbolide	Azadirachta indica	Limdo	–10.4	–9.4	–10.6	–11.2	–10.4	–10.2	–7.2	–12.0	–11.0	–9.8	–11.0	–11.3	–11.0	–10.4	–10.3	–10.1	–9.8
16.	Nimbolinin	Azadirachta indica	Limdo	–10.3	–8.6	–9.3	–12.4	–8.9	–9.8	–6.4	–13.2	–11.6	–11.6	–11.2	–12.1	–11.2	–10.7	–9.9	–9.7	–8.6
17.	Azadirachtin	Azadirachta indica	Limdo	–10.2	–9.4	–10.9	–11.8	–10.5	–10.0	–7.2	–10.4	–10.3	–9.8	–11.0	–10.6	–10.4	–10.3	–10.2	–10.3	–10.2
18.	Palmatoside	Tinospora cordifolia	Galo	–10.1	–8.8	–10.4	–11.6	–11.8	–9.9	–5.1	–12.1	–11.0	–10.8	–9.0	–9.4	–11.0	–10.6	–10.0	–10.0	–9.6
19.	Tinocordioside	Tinospora cordifolia	Galo	–10.0	–8.9	–9.9	–11.0	–10.0	–9.9	–5.4	–11.4	–11.2	–11.1	–10.7	–10.4	–10.9	–11.1	–9.8	–9.6	–9.2
20.	Ergosterol peroxide	Andrographis paniculata	Kariyatu	–10.0	–9.6	–10.0	–11.7	–10.4	–9.9	–7.2	–11.0	–10.3	–10.8	–9.8	–9.6	–9.2	–9.4	–10.1	–10.3	–1–10.6
21.	Picrinine	Alstonia scholaris R.	Saptaparna	–9.9	–9.3	–10.4	–11.0	–10.4	–9.9	–6.9	–10.3	–10.0	–9.1	–10.8	–10.2	–9.6	–9.6	–10.2	–10.4	–10.3
22.	Licuraside	Glycyrriza glabra	Yashtimadhu	–9.9	–9.1	–10.2	–10.3	–10.4	–9.4	–7.1	–10.8	–10.1	–10.3	–10.4	–10.6	–1–10.0	–10.5	–10.4	–9.9	–8.7
23.	Tinosporaside	Tinospora cordifolia	Galo	–9.8	–8.9	–9.7	–9–11.0	–9.7	–10.0	–5.4	–11.5	–9.4	–9.5	–11.0	–11.3	–10.7	–10.4	–10.0	–9.3	–9.6
24.	Nimbanadiol	Azadirachta indica	Limdo	–9.8	–8.7	–10.1	–10.3	–1–10.1	–9.6	–6.1	–10.7	–10.3	–10.2	–9.3	–10.3	–10.6	–10.5	–9.5	–9.7	–10.6
25.	Napelline	Aconitum ferox	Vachnag	–9.8	–10.5	–10.2	–11.7	–10.0	–9.9	–7.1	–10.2	–11.0	–11.8	–11.2	–9.8	–8.8	–8.6	−7.9	–9.0	–8.6

Molecular docking

Target proteins and ligands were docked using AutoDock Vina (The Scripps Research Institute, La Jolla, San Diego, USA) (Trott and Olson, 2010). A three-dimensional grid box to define binding pocket was set using MGL tools. Binding pockets were optimized for each protein. Docking of each ligand was done to each receptor with grid coordinates and grid boxes of certain sizes for each receptor using AutoDock Vina installed in the ToolShed Galaxy server (The Galaxy Team, 2014) was used for molecular docking. Ten binding poses were generated for each ligand and protein. Results were obtained and ranked by the binding affinity (Shukla et al., 2019), which was predicted as negative Gibbs free energy (ΔG) scores (kcal/mol). The docked ligand poses and protein were saved as a single file using PyMOL (DeLano Scientific LLC, Palo Alto, CA, USA and the interaction of proteins and ligands were analyzed using the BIOVIA Discovery Studio Visualizer. Protein-ligand docking poses were analyzed by studying the hydrogen-bond interactions, hydrophobic interactions, and bonding distances. The interaction of each ligand with the target proteins was studied and favorable confirmations were selected.

Molecular dynamics simulation

The protein-ligand complex structure of SARS-CoV-2 proteins and top 4 ranked ligands were subjected to MD simulation using Desmond Molecular Dynamics System, version 6.2, D. E. Shaw Research, New York (Bowers et al., 2006). MD simulation of Remdesivir with various receptors was also performed. A 100 ns simulation was carried out for each complex. Protein-ligand complexes were prepared for simulation by adding hydrogens, filling missing side chains or whole residues. Refinement of prepared structure file was done. In this step, protonation and tautomeric states of variable residues were changed manually, which was followed by the generation of a solvated system for simulation. TIP3P water model was used for the solvation of complexes followed by neutralization by adding ions. After neutralization, energy minimization was performed. OPLS_2005 force field parameters were used for simulations. To perform simulations, the temperature was 300 K, the pressure was 1.01325 bar, and cut off radius was 5 Å. Root mean square deviation (RMSD), root mean square fluctuation (RMSF), the number of hydrogen bonds, and the radius of gyration (Rg) were calculated.

Principal component analysis (PCA)

Protein-ligand confirmations and major global motions upon ligand binding were studied by principle component analysis (PCA). Bio3D, an R package was used to perform PCA. After removal of the rotational and translational movements, the positional covariance matrix of atomic coordinates and its eigenvectors were computed by overlaiding coordinates onto a reference structure followed by the diagonalization on the estimated symmetric matrix by creating the diagonal matrix of eigenvalues. Columns were the eigenvectors corresponding to the direction of motion relative to the initial coordinates in the diagonal matrix (Grant et al., 2006).

Binding free energy calculation

The molecular mechanics energies combined with the generalized Born and surface area continuum solvation (MM/GBSA) method was used to calculate the free energy of the binding of ligands to proteins. The free energy of solvation and MM potential energy were calculated using Prime 4.0 (Schrödinger, LLC, New York, NY). The last 10 ns of the MD trajectories were taken for MM/GBSA (Jacobson et al., 2002, 2004).

Results

Docking analysis generated binding affinities of compounds with SARS-CoV-2 proteins. A higher negative value of binding affinity was considered as the most favorable interaction. The compounds were showing varied docking scores with different target proteins, and the average binding affinity of compounds was calculated to get a holistic picture of the docking results (Related file, Docking results of compounds). Docking of drugs such as lopinavir/ritonavir, remdesivir, favipiravir, chloroquine and hydroxychloroquine, favipiravir, darunavir, ivermectin was also done. Among the drugs, darunavir showed the highest average docking score followed by azithromycin. Docking results of existing drugs are shown in Table 1.

On the basis of the docking results, the top 10 compounds were Liquiritic acid (–13.6), Terchenulin (–13), Liquorice acid (–12.7), Glabrolide (–11.6), Casuarinin (–11.6), Corilagin (–11), Chebulagic acid (–10.8), Neochebulinic acid (–10.9), Daturataturin A (–10.7) and Taraxerol (–10.6). Besides the top 10 compounds, numerous compounds have also shown a good binding affinity with target proteins ranging from –10.5 to –7.5. A list of the top 25 compounds having binding affinity ranging from –13.6 to –9.9 is given in Table 2.

Given the interactions of compounds with the individual proteins, Liquiritic acid showed the highest binding affinities with all the proteins selected for the study. Docking scores of Liquiritic acid ranged from –16 to –8.9, where the highest affinity was observed with NSP15 Endoribonuclease (–16) of SARS-COV-2 followed by the nucleocapsid protein (–15.1). Liquiritic acid showed good interactions with NSP12, main protease, spike glycoprotein, RNA-dependent RNA polymerase, NSP3 also. After Liquiritic acid, Terchebulin also showed good scores with target proteins. It showed the highest score with RBD/ACE2-B0AT1 complex (–15.5), main protease (–15), followed by spike glycoprotein (–14.1).

The ranking of plants based on the binding affinity of compounds is given in Table 3. Plant-wise sorting of docking results showed that the compounds from the plants i.e. Glycyrriza glabra (Yashtimadhu), Terminalia chebula (Harde), Tinospora cordifolia (Galo), Phyllanthus emblica (Amla), Semecarpus anacardium (Bhilamo), Andrographis paniculata (Kariyatu), Dhatura metel (Ganthalo Daturo), Alstonia scholaris (Saptaparna), Aconitum ferox (Vachnag), Azadirachta indica (Neem), Withania somnifera (Ashwgandha) were showing good binding affinities ranging from –13.6 to –9.6. Compounds from Tinospora cordifolia, Piper longum, Azadirachta indica, Phyllanthus emblica, Oroxylum indicum, Aegele marmelos and Ocimum sanctum were also showing good interactions with proteins in terms of binding affinity.

Table 3.

Classifying plants based on higher docking scores of compounds.

Ranking and docking scores	Plants
1–25 (–13.6 to −9.6)	Glycyrrhiza glabra (Yashtimadhu), Tinospora cordifolia (Galo), Phyllanthus emblica (Amla), Semecarpus anacardium (Bhilamo), Terminalia chebula, (Harde), Andrographis paniculata (Kariyatu), Datura metel (Ganthalo Daturo), Alstonia scholaris (Saptaparna), Azadirachta indica (Neem), Aconitum ferox (Vachnag)
26–50 (–9.6 to –8.8)	Piper longum (Lindi Pippar), Oryza sativa (Padad), Oroxylum indicum (Tetu)
51–75 (–8.8 to –8.4)	Withania somnifera(Ashwagandha), Tribulus terrestris (Betha Gokharu), Solanum xanthocarpum (Bethi Bhoy Ringani), Aegele marmelos L. (Bili), Picrorhiza Kurroa (Katuki), Acorus calamus (Godavaj), Curcuma longa (Turmeric), Gmelina arborea L. (Shevan)
76–100 (–8.4 to –8.1)	Desmodium gangeticum (Shalparni), Ginkgo biloba(Maidenhair Tree), Ocimum sanctum(Tulsi), Swertia chirata (Chitara), Caesalphinia crista (Kuberaksha)
101–188 (–8.0 to –6.6)	Terminalia bellirica (Baheda), Rubia cordifolia (Majethi), Uraria picta (Pilo-samervo), Moringa oleifera (Sargavo), Solanum xanthocarpum (Ubhi Bhoin Ringani), Embelia Ribes (Vavding)

Formulation-wise segregation of the docking results showed that Yashtimadhu, Pathyadi kwath, Sanjeevani vati, Septillin, and Tribhvan Keerti ras have plants and phytochemicals in formulations showing higher binding scores in silico. A list of formulations containing the top 10 compounds is given in Table 4.

Table 4.

Formulations containing top 10 compounds based on higher docking scores.

Formulation	Plant	Common name	Compound	Average Docking score
Yashtimadhu	Glycyrriza glabra	Yashtimadhu	Liquiritic acid	–13.6
Pathyadi kwath	Terminalia chebula	Harde	Terchebulin	–13.0
Sanjeevani Vati	Terminalia chebula	Harde	Terchebulin	–13.0
Septillin	Terminalia chebula	Harde	Terchebulin	–13.0
Yashtimadhu	Glycyrriza glabra	Yashtimadhu	Liquorice acid	–12.7
Yashtimadhu	Glycyrriza glabra	Yashtimadhu	Glabrolide	–11.6
Pathyadi kwath	Terminalia chebula	Harde	Casuarinin	–11.6
Sanjeevani Vati	Terminalia chebula	Harde	Casuarinin	–11.6
Septillin	Terminalia chebula	Harde	Casuarinin	–11.6
Pathyadi kwath	Phyllanthus emblica	Amla	Corilagin	–11.0
Sanjeevani Vati	Phyllanthus emblica	Amla	Corilagin	–11.0
Bhoomyamalaki	Phyllanthus niruri	Bhoamli	Corilagin	–11.0
Pathyadi kwath	Terminalia chebula	Harde	Neochebulinic acid	–10.9
Sanjeevani Vati	Terminalia chebula	Harde	Neochebulinic acid	–10.9
Septillin	Terminalia chebula	Harde	Neochebulinic acid	–10.9
Pathyadi kwath	Terminalia chebula	Harde	Chebulagic acid	–10.8
Sanjeevani Vati	Terminalia chebula	Harde	Chebulagic acid	–10.8
Septillin	Terminalia chebula	Harde	Chebulagic acid	–10.8
Tribhuvan Keerti ras	Dhatura metel	Dhaturo	Daturataturin A	–10.7
Yashtimadhu	Glycyrriza glabra	Yashtimadhu	Hispaglabridin A	–10.5

Total eight ligand-protein complexes (including reference ligand Remdesivir) were subjected to MD simulation. A box of water with Na+ and Cl- ions were used to prepare protein-ligand complexes for the simulation. The stability and confirmation of the protein-ligand complexes were assessed by simulating for 100 ns and analyzed by plotting the RMSD, RMSF, a radius of gyration, and H-bonds graph for 100 ns.

RMSD plots of ligands and receptors showed that the ligands were having multiple binding orientations. Figure 1 shows the RMSD values of a protein (left Y-axis), and ligand (right Y-axis). Lig fit Prot and Lig fit Lig show the RMSD of a ligand when the protein-ligand complex is first aligned on the protein backbone of the reference and its reference confirmation, respectively. The RMSD plot of Liquiritic acid-6VSB complex showed the changes in RMSD values larger than 3 Å, which indicates the large conformational change in the protein structure during the simulation. However, the RMSD values of protein and ligand remained intact between 0–15 ns, 20–40 ns, 60–70 ns, and 85–90 ns of the simulation (Figure 2(A)). The RMSD plots of Tinocordioside-6M3M, Daturataturin-6NUR, and Corilagin-6VWW complexes showed that the complexes equilibrated after 80, 60, and 50 ns respectively (Figure 2(B–D)). Noticeably, in all the complexes except liquiritic acid-6VSB complex, stabilized RMSD values indicated that the system had equilibrated after 100 ns of simulation.

Figure 2.

RMSD plots of ligands docked with SARS-CoV2 proteins throughout the MD simulation. (A) Liquiritic acid-6VSB complex, (b) Tinocordioside-6M3M complex, (C) Daturataturin-6NUR complex, (D) Corilagin-6VWW complex.

RMSD plots of reference compound (Remdesivir) with the selected proteins (Figure 3), indicated that the reference compound formed more stable complexes as compared to ligands. However, in the case of Remdesivir-6VSB complex, RMSD values were found to be increasing even after 100 ns simulation.

Figure 3.

RMSD plots of reference ligand docked with SARS-CoV2 proteins throughout the MD simulation. (A) Remdesivir-6VSB complex, (B) Remdesivir-6M3M, (C) Remdesivir-6NUR complex, (D) Remdesivir-6VWW.

The RMSF of the protein-remdesivir complexes and the protein-ligand complexes were calculated and shown in Figure 4. RMSF plots were generated for all the complexes to analyze the residual atomic fluctuations of the protein atoms in the presence of ligands. In RMSF plots, fluctuations were observed during the simulation which was indicated by peaks. It is commonly known that the helices and sheets show less flexibility, whereas loops and turns show higher flexibility (Shukla et al., 2019). Here in the present study, the average RMSF value for Remdesivir-6VSB, 6NUR, 6M3M, and 6VWW complex was 6.57, 6.09, 20.79, and 5.26 Å, respectively, whereas the average RMSF value for liquiritic acid-6VSB, daturataturin-6NUR, tinocordioside-6M3M, and the colilagin-6VWW complex was 12.41, 2.63, 12.73, and 4.62 Å, respectively. The RMSF values indicate that Remdesivr showed lesser fluctuations than daturataturin and corilagin. RMSF plots of Daturataturin-6NUR and Corilagin-6VWW complexes were found to fluctuate more from their mean structure. However, Corilagin-6VWW complex was stable after the 50 ns simulation, whereas initially RMSF value of Liquiritic acid-6VSB complex was decreased after 35 ns and was found to be almost stable between 60 ns simulation.

Figure 4.

RMSF plots of reference ligand and phytocompounds docked with SARS-CoV2 proteins though out the MD simulation. (A) Remdesivir- and Liquiritic acit-6VSB complex, (B) Remdesivir- and Tinocordioside-6M3M complex, (C) Remdesivir- and Daturataturin- 6NUR complex, (D) Remdesivir- and Corilagin- 6VWW complex.

Intermolecular H-bonds are also important factors to be considered to check the stability of protein-ligand complexes during the simulations. The number of H-bonds was also analyzed using simulation trajectories. H-bonds of ligands (reference and phytochemicals) are shown in Figure 5. Among the four ligands (liquiritic acid, corilagin, tinocordioside, daturataturin), corilagin achieved the highest H-bonds (average ∼4). Reference ligand showed the lesser number of H-bonds as compared to phytochemicals. In total, the constant RMSD of protein-ligand complexes during the simulation period was achieved by the contribution of the number of H-bonds. From the results, it can be said that ligand-binding affinity alters with the change of the number of H-bonds.

Figure 5.

Assessment of no. of H-bonds during the MS simulation. (A) Remdesivir and Liquiritic acid, (B) Remdesivir and Tinocordioside, (C) Remdesivir and Daturataturin, (D) Remdesivir and Corilagin.

In addition to H-bonds, hydrophobic and ionic interactions also play an important role in protein-ligand interactions. A schematic of detailed ligand atom interactions with the protein residues is given in Figure 6. Liquiritic acid docked with 6VSB receptor showed the interactions to 3 residues THR302 (47%), LEU303 (48%), LYS304 (45%), whereas remdesivir showed H-bond with ILE312 (65%). In the case of 6M3M receptor, Tinocordioside showed interaction with VAL159 with 32% occupancy, whereas Remdesivir showed interaction with ASN76, SER79 and ILE75. In the case of 6NUR, Daturataturin showed a better interaction as compared to remdesivir.

Figure 6.

Ligand-protein interactions.

The radius of gyration (Rg) values of the protein and ligand complexes were also calculated to analyze the compactness of protein-ligand complexes (Figure 7). Figure 7(A) shows the Rg values of reference compounds. The average Rg values of remdesivir with 6VSB, 6M3M, 6NUR, and 6VWW were 0.52, 0.45, 0.48, and 0.46, respectively. The reference complexes showed less Rg values compared to phytocompound-protein complexes, suggested that it generates more compact complexes compared to the phytocompounds.

Figure 7.

Assessment of radius of gyration during MD simulation.

In the case of 6VSB and liquiritic acid, the Rg value of protein and ligand was found to be 0.90 and 0.88 nm respectively after 40 ns simulation. Protein 6VSB showed fluctuation from ∼1.43 to ∼0.91 nm, whereas ligand liquiritic acid showed fluctuation from ∼1.36 to ∼0.91 nm (Figure 7(B)). In the case of 6M3M and tinocordioside, Rg values were found to be increased after 100 ns simulation from ∼1.0 to ∼1.20 nm (Figure 7(C)), whereas, in the case of 6NUR and daturataturin, the Rg value of protein and ligand was found to be 1.87 and 1.86 nm respectively after 100 ns simulation. Both the protein 6NUR and ligand daturataturin showed fluctuation from ∼1.87 to ∼1.80 nm (Figure 7(D)). Rg values of 6VWW and corilagin were found to fluctuate from ∼1.63 to ∼1.67 nm (Figure 7(E)). However, the fluctuation pattern of the Rg values was the same for ligands and their complex proteins and overall Rg plot analysis showed that ligands Rg curves significantly followed the fluctuation pattern of the respective proteins, indicating the ability of the ligand to stabilize the proteins.

Principle component analysis (PCA)

PCA method was used to study the concerted motions after ligand binding. PCA plot and Porcupine plot of reference complex and liquiritic acid-6VSB complex are given in Figure 8. It shows the statistically meaningful confirmations in both the complexes. The principal motions and the vital motions required for conformational changes were identified in both the complexes. Two different groupings along the PC1 plane were observed, whereas the groupings observed along the PC2 and PC3 planes were not separated. These results indicate the presence of a non-periodic conformational change and periodic global motions. However, the porcupine plot showed the difference between the two complexes more accurately. The direction and magnitude of the motion were presented by the arrows on the protein-ligand complexes. Differences in the motions were observed in the Remdesivir-6VSB complex as compared to Liquiritic acid-6VSB complex. Prominent motions were observed in the Remdesivir-6VSB complex, whereas a large deviation pattern was observed in the Liquiritic acid-6VSB complex.

Figure 8.

Motion analysis of MD trajectory. (A) 6VSB-Luquiritic acid complex (i) PCA plot, (ii) Porcupine plot. (B) 6VSB-Remdesivir complex (i) PCA plot (ii) Porcupine plot.

The principal motions and the vital motions required for conformational changes in other protein-ligand complexes are shown in PCA plots (Figure 9). The graph showed that Remdesivir-protein complexes showed a much stable cluster as compared to all complexes. PCA plots were also generated for the rest of the six complexes. Periods jumps were observed in almost all complexes. 6NUR_Daturataurin and 6NUR_Remdesivir were showing positive distinct correlated motions (blue), may be defined these complexes have stable inhibitors for respected protein.

Figure 9.

PCA plots of MD trajectory. (A) 6NUR-Daturataturin complex, (B) 6NUR-Remdesivir complex, (C) 6VWW-Corilagin complex, (D) 6VWW-Remdesivir complex, (E) 6M3M-Tinocordioside complex, (F) 6M3M-Remdesivir complex.

Binding free energy analysis

All the protein-ligand complexes with reference ligand were subjected for calculation of binding free energy by using MM‐GBSA method. All the complexes including phytocompound-protein complexes and Remdesivir-protein complexes showed a negative binding affinity, indicating favorable interactions. Additionally, other interactions such as electrostatic energy, Van der Waals interactions, nonpolar solvation energy, etc. were also calculated and presented in Table 5. In all the cases, the reference compound showed a better binding affinity as compared to phytocompounds. Other interactions i.e. solvation energy, contributed to the total interaction energy adversely, whereas electrostatic(coloumb), van der Waals energy contributed to the binding energy favorably. However, the results obtained from binding free energy analysis of phytocompounds-protein complexes indicate that the phytocompounds from various formulations can bind the SARS-CoV2 proteins, and they can generate stable complexes with proteins. Figure 10, showing electrostatic interactions which were calculated using the PRIME energy calculation module of Schrodinger, whose surface depicting electrostatic energy ranging from values –30 kcal/mole (blue) to +30 kcal/mole(red). In 6VSB-Remdesivir overall electrostatic energy was contributing more compare to 6VSB_Liquiritic acid. Where ligand seems to contribute electrostatic potential positive in both cases. Minimized structure obtained through MMGBSA, as shown in Figure 10. 6VSB_Liquiritc acid seems to have more hydrogen bond interactions with amino acids GLN957, ASN953, and THR827 compare to remdivisir_6vsb showing interaction with only ILE-312.

Table 5.

The relative free energies (kcal/mol) obtained by Prime MM-GBSA.

Ligand-protein complex	MMGBSA-dG-Binding energy	MMGBSA-dG-Bind in Coulomb	MMGBSA-dG-Bind-Covalent	MMGBSA-dG-Bind-Hbond	MMGBSA-Dg-Bind-Lipo	MMGBSA-dG-Bind-Packing	MMGBSA-dG-Bind-SelfCont	MMGBSA-dG-Bind-Solvation	MMGBSA-dG-Bind-van der Walls
Remdesivir-6VSB	–264.385	–147.499	12.51198	–10.9603	–68.0484	–5.37883	0.382072	127.9998	–173.392
Liquiritic acid-6VSB	–56.429	–28.224	−0.05055	–5.7133	–7.32536	–1.51721	0	24.49377	–38.0924
Remdesivir-6M3M	–42.0687	–16.1653	4.135494	–1.23813	–10.0852	–0.97695	0	20.04884	–37.7875
Tinocordioside-6M3M	–32.0858	–6.06665	3.725588	–1.05783	–12.5111	–1.35488	0	16.52942	–31.3503
Remdesivir-6NUR	–95.4194	11.88252	3.622968	–2.27857	–35.7673	–2.29073	−4.6056	–3.22673	–62.7559
Daturataturin-6NUR	–93.3872	3.363008	–1.33457	–1.46458	–30.2135	–1.20854	0.029084	5.614364	–68.1724
Remdesivir-6VWW	–67.784	–33.5516	5.603524	–4.19151	–12.8151	–0.00081	0	33.09846	–55.9269
Corilagin-6VWW	–47.9162	–28.761	–0.11962	–4.2377	–4.69919	–0.99414	0	29.49081	–38.5953

Figure 10.

Intetaction of refernce compound and liquiritic acid with 6VSB receptor.

Discussion

Despite extensive advancement in modern medicine, the effectiveness of drugs such as remdesivir, hydroxychloroquine, favipiravir was found to be decreased due to the rapid emergence of SARS-CoV-2 drug resistance. These drugs can target cellular functions such as viral replication, terminal glycosylation of ACE2, etc., and affect the viral life cycle. Moreover, Padhi et al. (2021) performed a mutational mapping and provided insight into the functional outcomes of mutations in the remdesivir-binding site in the nsp12 subunit of RdRp. Their study suggested that very few mutations in nsp12 of SARS-COV-2 can lead the resistance against remdesivir. The drug resistance takes place due to the evolutionary pressure on the viral population, and the pressure provides a favorable condition to the subpopulation that has a relatively better fitness with reduced fidelity and emergence of mutation. In the presence of evolutionary pressure, the mutated viral subpopulation will become predominant. In face of such difficulties in the management of diseases, currently, a therapeutic requirement of alternative treatment is high. In this line, natural products can be trusted due to their chemical diversity, biological diversity, and drug-like properties. The traditional medicinal system uses natural products, which have been practiced around the globe for several years. In India, Ayurveda is one of the traditional medicinal systems, in which single or multiple herbs are being used for the treatment. A single herb or polyherbal formulation may contain more than one phytoconstituent which may act synergistically to achieve extra therapeutic effectiveness by maintaining healthy conditions and improving immunity. Simultaneously these compounds may target multiple functions of pathogens at the various stages of the infection. In total, traditional medicine can use the “multi-drugs and multi-targets” mode for the treatment of complex diseases which may be more effective than the individual drug (Padhi et al., 2021; Parasuraman et al., 2014).

As finding solutions for the management of COVID-19 disease is the highest priority in the current pandemic situation, molecular docking studies were planned. SARS-CoV-2 genome encodes for 14 open reading frames and 16 protein replicase-transcriptase that consists of multiple enzymes essential for attachment and entry of the virus into host cells, replication and pathogenicity, virus transcription, infection to human cells as well as transcription process (Cheng, 2007; Dinesh et al., 2020; Gao, 2020; Kirchdoerfer & Ward, 2019; Yan, 2020). In this study, we aimed to evaluate the anti-COVID-19 potential of ayurvedic polyherbal formulations that contain several natural compounds and plant extracts that have been reported for their antiviral potential (Churiyah et al., 2015; Ganjhu, 2015; Saha & Ghosh, 2012; Tan et al., 2013). In the present study, the selection of polyherbal formulations was done based on the advisory of the Ministry of AYUSH for the management of COVID-19 spread in India. Various interventions for the prevention of disease have been proposed by AYUSH which are being used as prophylactic and symptomatic management for disease management (Guidelines for Ayurveda practitioners for covid 19; WHO, 2020).

The objective was to study the interactions of ligands present in polyherbal formulations with the proteins involved in COVID-19 and to analyze the stability of protein-ligand complexes. Compounds from polyherbal/polycomponent formulations were docked into the active site of the SARS-CoV-2 target proteins. Amongst the top five compounds, liquorice acid has been used in the treatment of hepatic disease for more than 40 years in Japan (Li et al., 2014). It has also been reported for its potential to inhibit SARS-coronavirus (SARS-CoV) replication (Hoever, 2005). Another compound casuarinin was also found to be effective against herpes simplex type 2 (HSV-2) in vitro (Cheng et al., 2002).

Plant-wise sorting indicated that compounds from G. glabra, T. chebula, T. cordifolia, P. emblica were found to be effective in silico. In Ayurveda, these plants are extensively used for the treatment of various disorders. Yashtimadhu has been prescribed for respiratory and digestive disorders, and to improve immunity. Glycyrriza glabra reported to have activity against herpes simplex, Varicella zoster, Japanese encephalitis, influenza virus, vesicular stomatitis virus, type A influenza virus (Adam, 1997; Pompei et al., 1979, 1980). In addition to that, authors (Cinatl et al., 2003) demonstrated in vitro effect of compounds present in Glycyrrhiza glabra against two clinical isolates of coronavirus (FFM-1 and FFM-2) from patients with SARS admitted to the clinical center of Frankfurt University, Germany. Likewise, T. chebula has also been reported for various biological properties such as anticancer, anti-inflammatory, antioxidant, anti-protozoal, antimicrobial activity. Badmaev and Nowakowski (2000) reported the protective activity of T. chebula against influenza A virus. Various studies reported the antiviral activity of T. chebula on different viruses (Ajala et al., 2014; Kolla et al., 2017; Oyuntsetseg et al., 2014). Chebulagic acid from Terminalia chebula is reported to have potential as broad-spectrum antivirals for controlling viral infections which engage host cell surface glycosaminoglycans that play a role in initial binding to the host cell (Lin, 2013). Plants such as Tinospora cordifolia, Piper longum, Azadirachta indica, Phyllanthus emblica, Oroxylum indicum, Aegele marmelos, and Ocimum sanctum have been used in various polyherbal formulations of Ayurveda and also reported for their antiviral properties against numerous viruses (Ghoke et al., 2018; Jiang, 2013; Lalita, 2002; Tiwari et al., 2010; Vellingiri et al., 2020).

The results of formulation-wise sorting demonstrate the effectiveness of ayurvedic formulations. Few of these formulations and the plants used in formulations have been also mentioned in the Indian System of Medicine such as Siddha and Unani. Overall, the study demonstrates the efficacy of polyherbal formulations using modern tools.

In the context of epidemics and pandemics, complementary and alternative medicine recommends preventive and prophylactic measures to enhance immunity. The traditional medicinal systems recognize the effect of ecological and environmental conditions on the state of health of humans and make a great effort on the prevention of diseases.

Few plants such as Glycyrriza glabra (Yashtimadhu), Aegele marmelos (Bilva), Terminalia chebula (Harde), Tinospora cordifolia (Guduchi), etc. also reported having anti-inflammatory activity (Bag, 2013; Nirmala & Selvaraj, 2011; Rajaram, 2018; Upadhyay et al., 2010). In severe COVID-19 patients, the uncontrolled release of pro-inflammatory cytokines is very common and the progression of the disease results in the loss of immune regulation due to exacerbation of the inflammatory components (Garcia, 2020). Zhang et al. (2020) shared their experience about the anti-inflammation treatment of patients with severe COVID-19 and described the importance of anti-inflammatory drugs for the treatment in patients suffering from symptoms including multiple organ failure and acute respiratory distress syndrome (ARDS).

In the current outbreak of COVID-19, one other major concern is the D-dimer production. Few studies revealed the increase of D-dimer levels in severe COVID-19 patients (Tang et al., 2020; Zhou et al., 2020). Wichmann et al. (2020) performed autopsies on 12 patients with COVID-19. They showed deep vein thrombosis and pulmonary embolism (PE) in seven and four patients respectively. In this context, Ayurveda has also described the management of D-dimer formation. Herbs of the formulations such as Yashtimadhu, Pathyadi Kwath, Sanjeevani Vati, Samshamni Vati were reported to have antithrombotic activity. Additionally, compounds such as Glycyrrhizin, Corilagin, Withanoferin A also reported for antithrombotic activity (Ku & Bae, 2014; Lugun, 2018; Mendes-Silva, 2003; Ri & Ho Ju, 2018; Saleem, 2019).

Molecular dynamics simulation results indicate that RMSD values of proteins and ligands followed a similar pattern except for liquiritic acid-6VSB complex. Overall, protein-ligand complexes were stable and RMSD fluctuations were observed due to the small size of ligands and partial binding pocket occupancy. The number of H-bonds indicated the higher stability of protein-ligand binding conformations (Londhe et al., 2019). Peele et al. (2020) carried out 20 ns MD simulations on the main protease and inhibitors complex. Their study validated the stability of lopinavir, amodiquine, and theaflavin digallate in the protein binding pocket as potent binders. Similar types of studies have been reported about the screening of natural compounds having anti-COVID potential from Traditional Chinese Medicine (Romeo et al., 2020; Selvaraj et al., 2020).

Overall, a collective approach of docking and MD simulation is beneficial for repurposing ayurvedic formulations for the management of the disease. Docking results and their correlation with in vitro studies reported in the literature generate reliable evidence to use these ayurvedic formulations for the management of COVID-19. Ancient literature such as Siddha, Unani, Ayurveda has the potential to prevent and treat COVID-19 which can be further tested. According to guidelines published by AYUSH, these formulations can be used for the management of the disease. This study provides evidence of the interaction of the natural compounds with the proteins involved in the disease using computational tools.

According to the present study, the formulations tested in this study had significantly higher binding efficacy against their SARS-CoV-2 targets. In silico results reveal that these formulations may be effective inhibitors of SARS-CoV-2 through their binding to the spike glycoprotein, RNA-dependent RNA polymerase, protease, which can be further studied in vitro. Interactions of compounds with the target proteins suggest that these compounds will be efficacious in preventing both viral attachment and replication, as well as these formulations, can be directly used for patients having mild to moderate symptoms of COVID-19. SARS-CoV-2 is reported to have a high predilection in the pharyngeal epithelial cells. As the deliverability of these extracts to the pharyngeal regions is feasible by suitable oral formulations, they will be useful in the medical management of SARS-CoV-2 infections.

Click here for additional data file.