Molecular docking analysis of selected phytochemicals against SARS-CoV-2 Mpro receptor.

Presently world is on a war with the novel coronavirus and with no immediate treatments available the scourge caused by the SARS-CoV-2 is increasing day by day. A lot of researches are going on for the potential drug candidate that could help the healthcare system in this fight. Plants are a natural data bank of bioactive compounds. Many phytochemicals are being studied for various ailments including cancer, bacterial and viral infections, etc. The present study aims to screen 38 bioactive compounds from 5 selected plants viz., Azadirachta indica, Curcuma longa, Zingiber officinale, Ocimum basilicum and Panax ginseng against SARS-CoV-2. Lipinski's rule was taken as the foundation for initial screening. Shortlisted compounds were subjected to molecular docking study with Mpro receptor present in SARS-CoV-2. The study identified that gedunin, epoxyazadiradione, nimbin and ginsenosides have potential to inhibit Mpro activity and their binding energies are - 9.51 kcal/mol, - 8.47 kcal/mol, - 8.66 kcal/mol and - 9.63 kcal/mol respectively. Based on bioavailability radar studies gedunin and epoxyazadiradione are the two most potent compounds which are used for molecular dynamics simulation studies. Molecular dynamics studies showed that gedunin is more potent than epoxyazadiradione. To find the effectiveness and to propose the exact mechanism, in-vitro studies can be further performed on gedunin. © Society for Plant Research 2020.

Introduction

As per World Health Organization (WHO), on August 8, 2020, there were over 18.9 million confirmed cases of COVID-19 worldwide with the total death toll of over 709 thousand. Following the USA and Brazil, India has the third-highest confirmed COVID-19 cases with toll going over 1.7 million cases (https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports/). In December 2019, an outbreak of pneumonia cases in Wuhan, Hubei Province, PRC (People’s republic of China) was recorded. The unknown pathogen entity was soon identified as a novel coronavirus. Therefore, pneumonia caused by the virus was termed as Novel Coronavirus Infected Pneumonia (NCIP). Initially, the virus was called 2019-nCoV, which was later changed to SARS-CoV-2 by the International Committee on the Taxonomy of Viruses. The taxonomy was soon followed by the declaration of a pandemic on March 11, 2020 by the WHO (Huang et al. 2020).

Belonging to the family of Coronaviridae, coronaviruses are enveloped viruses having non-segmented RNA as their genetic material. The 2019-nCoV was isolated and identified from bronchoalveolar lavage fluid samples from a patient in Wuhan district on January 3, 2020. Full genome sequencing and phylogenetic analysis of the virus reported that 2019-nCoV is a distinct clade of beta-coronaviruses and related to SARS-CoV and MERS-CoV (Wang et al. 2020). It was also reported that 2019-nCoV shares 89.1% similarity with SARS-CoV (Wu et al. 2020). Hence, it was renamed to SARS-CoV-2. The structure of Mpro from SARS-CoV-2 was given by Jin et al. (2020). Mpro (6LU7) is an inclusive and important part of the viral replication machinery and hence, making it a potential target receptor to combat against COVID-19 pandemic (Jin et al. 2020).

The most common symptoms experienced by the patient on the onset of disease are cough, fever, and myalgia or fatigue while less common symptoms are headache, sputum production, haemoptysis, and diarrhoea (Huang et al. 2020). Organ dysfunction such as acute respiratory distress syndrome (ARDS), shock, acute cardiac injury, or even death can occur in severe cases (Wang et al. 2020). The SARS-CoV-2 virus has also been shown to disrupt the normal immune responses of the body which further prompts an uncontrolled inflammatory response in severe cases. Hence, the patients with severe cases manifest lymphopenia, lymphocyte activation and dysfunction, granulocyte and monocyte abnormalities, and an increase in immunoglobin G (IgG) (Yang et al. 2020).

The statistics shows that despite of the measures taken by different countries the daily number of cases are still rising. To control the spread there is an urgent need for a novel or repurposed drug to combat the pandemic. Bioinformatics significantly reduces the time and monetary costs that go into designing experiments, their execution, and laboratory trials. Molecular docking applies an innovative method by combining physiochemical principles with the complex scientific calculation algorithms. Docking assists in the drug development based on the characteristic interaction between the receptor and ligand/drug molecule. It helps industry to focus on the compounds with high potential against the particular target (Gao et al. 2020).

Polysaccharides derived from leaves of Azadirachta indica show anti-viral properties against poliovirus, anti-bovine herpes virus type 1, and duck plague virus (Kumar and Navaratnam 2013). Rhizome extracted from Curcuma longa has inhibitory activity against influenza A neuraminidases (Dao et al. 2012) and H5N1 avian influenza virus (Sornpet et al. 2017). Fresh Zingiber officinale was effective against the human respiratory syncytial virus (Chang et al. 2013) and chikungunya virus (Kaushik et al. 2020). Ethanolic extracts and purified compounds from Ocimum basilicum exhibit a wide variety of anti-viral activities against hepatitis B antigen, HSV-1, and ADV-8 to name a few. Panax ginseng has been proven to inhibit the replication of human gammaherpesviruses (Kang et al. 2017). Therefore, in this study, we have selected these five plants i.e. Azadirachta india (Indian lilac/Neem), Curcuma longa (turmeric), Zingiber officinale (Ginger), Ocimum basilicum (Basil) and Panax ginseng (Asian ginseng) to find a potential therapeutic natural compound against Mpro receptor of SARS-CoV-2. The results obtained would help in developing an effective treatment for COVID-19 from plant source.

Material and methods

Protein/macromolecule

Main protease of SARS-CoV-2 was used in this study. The 3-dimensional structure was extracted in PDB format from the RCSB PDB data repository. PDB id given to the structure is 6LU7 and primarily structure is a homodimer having two A chains composed of 306 amino acids and N3 molecule acting as its inhibitor (Fig. 1).

Fig. 1

3D structure of SARS-CoV-2 Mpro with native inhibitor N3

Ligands

A total of 38 bioactive compounds from five different plants including Azadirachta indica, Curcuma longa, Zingiber officinale, Ocimum basilicum and Panax ginseng were selected as ligands and structures were obtained from PubChem databank in.sdf format (Table 1). For the docking purpose, all the ligands were converted into.pdb file format using Biovia Discovery Studio Visualizer.

Table 1

ADME analysis of 38 phytochemicals

S. No	Compound name	PubChem ID	Compound structure	Analysis
1.	Nimbin	108058		Molecular weight (< 500 Da)	540.6
				Lipophilicity (LogP < 5)	2.28
				H bond donor (< 5)	0
				H bond acceptor (< 10)	9
				Violations	1
2.	Nimbolide	12313376		Molecular weight (< 500 Da)	466.5
				Lipophilicity (LogP < 5)	2.2
				H bond donor (< 5)	0
				H bond acceptor (< 10)	7
				Violations	0
3.	Azadirachtin	5281303		Molecular weight (< 500 Da)	720.21
				Lipophilicity (LogP < 5)	1.09
				H bond donor (< 5)	3
				H bond acceptor (< 10)	16
				Violations	2
4.	Salannin	6,437,066		Molecular weight (< 500 Da)	596.71
				Lipophilicity (LogP < 5)	3.93
				H bond donor (< 5)	0
				H bond acceptor (< 10)	9
				Violations	1
5.	Gedunin	12004512		Molecular weight (< 500 Da)	482.52
				Lipophilicity (LogP < 5)	4.22
				H bond donor (< 5)	0
				H bond acceptor (< 10)	7
				Violations	0
6.	Epoxyazadiradione	49863985		Molecular weight (< 500 Da)	466.57
				Lipophilicity (LogP < 5)	4.24
				H bond donor (< 5)	0
				H bond acceptor (< 10)	6
				Violations	0
7.	Curcumin	969516		Molecular weight (< 500 Da)	368.4
				Lipophilicity (LogP < 5)	3.2
				H bond donor (< 5)	2
				H bond acceptor (< 10)	6
				Violations	0
8.	Oleoresin	11979920		Molecular weight (< 500 Da)	1015.1
				Lipophilicity (LogP < 5)	9.69
				H bond donor (< 5)	6
				H bond acceptor (< 10)	15
				Violations	4
9.	Turmerone	14367555		Molecular weight (< 500 Da)	218.33
				Lipophilicity (LogP < 5)	3.3
				H bond donor (< 5)	0
				H bond acceptor (< 10)	1
				Violations	0
10.	Zingiberene	92776		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	5.2
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	1
11.	Curlone	196216		Molecular weight (< 500 Da)	218.33
				Lipophilicity (LogP < 5)	4
				H bond donor (< 5)	0
				H bond acceptor (< 10)	1
				Violations	0
12.	α-Curcumene	92139		Molecular weight (< 500 Da)	202.33
				Lipophilicity (LogP < 5)	5.4
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	1
13.	α-Santalene	94164		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	5.2
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	1
14.	6-Gingerol	442793		Molecular weight (< 500 Da)	294.4
				Lipophilicity (LogP < 5)	2.5
				H bond donor (< 5)	2
				H bond acceptor (< 10)	4
				Violations	0
15.	Zingerone	31211		Molecular weight (< 500 Da)	194.23
				Lipophilicity (LogP < 5)	0.8
				H bond donor (< 5)	1
				H bond acceptor (< 10)	3
				Violations	0
16.	Gingerenone-A	5281775		Molecular weight (< 500 Da)	356.4
				Lipophilicity (LogP < 5)	3.7
				H bond donor (< 5)	2
				H bond acceptor (< 10)	5
				Violations	0
17.	6-Shogaol	5,281,794		Molecular weight (< 500 Da)	276.4
				Lipophilicity (LogP < 5)	3.7
				H bond donor (< 5)	1
				H bond acceptor (< 10)	3
				Violations	0
18.	6-Dehydrogingerdione	22321203		Molecular weight (< 500 Da)	290.4
				Lipophilicity (LogP < 5)	4.2
				H bond donor (< 5)	2
				H bond acceptor (< 10)	4
				Violations	0
19.	β-Bisabolene	10104370		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	5.2
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	1
20.	α-Farnesene	5,281,516		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	6.1
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	1
21.	Apigenin	5280443		Molecular weight (< 500 Da)	270.24
				Lipophilicity (LogP < 5)	1.7
				H bond donor (< 5)	3
				H bond acceptor (< 10)	5
				Violations	0
22.	Linalool	6549		Molecular weight (< 500 Da)	154.25
				Lipophilicity (LogP < 5)	2.7
				H bond donor (< 5)	1
				H bond acceptor (< 10)	1
				Violations	0
23.	Eugenol	3314		Molecular weight (< 500 Da)	164.2
				Lipophilicity (LogP < 5)	2
				H bond donor (< 5)	1
				H bond acceptor (< 10)	2
				Violations	0
24.	β-Elemene	6918391		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	6.1
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	1
25.	α-Bergamotene	86608		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	4.8
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	0
26.	α-Guaiene	5317844		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	4.6
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	0
27.	Germacrene D	5317570		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	4.7
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	0
28.	Cubenol	519857		Molecular weight (< 500 Da)	222.37
				Lipophilicity (LogP < 5)	3.7
				H bond donor (< 5)	1
				H bond acceptor (< 10)	1
				Violations	0
29.	τ-Cadinol	160799		Molecular weight (< 500 Da)	222.37
				Lipophilicity (LogP < 5)	3.3
				H bond donor (< 5)	1
				H bond acceptor (< 10)	1
				Violations	0
30.	α-Gurjunene	15560276		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	4.1
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	0
31.	δ-Cadinene	441005		Molecular weight (< 500 Da)	204.35
				Lipophilicity (LogP < 5)	3.8
				H bond donor (< 5)	0
				H bond acceptor (< 10)	0
				Violations	0
32.	Ginsenosides	3086007		Molecular weight (< 500 Da)	444.7
				Lipophilicity (LogP < 5)	8.5
				H bond donor (< 5)	2
				H bond acceptor (< 10)	2
				Violations	1
33.	Epigallocatechin gallate	65,064		Molecular weight (< 500 Da)	458.4
				Lipophilicity (LogP < 5)	1.2
				H bond donor (< 5)	8
				H bond acceptor (< 10)	11
				Violations	2
34.	Theaflavin digallate	135403795		Molecular weight (< 500 Da)	868.7
				Lipophilicity (LogP < 5)	4.7
				H bond donor (< 5)	13
				H bond acceptor (< 10)	20
				Violations	3
35.	Genistein	5280961		Molecular weight (< 500 Da)	270.24
				Lipophilicity (LogP < 5)	2.7
				H bond donor (< 5)	3
				H bond acceptor (< 10)	5
				Violations	0
36.	Hesperidin	10621		Molecular weight (< 500 Da)	610.6
				Lipophilicity (LogP < 5)	− 1.1
				H bond donor (< 5)	8
				H bond acceptor (< 10)	15
				Violations	3
37.	Neohesperidin	442439		Molecular weight (< 500 Da)	610.6
				Lipophilicity (LogP < 5)	− 0.5
				H bond donor (< 5)	8
				H bond acceptor (< 10)	15
				Violations	3
38.	Diosmin	5281613		Molecular weight (< 500 Da)	608.5
				Lipophilicity (LogP < 5)	− 0.8
				H bond donor (< 5)	8
				H bond acceptor (< 10)	15
				Violations	3

ADME analysis

For the initial screening purposes, a web-based tool named SwissADME (https://www.swissadme.ch/) was used to eliminate a few compounds according to Lipinski’s rule of five parameters. For a compound to qualify as ligand it should have < 500 Da molecular weight, a high lipophilicity i.e. value of Log P being less than 5, hydrogen bond acceptors being less than 10 and H-bond donors less than 5. Any compound with 2 or more violations was ruled out for further study (Lipinski 2004).

Molecular docking

To obtain protein–ligand docked complex Autodock 4.2 was utilized. The downloaded structure of 6LU7 and each ligand was optimized prior to docking. From the protein 3D structure, water molecules and the inhibitor N3 molecule were removed. Addition of polar hydrogen bonds, Kollman charges and Gasteiger charges summed up the protein and ligand optimization. A grid box of 60 × 60 × 60 was prepared around the binding site of the protein with 0.375 Å spacing. Genetic algorithm was set as the search parameter and output was handled in Lamarckian GA run and docking log file (DLG) were obtained for further analysis of binding energy. The analysis of DLG file revealed a total of 10 conformations for each ligand. The conformation with highest negative binding energy was selected and docked complex was converted to a 2D structure to examine the interactions formed at binding site of 6LU7 with ligand.

Bioavailability radar

Drug-likeliness of drug candidates with binding energy less than the control one was analyzed in a comprehensive way taking 6 physiochemical properties into consideration and forming a bioavailability radar using the SwissADME tool (https://www.swissadme.ch/). Six parameters considered: solubility, size, polarity, lipophilicity, flexibility and saturation. The pink shaded region defines the optimal values of the 6 parameters and deviation from these parameters on a large scale is suggestive of the ligand not being orally bioavailable (Diana et al. 2017).

Molecular dynamics

The receptor-ligand complex obtained after docking were subjected to molecular dynamic simulation using the Desmond-Maestro module 2020. The software provides high performance algorithms as its default settings and all the default setting were used to obtain high speed and precise results. The docked complex was first submerged in a TIP3P water model in an orthorhombic shape. The whole system was neutralized by adding 3 sodium ions at 0.15 M concentration. All the atoms were aligned using optimized potentials for liquid simulations-AA (OPLS-AA) 2005 force field. NVT was used as the ensemble class with SHAKE/RATTLE algorithm to limit the moment of the covelntly bonded atoms. The simulation was started at 300 K and 1 bar pressure for 50 ns. RESPA integrator was utilized to combine all the parameters of the dynamic simulation. A trajectory for 50 ns was set to show in 1000 frames to analyze the dynamic nature of the interaction and component stability.

Result and discussion

Lipinski’s rule of five was applied to estimate the drug likeliness of the all selected 38 candidates. This comparative method helps us to rule out few compounds according to their physiochemical properties. Compounds violating two or more parameters were out listed and rest of the compounds were considered to be ligands for the docking study. Out of 38 phytochemicals, 7 compounds: azadirachtin, oleoresin, epigallocatechin gallate, theaflavin digallate, hesperidin, neohesperidin and diosmin violated more than 1 parameter henceforth, remaining 32 compounds were subjected to docking studies (Table 1).

All the filtered ligands from the ADME analysis were subjected to molecular docking analysis. Molecular docking is an essential computational tool in the drug discovery domain. It is done to further select the potential compounds and study the bond formation in the protein–ligand complex at the binding site. Figure 1a represents all 10 residues namely: THR24, THR26, ASN142, CYS145, PHE140, HIS163, HIS164, GLY143, GLU166, HIS172 which are present in at the active site of the Mpro protein. N3 (native inhibitor) was taken as a control and comparative study of the docking results of all 31 ligands (Table 2) with the control revealed that four compounds having better binding energy as compared with the binding energy of N3 (-8.15 kcal/mol) (Fig. 6).

Table 2

Molecular docking results of 31 ligands with 6LU7

S. No	Ligands	Binding energy (ΔG) (kcal/mol)	Ligand efficiency	Inhibition constant (µM)	Intermolecular energy	Vdw H-bond desolvation
1.	Nimbin	− 8.66	− 0.22	0.44923	− 11.05	− 11.05
2.	Nimbolide	− 8.09	− 0.24	1.17	− 9.29	− 9.35
3.	Salannin	− 7.45	− 0.17	3.49	− 10.13	− 9.93
4.	Gedunin	− 9.51	− 0.27	0.10606	− 10.41	− 10.31
5.	Epoxyazadiradione	− 8.47	− 0.25	0.61357	− 9.37	− 9.31
6.	Curcumin	− 7.57	− 0.28	2.82	− 10.55	− 10.4
7.	Turmerone	− 6.43	− 0.4	19.4	− 7.62	− 7.57
8.	Zingiberene	− 6.8	− 0.45	10.43	− 7.99	− 8.00
9.	Curlone	− 7.05	− 0.44	6.79	− 8.24	− 8.19
10.	α-Curcumene	− 6.42	− 0.43	19.7	− 7.61	− 7.61
11.	α-Santalene	− 6.62	− 0.44	14.09	− 7.51	− 7.51
12.	6-Gingerol	− 5.65	− 0.27	72.02	− 9.23	− 9.12
13.	Zingerone	− 5.28	− 0.38	135.49	− 6.77	− 6.67
14.	Gingerenone-A	− 7.00	− 0.27	7.36	− 10.28	− 10.10
15.	6-Shogaol	− 6.61	− 0.33	14.33	− 9.59	− 9.43
16.	6-Dehydrogingerdione	− 6.97	− 0.33	7.79	− 9.95	− 9.61
17.	β-Bisabolene	− 4.98	− 0.33	222.12	− 6.18	− 6.16
18.	α-Farnesene	− 6.37	− 0.42	21.3	− 8.16	− 817
19.	Apigenin	− 7.36	− 0.38	2.88	− 8.75	− 8.61
20.	Linalool	− 5.15	− 0.47	167.82	− 6.64	− 6.6
21.	Eugenol	− 4.96	− 0.41	230.21	− 6.16	− 6.13
22.	β-Elemene	− 6.49	− 0.43	17.57	− 7.38	− 7.38
23.	α-Bergamotene	− 7.27	− 0.48	4.71	− 8.16	− 8.17
24.	α-Guaiene	− 6.88	− 0.46	9.09	− 7.18	− 7.18
25.	Germacrene D	− 6.42	− 0.43	19.7	− 6.72	− 6.71
26.	Cubenol	− 6.99	− 0.44	7.57	− 7.58	− 7.49
27.	τ-Cadinol	− 7.24	− 0.45	4.97	− 7.83	− 7.82
28.	α-Gurjunene	− 6.59	− 0.44	14.78	− 6.59	− 6.59
29.	δ-Cadinene	− 6.63	− 0.44	13.71	− 6.93	− 6.94
30.	Ginsenosides	− 9.63	− 0.3	0.08745	− 11.42	− 11.28
31.	Genistein	− 7.02	− 0.35	7.14	− 8.21	− 8.1
32.	N3 (Control)	− 8.13	− 0.17	1.07	− 13.54	− 13.34

Fig. 6

Molecular docking results of 31 ligands with Mpro of SARS-CoV-2

Among the 10 conformation of Nimbin, − 8.66 kcal/mol was the least binding energy obtained. Five different types of interaction were observed including van der waals, H-bond, alkyl, pi-alkyl and carbon hydrogen bond (Fig. 2). His163 and Cys145 forming the conventional H-bond while Pro168 and Met165 were engaged with a pi-alkyl and alkyl bond respectively. Glu166 and Leu167 were interacting with the ligand using carbon hydrogen bond and remaining residues weakly interact with the ligand via van der waals bond formation.

Fig. 2

Interaction of Nimbin with Mpro

Gedunin-Mpro complex (Fig. 3) had − 9.51 kcal/mol as the minimum binding energy. A sum total of 7 types of bond formation was observed. Glu166 formed convention H-bond, pi-anion bond and carbon hydrogen bond with gedunin. His41 forms a pi-sigma as well as pi-alkyl bonds while Cys145, Met49, and Leu27 formed alkyl bonds. His172 and Phe140 both formed carbon hydrogen bond. His143 contributed to stabilization by forming an additional pi-alkyl bond remaining all the residues were attracted to the ligand by van der waals bond.

Fig. 3

Interaction of Gedunin with Mpro

Epoxyazadiradione was the other compound with binding energy greater than of control. It showed a binding energy of − 8.47 kcal/mol. Four types of interaction can be observed. Along with numerous residues involved in weak van der waals interaction, Cys145, Gly143 and Thr26 forms H-bond. Cys145 and Met165 formed alkyl bonds and Thr25 forms a carbon hydrogen bond (Fig. 4).

Fig. 4

Interaction of Epoxyazadiradione with SARS-CoV-2 Mpro

Ginsenosides-Mpro complex showed the minimum binding energy of -9.63 kcal/mol among all the conformations and ligands. A total of five different type of stabilizing interactions were observed. A conventional H-bond formation was done by Thr26 and Glu166 residues. His41 formed pi-sigma and pi-alkyl with 3 atoms of the ligand. Leu27, Cys145, Met165 and Met49 helped stabilizing complex via alkyl bond formation. Nine more residues can be observed around the ginsenosides interacting via van der waals forces (Fig. 5).

Fig. 5

Interaction of Ginsesnosides with 6LU7 protein

Mpro belongs to a protease class of proteins and this particular protease is seen to play essential role in replication process of many viruses (Chen 2020). The main protease of SARS-CoV-2 was recently positioned in the pdb database as 6LU7 and provided the first potent target for the drug development process. The structural comparison revealed a high percentage (96%) of similarity between the SARS-CoV-2 and SARS virus main protease. The catalytic site is conversed in SARS-CoV-2 main protease and it is formed of two amino acids: His41 and Cys145 forming a catalytic dyad (Mirza and Froeyen 2020). Catalytic dyad serves its purpose in maturation of the virus thus, making it an active target site in the protein (Chang 2010). It can be observed that except nimbin all the other three ligands interact with the catalytic dyad residues in different fashion. Among the three, gedunin and ginsenosides had direct bond formation with both residues while epoxyazadiradione forming H-bond with only Cys145 and having weak Van der wall interaction with His41. Thus, all the four compounds can hold potential and might help in the process of finding the right cure to the ailment (Fig. 6).

Further analysis included four selected ligands viz., nimbin, gedunin, epoxyazadiradione and ginsenosides. A more illustrated and comprehensive study was done using bioavailability radar. Bioavailability radar is a descriptive tool to look into the drug-likeliness of the ligands based on six physiochemical properties. Study found gedunin and epoxyazadiradione to be orally bioavailable as the ligand radar entirely fit the pink shaded area. Nimbin was not an orally bioavailable ligand as it disobeyed the size parameter. Along with nimbin, ginsenosides having poor solubility in water and high value of lipophilicity was also predicted as not orally bioavailable candidate (Fig. 7).

Fig. 7

Bioavailability radar of Gedunin, Epoxyazadiradione, Nimbin and Ginsenosides. The pink shaded zone is an estimated physicochemical space for oral bioavailability. LIPO (Lipophility): − 0.7 < XLOGP3 <  + 5.0. SIZE: 150 g/mol < MV < 500 g/mol. POLAR (Polarity): 20Å2 < TPSA < 130Å2. INSOLU (Insolubility): 0 < LogS (ESOL) < 6. INSATU (Insaturation): 0.25 < fraction Csp3 < 1. FLEX (Flexibility): 0 < Number of rotatable bonds < 9

Two ligands (which completely satisfy bioavailability radar) i.e., gedunin and epoxyazadiradione were subjected to molecular dynamic simulation for 50 ns. − 117,519 kcal/mol and − 117,476 kcal/mol which obtained potential energies for gedunin-6LU7 and epoxyazadiradione-6LU7 respectively. Total energy was another parameter obtained at the end of the 50 ns simulation. Total energy of − 95403.378 kcal/mol was found for gedunin-6LU7 and − 95365.196 kcal/mol for epoxyazadiradione-6LU7.

Structural deviation and compactness

The mean Root-mean-square deviation (RMSD) of gedunin-Mpro complex and epoxtazadiradione-Mpro complex was 2.017Å and 2.023Å respectively. Figure 8a indicates that RMSD for the gedunin-Mpro complex was increasing for a brief period of time, however, the complex can be seen stabilizing toward the end of the simulation. The RMSD for epoxyazadiradione-Mpro complex can be seen increasing throughout the simulation period. Small crust and troughs in RMSD can be seen scattered throughout the trajectory but they are due to the internal vibrations of the system.

Fig. 8

a RMSD plot of 6LU7 with gedunin and epoxyazadiradione b RMSF plot of 6LU7 with gedunin and epoxyazadiradione c Radius of gyration of 6LU7 with gedunin and epoxyazadiradione d SASA plot of 6LU7 with gedunin and epoxyazadiradione (Blue plot: Gedunin; Orange plot: Epoxyazadiradione)

Root mean square fluctuation (RMSF) plot obtained assisted in understanding the stability of residues in the protein. Residues packed to form secondary structures like beta-sheet and alpha-helix confer less fluctuation thus have less RMSF while the loops and terminal residues have higher RMSF. Figure 8b reveals that both complex’s residues follow the above stated phenomenon. However, most of the residues binding with epoxyazadiradione were observed to have higher RMSF value suggesting instability at the binding void.

The average radius of gyration (Rg) for gedunin-Mpro complex and epoxyazadiradione-Mpro complex were 4.099Å and 4.0813Å respectively. From Fig. 8c it can be seen that Rg value for gedunin-Mpro complex underwent some minor fluctuations, up to 0.1Å, which could be possible due to the packaging of the complex while epoxyazadiradione-Mpro complex underwent some fluctuations, up to 0.15Å. Towards the end of the simulation, gedunin-Mpro complex attained Rg value near to its average value. Epoxyazadiradione-Mpro complex also obtained a value near its average Rg but it experiences more frequent fluctuations throughout the simulation.

The mean solvent-accessible surface area (SASA) for gedunin-Mpro complex and epoxyazadiradione-Mpro complex came out to be 171.412 ± 30.815 Å2 and 269.354 ± 49.952 Å2 respectively. Evident from Fig. 8d, the gedunin-Mpro complex had attained a value about the mean value of SASA in the middle of the simulation. Epoxyazadiradione-Mpro complex showed a very high peak in SASA towards the end of the simulation signifying interaction of hydrophobic and internal residues with the solvent which destabilizes Mpro and is suggestive of greater chances of unfolding of the protein structure.

Interaction dynamics and secondary structural analysis

By analyzing the interactions of gedunin-Mpro complex, all residues observed in docking can be seen interacting with the ligand. Few more additional residues were found interacting with gedunin line THR26, THR24, ASN142, etc. THR26 and GLU166 were interacting with the ligand 100% of the simulation time with two types of bonds namely, H-bond and water bridges thereby, holding the ligand in the binding pocket. Other residues like GLN189, SER144, HIS172 formed H-bond but for a very brief period of time. Other types of interactions were also visible accounting for additional stability of the ligand in the binding void (Fig. 9a). Study of interaction of epoxyazadiradione-Mpro revealed that all the residues observed in docking results showed interaction with the ligand. THR26 and HIS41 were forming H-bond for interaction. THR24 formed around 40% of the time waster bridges. Other residues like MET49, GLN189, GLN306, GLU166 and many more added to the stability of the ligand in the void. Comparing the total number of contacts, the protein molecule was able to form number of bonds with gedunin as compared to epoxyazadiradione. During the 50 ns simulation time with gedunin only one instance was observed devoid of any contacts while with epoxyazadiradione lot of instances can be observed having zero contacts (Fig. 9b). Computing and analyzing secondary structural analysis can be used to understand the protein packing and characteristics of folding of the protein with different ligands. Table 3 shows the percentage number of secondary structures of 6LU7 protein with gedunin epoxyazadiradione ligands. With gedunin 40.23% of total residues formed secondary structure while binding of epoxyazadiradione with the protein led to dropping of number to 40.09%.

Fig. 9

Protein–Ligand contact histogram and timeline of particular interaction a gedunin b epoxyazadiradione

Table 3

Protein secondary structure element (SSE) analysis

Compound	% Total SSE	% Alpha helix	% Beta strand
Gedunin	40.23	16.27	23.97
Epoxyazadiradione	40.09	16.41	23.68

Conclusion

In the present study 38 compounds were selected from five plants. These compounds were screened using Lipinski’s rule of five and determined drug-likelihood of the compound. 31 compounds were drug-likeable which were subjected to molecular docking. Docking results were compared with the native inhibitor of Mpro i.e., N3 and four ligands were found more potent in terms of binding energy. All four ligands were studies using bioavailability radar. Our results proposed gedunin and epoxyazadiradione to be potent inhibitors of the 6LU7 along with nimbin but nimbin was not found orally bioavailable. Similarly, ginsenosides showed the best docking result but not the bioavailability radar result. Further molecular dynamics showed that gedunin is more potent than epoxyazadiradione. To find the effectiveness and to propose the exact mechanism in-vitro studies can be encouraged on gedunin, epoxyazadiradione and nimbin to find a potent cure for the COVID-19.