Identification of phytocompounds as newer antiviral drugs against COVID-19 through molecular docking and simulation based study.

COVID-19 pandemic has emerged as a global threat with its highly contagious and mutating nature. Several existing antiviral drugs has been worked on, without proper results and meanwhile the virus is mutating rapidly to create more infectious variant. In order to find some alternatives, phytocompounds can be opted as good one. In this study, three hundred phytocompounds were screened virtually against two viral proteins namely main protease and spike protein. Molecular docking and dynamic simulation study was used to find binding affinity, structural stability and flexibility of the complex. Pharmacokinetic properties were studied through ADMET analysis. To understand energy variation of the complex structure free energy landscape analysis was performed. Among three hundred phytocompounds virtual screening, three phytocompounds were selected for detailed molecular interaction analysis. Oleanderolide, Proceragenin A and Balsaminone A, showed strong binding affinity against both the target proteins and reflected conformational stability throughout the MD run. Oleanderolide, proceragenin A and balsaminone A has docking score -9.4 kcal/mol, -8.6 kcal/mol, and -8.1 kcal/mol respectively against main protease and same -8.3 kcal/mol docking score against spike protein. These three phytocompounds has high gastrointestinal absorption capacity. They were unexplored till now for their antiviral activity. Their promising in silico results suggests that they can be promoted in the long run for development of new antiviral drugs.

Introduction

The coronavirus disease (COVID-19) is a respiratory illness disease caused by SARS-CoV2, was first reported in Wuhan, China in December 2019 [20]. Since then more than 4 million people were died because of the disease and as of 3rd September 2021 the number of active cases around the world is more than 2 million (https://www.worldometers.info/coronavirus/).

This worldwide spreading novel virus is named SARS-CoV-2 belongs to the family coronaviridae. It has a genomic similarity of more than 80% with previously reported SARS-CoV virus. The coronavirus family (CoV) has single-stranded positive sense RNA as genetic material which is further surrounded by crown like spike protein (in Latin corona) [23]. The spike protein is responsible for viral entry into the host cell.

It is believed that bats are the most likely natural host for CoVs. Historically it was believed, that coronavirus was causing respiratory and enteric illness only in farm animals, but the emerging coronaviruses (beta-coronavirus) are proved to be pathogenic to humans as well, causing severe respiratory and gastrointestinal symptoms [3]. Unlike the previous pathogenic CoVs, the global pandemic causing SARS-CoV2 has a range of 0.2–9% with an average of 2% case fatality ratio, which also varied in different countries, but the highly contagious nature of the virus had made it the most consequential global health crisis. The pandemic has exposed the bottlenecks in global healthcare system to handle sudden outbreaks and the global economic structure is also shattered because of the prolonged shutdown [13].

The most common symptoms of COVID-19 are shortness of breath, cough, fever, fatigue, body ache, and in extreme cases it may result in lungs infection, kidney and multi-organ failure. Cohort studies have demonstrated that patients with co-morbidity have increased risk of death [38]. In this extended period of the pandemic, COVID-19 has also developed several variants with increased transmissibility or virulence which has increased concern for the medical world. Although several studies were undertaken to find inhibitors against COVID-19, but none of them demonstrated much effectiveness.

In search for inhibitors against COVID-19, computer aided drug discovery (CADD) has emerged as a powerful tool in the drug discovery process, which reduces the time of screening of potential hit compounds in the task of finding lead drug candidates. In this study, we have selected the main protease (Mpro) protein because it is highly conserved in coronaviruses and plays important role in viral replication. Besides that, no proteases have been reported in humans with the same cleavage sites, so finding an inhibitor against Mpro is safe for human [37]. Spike protein has been chosen for another drug target because of its role in both receptor recognition (through S1 subunit) and membrane fusion (through S2 subunit). Spike protein binds with the acetyl-choline esterase (ACE2) receptor in the host cell and their increased affinity is the reason for the contagion nature of COVID-19 [18].

In this study, we have targeted two major viral proteins (Mpro and spike) and performed virtual screening; molecular docking and all-atoms molecular dynamics simulations with phytocompounds reported from twenty-five Indian medicinal plants and performed a detailed interaction study, which may lead us to select the potent compounds for further the clinical trial.

Materials and methods

To start with the virtual screening process, the crystallographic structures of the target proteins (6LU7 and 6M0J) were obtained from RCSB Protein Data Bank (PDB). COVID-19 main protease (Mpro) with PDB ID- 6LU7 and spike receptor binding domain (RBD) with PDB ID-6M0J were downloaded for further experimental analysis.

Receptor protein preparation

The X-ray crystallographic structures of the target proteins were selected based on their Resolution (6LU7:2.16A°, 6M0J: 2.45A°), R-value (6LU7:0.235A°, 6M0J: 0.227A°), and no mutation in the structure [19]. The protein structures were first checked through Swiss PDB viewer for any missing atoms, then further optimized using AutoDock tool by deleting all water molecules, heteroatom's and inhibitors followed by addition of polar H atom, Kollman charge and computation of Gasteiger charge. Finally, the 3D coordinates of the optimized protein structures bearing partial charges were saved in PDBQT format for the next step of virtual screening using AutoDock Vina [35].

Preparation of ligands

In consideration of the significant importance of Indian medicinal plants in controlling diseases, we have selected twenty-five commonly found Indian medicinal plants for screening their activity against COVID-19. We have searched phytocompounds from the respective plants through the IMMPAT database which is the largest repository of Indian medicinal plants [25]. The three-dimensional structures were saved as SDF files, which were further processed through PyMOL to convert into PDB files and optimized using AutoDock tools. For the compounds with only two-dimensional structures (available in PubChem), the three-dimensional structures were created through VegaZZ software and followed the previous process of optimization [27].

Docking of phytocompounds against target proteins

Molecular docking was performed to assess the binding energy of the protein-ligand interaction and to provide initial coordinates for the molecular dynamic simulation. The saved PDBQT proteins files were separately incorporated into AutoDock tools 1.5.6. The information regarding active site residues was collected from a previous literature study [31]. The active site residues were marked to draw the grid box/search area with the spacing of 0.375 Å. For main protease (Mpro), a grid box with dimension 66X66X66 Å3 was centred at the coordinates X = −13.099, Y = 15.875, Z = 70.369 and for the spike protein RBD (6M0J) a grid box with dimension 50X50X50 Å3was centred at the coordinates X = −32.231, Y = 15.038, Z = 28.403. For both the proteins, the selected ligands were incorporated individually and with the in-house scripts that were incorporated in the Autodock Vina program, a virtual screening process was performed. The 300 phytocompounds from twenty-five plants were screened and ranked according to the binding energies. The docking results were confirmed by re-docking the ligands with corresponding proteins. After docking, the structures were visualized using PyMOL (www.pymol.org). Antiviral drugs Flavipiravir, Ritonavir were used as reference standards for selecting phytocompounds [2]. For best 1% (30) compounds, docking scores and interactive amino acid residues were listed in Table 1 .

Table 1

List of various interaction and interacting residues in the active site of Mpro and RBD with potential natural compounds, were listed from docked complexes.

Sl No	Compound Name	Main protease (6lu7)			RBD (6M0J)
		Docking score (kcal/mol)	H bond	Hydrophobic interactions	Docking Score (kcal/mol)	H bond	Hydrophobic interactions
1.	Flavanomarein	−8	T29,T190,Q192, Q189,R188	L27,L141,G143, E166,M165,N142	−7.5	C336,G339,F436	N343,F338,F342,L368,
							S371,L335,V367
2.	Pindrolactone	−8.3	R131,D289,T199	M270,G275,Y239, L271,L287,L280	−8.1	W436	F374,S371,S373
3.	Cyanidin-3-glucoside	−8.3	T26, Y54, F140, L141,S144,C145, H163	N142,E166,M165,D187,Q189,R188,G143,M49,H41,L27	−7.3	R355,D428,T430, P463,F525	P426,F429,S514,F464
4.	Diosgenin	−8.3	S144,L141	A191,P168,D166,N142,H163,C145,M165,Q189,T190	−8.1	–	V362,D364,L335,V367, C336,A363,L368,F342, S371,N343,S373
5.	β-Amyrin	−7.9	K137	D197,T198,N238,T199,Y237,L272,L271,L287, Y239	−8.2	–	V367,G339,F342,N343,W436,L441,S371,S373,F374,L368
6.	Fumarprotocetraric acid	−8.3	G143,C145,E166,D187	L141,H164,N142,Q189,M165,H41,P52,M49,Y54,C44,R188	−7.6	F347,S349,Y351,N354,N450	A352,Y451,A348,R346, DV341,S399,A344
7.	α - Amyrin	−8	K137, D197	L271,L286,L287,L272,Y237,T198,T199,Y239	−8.5	–	V367,L368,G339,S371,N343,W436,S373,N440,L441,F374,F342
8.	Bryonolic acid	−7.9	K137,D289,M276	R131,T199,Y239,L272,L287,L286,G275	−8.2	T430,F515	R355,P463,F464,P426,F429,S514,D516
9.	Citrullonol	−9.1	R105	V104,D153,I152,F294,N151,Q110,F8,S158	−8.1	F515	D428,P426,T430,S524,F464,G516,R355
10.	Quassinoid analog	−8.4	T111,N151,S158	I152,D153,I106,Q110,V104,F294	−7.3	–	D428, P426, T430,S514,F464,F515,E516,R355,Y396
11.	Astragalin	−8.4	Y54,L141,S144,H163,D187,R188,T190,Q192	P168,E166,E89,H41,C145,M165,F140	−6.7	E516,L517,F515,P463	L518,F464,P426,D428,F429,T430,S514
12.	Cynaroside	−8.1	T24,T25,T26,G143,S144,R188	M49,T45,Q189,E166,D187,M165,H164,N142,C145,H41	−7.4	R355,T430,R466,F515	D516,F464,P463,F429,S514,F429,W353
13.	Epicatechin gallate	−8.2	F464,S514,F429,P426	R355,F515,T430,D428,P463	−7.0	E166,F140	H172,H163,N142,S144,H164,M165,R188,G143,L141,T26,H41,L27,D187,Q189,T25,M29
14.	Astilbin	−8.4	N142,G143,Y54,E166,D187	M49,H41,R188,Q189,C145,M165,L141	−6.9	T430,P463	F464,D428,P426,S514,F525,F429
15.	Baicalin	−8.3	T24,T25,T26,C44,S144,C145	S46,T45,G143,N142,M165,Q189,M49	−7.7	E340,A348,N354	A352,R346,F347,S399,A344,V341
16.	Guajavrin	−8.3	L141,S144,C145,H163,E166	H41,D187,N142,M165,Q189,R188,G143,M49	−7.1	N343,A344,S373,R509	S371,W436,N440,F342,F374,L441
17.	Balsaminone A	−8.1	L287	L272,T199,Y237,Y239,R131,L286,K137,D287	−8.3	-	S373,F374,W436,S371,F338,N343,F342,V367,L368,G339
18.	Canophyllic acid	−8	R131,M276	A285,L271,T199,L287,K137,D286,D289,G275	−8.2	T430,D428	R355,F464,P426,S514,F515,F429,P463
19.	Oleanderolide	−9.4	L271,L272,L287	D289,L286,Y237,Y239,T199,G275	−8.3	K378	G504,V503,D405,R408,V407,T376,G404
20.	Oleanolic acid	−8.5	R131,D289	K137,L287,L286,Y237,Y239,T199,L272	−8.2	S373	V367,S371,L368,F342,N343,F374,W436,L441
21.	Proceragenin A	−8.6	L271	Y237,Y239,L272,T199,L287,L286,K137	−8.3	R355,T430	D428,F429,D516,L517,Y396,F464
22.	Vitexin	−8	Y54,D187	H41,M49,H164,M165,D166,R188,Q189,T190	−7.2	R355,F515	T430,P426,P463,F429,F464,Y396
23.	Betulic acid	−7.7	T199,N238	L271,L272,L286,L287,D197,R131,T198,Y237	−7.4	Y508	A411,Y380,G404,T376,407,S375,V503,D405,R408
24.	Avicularin	−8.8	Y54,L141,N142,E166,D187	G143,Q189,H41,F140,C145,R188,M49,H163,L27,H172,M165	−7.2	D428,T430,P463,F515,E516	R355,Y396,F464,F429,S514,P426
25.	Ursolic acid	−7.7	T199,M276	G275,A285,N238,K137,Y239,L287	−8.3	S373	N343,L441,N440,L368,G339,F342,F338,W436,S371,F374,V367
26.	Subulin	−7.9	T24,T26,S46,G143,E166,T190	T25,M49,N142,Q189,P168,R188,Q192,H163,M165	−7.8	L335,C336,G339,N343,S371,S373	F374,L368,V367,F338,P337,W436,L441,F342
27.	Silybin	−7.9	G143,H164,T190	M165,E166,S144,C145,L141,N142,P168,Q189,R188,H41	−7.6	D428,T430,R466,F515	E465,F464,R355,E516,S514,P426
28.	Hexanorcurcubitacin	−7.7	T111	Q110,N151,F294,S158,I106,V104	−8.8	R355,T430,S514	F429,Y396,E516,D428,P426,F464
29.	Aromadendrin	−7.6	Y54,L141,G143,S144,E166,D187	N142,C145,H164,H41,R188,M49,Q189,M165	−6.5	R355,T430,S517	F464,P426,E516,F515
30.	Clausenin	−6.8	Q107, S158	F8,I106,V104,Q110,N151,D153,F294	−6.9	–	F342,N343,S371,F374,L368,W436,L441
31.	Ritonavir	−6.4	R131,Y237	N238,D289,T199,T198,D197,L272,Y239,L286,L287	−6.1	S373,N343,G339	E340,P337,L335,V367,D364,C336,L368,F342,F338
32.	Favipiravir	−5.7	Q110,N151	F294,T111,T292	−5.3	A348,S399	V341,R346,E340,A344,N354,F347

ADMET calculation

The previously selected thirty phytocompounds were also evaluated for their ADMET (Absorption, Distribution, Metabolism, Excretion, Toxicity) properties through SwissADME and pkCSM server [10,28]. Lipinski Rule of five (R05) violation, GI absorption, BBB permeability, hepatotoxicity, oral rat acute toxicity (LD50), synthetic and bioavailability scores were into consideration for selection of the compounds for MD simulation. The pharmacokinetic parameters help to understand the reaction after drug entry into the body. ADMET analysis helps to narrow down the search process for selecting candidate drugs.

Molecular interaction analysis

The selected phytocompounds were studied their protein-ligand interacting atoms using LigPlot+ v1.4.5 [36] and PLIP [1]. The structures were visualized and illustrated through PyMOL and Adobe Illustrator. Non-bonded interactions like hydrophobic interactions, hydrogen bonds, salt bridges were marked for each docked complex. These results help us to identify important interacting amino acid residues and the mode of interaction.

Molecular Dynamic (MD) simulation

From the ADMET analysis, we have selected three compounds were effective against both Mpro (6LU7) and spike (6M0J) proteins. The three selected phytocompounds along with the plants’ name, from where they have been collected are – balsaminone A (Impatients balsamina), oleanderolide (Nerium oleander),proceragenin A(Calotropis procera).

The Molecular dynamic simulation was performed with those six complexes to study the dynamic behavior, stability, and conformational flexibility using Gromacs 2018.3 version [4,29]. Charmm36 force field was applied to generate topology for proteins and ligands [22]. A cubic water box with TIP3P water model was used for structure solvation. The complex structure was electro-neutralized by adding 0.15 M NaCl. The structure was energy minimized using the steepest descent algorithm in 5000 steps. After energy minimization, the structure was equilibrated through NVT and NPT ensembles for 5ns and 10ns respectively. Finally, MD simulation was performed for 100ns at 300K with 2fs time step using Leap-Frog integrator [8].

MD trajectory analysis

The quality of the trajectories was analyzed using GROMACS utility toolkit through determining conformational flexibility, stability of the trajectories. Quality assurance parameters like root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and intermolecular H Bond, were calculated using gmx rms, gmx rmsf, gmx gyrate, gmx sasa, gmx hbond built-in commands in Gromacs respectively. XmGrace, PyMOL, BIOVIA Discovery Studio Visualizer were used for visualization of 2D graphs and structures respectively and the plots were created in Adobe Illustrator.

Principal Component (PC) analysis

Principal component analysis (PCA) is a multivariate statistical analysis tool that is used to reduce the dimensionality of a dataset while preserving the variability as much as possible. This method is used to understand the correlated motion of the residues to a set of linearly uncorrelated variables which are called (PC) principal components [12]. To analyze principal component, gmx covar and gmx anaeig built-in modules of GROMACS were used. In GROMACS, gmx covar tool is used to construct the covariance matrix. The matrix generates diagonal eigenvectors (PCs) which represent the correlated motion of the protein and eigenvalues indicate atomic contribution of the motion of protein-ligand complex system. PC Analysis was done by calculating eigenvectors along with eigenvalues and associated projections of the first two principal components where gmx anaeig tool was used [34]. Here covariance matrix was constructed from the last 50ns of MD trajectory of the protein-ligand complex system and the first 20 eigenvectors were used for representation. The porcupine plot was generated using modevectors.py script using PC1 and PC2 as parameters.

Free energy landscape analysis

Free energy landscape (FEL) is a popular tool in bioinformatics to analyze the folding or aggregation of the complex after MD simulation, which in terms gives an idea about the stability of the protein [15]. Principal component values were used as order parameters for the calculation. The gmx sham module was used to obtain the FEL of each complex using the whole MD trajectories. In this study, the free energy landscape of protein-ligand complexes was analyzed using Gibbs free energy calculation equationGi = - KWhere KB stands for Boltzmann's constant, T for temperature (300K), Ni for population bin I, and Nmax depicts the most inhabited population [11]. A color model was used for representation of the energy level in plot, where red and blue showed the highest and lowest energy level respectively.

Cluster analysis

Geometric clustering (gmx cluster) was performed to explore heterogeneity in the protein structure after performing MD simulation. Gromos clustering algorithm by Daura et al., 1999 was employed. During analysis, Cα RMSD cut off 0.2 nm was set for identifying structurally similar clusters and the top 2 ranked clusters with maximum screenshots were analyzed in detail and compared with docked conformations [30].

Results and discussion

The COVID-19 pandemic has become a global threat that impacted human health care facilities and the socio-economic structure of the world. Studies have shown that, during the first wave of COVID-19, there was increased use of azithromycin for treatment, but later it was declared by WHO and CDC that, antibiotics have no direct impact on COVID-19. Although some vaccines have arrived in the market, but the effectiveness of the vaccines and fighting the pandemic through only vaccination is not recommended. Besides in largely populated countries complete vaccination will take a lot of time and the virus is evolving rapidly by developing variants, which are proved to be more deadly in nature [16]. Besides that, with COVID-19 prevalence the intake of antibiotics has also increased tremendously [24].

So it is high time to find some drugs for treatment of the infected persons with confidence, so that the spread of the virus can be blocked. In this present study, we have selected Mpro and spike protein (RBD) as drug targets, because both of them are essential for viral entry and replication. Case studies reported effectiveness of traditional Chinese medicines in the containment and prevention of COVID-19 infection in China [9]. Indian medicinal plants are also well known for their medicinal properties. Dual targeting of proteins can be a good way to select phytocompounds for next step of clinical trial. In virtual screening, molecular docking followed by MD simulation has become a classic way to screen the inhibitors. In MD simulation, various utility tool kit has been developed, which help to understand detailed intermolecular interaction, conformational flexibility, energy variation of the system, that gives a very clear prospect of the suppose mechanism of action of the compound against target protein.

Screening of phytocompounds against SARS CoV-2 viral proteins

Molecular docking between protein and ligand reflects binding affinity between them and the lowest docking score is considered to have the highest binding affinity. In our study, 300 phytocompounds were categorized into three categories for each protein – category I (−3 to -5 kcal/mol), category II (−5.1 to -7 kcal/mol), and category III (−7.1 kcal/mol to above) (S1).

Among 300 phytocompounds, 161 phytocompounds for Mpro protein and 150 phytocompounds for spike protein (RBD) belong to category II. The docking scores of ritonavir against Mpro and spike protein (RBD) were −6.4 kcal/mol and −6.1 kcal/mol respectively. Another antiviral drug favipiravir has a lower binding affinity than ritonavir. The docking scores of favipiravir against Mpro and spike protein were −5.7 kcal/mol and −5.3 kcal/mol respectively. Among 300 compounds, we have selected the best 30 compounds for further analysis. The docking scores and ligand efficiency values were listed in Table 1.

All of these compounds have shown higher binding affinity than the standard reference antiviral drugs. Among those thirty compounds, Oleanderolide (−9.4 kcal/mol) and α -Amyrin (−8.5 kcal/mol) has shown the highest binding affinity against Mpro and RBD protein respectively. In comparison with other phytocompounds, Clausenin (−6.8 kcal/mol) has low binding affinity against Mpro and Aromadendrin (−6.5 kcal/mol) has the lowest binding affinity against RBD. But among the above listed compounds, 8 phytocompounds were found to have docking score ≤ −8.0 kcal/mol.

Evaluations of ADMET properties

The results of ADMET analysis for the best 30 phytocompounds were listed in Table 2 . It was stated that any drug which violates two or more rules of Lipinski Rule of Five (R05) will not be orally active [17]. Six compounds have been found to fulfil Lipinski's Rule of five. Quassinoid analog, Pindrolactone, Bryonolic acid, Citrullonol, Canophyllic acid, Oleanolic acid, Ursolic acid, Betulic acid and Clausenin are hepatotoxic and rests three Aromadendrin, Guajaverin, Subulin has very low intestinal absorption capacity. Based on the high binding affinity against the target proteins and non-hepatotoxicity, high gastrointestinal absorption (>90%), the three phytocompounds were chosen for further analysis.

Table 2

ADMET analysis of the phyto-compounds employed in this study.

Sl No	Compound name	Chemical formula	nVio	Hepatotoxicity	LD50 (mol/kg)	GI absorption (%)	BBB	SA	Bioavailability Score
1	Flavanomarein	C21H22O11	2 vio	no	2.73	40.3	no	5.06	0.17
2	Pindrolactone	C30H42O3	1 vio	yes	1.84	95.8	no	6.33	0.55
3	Cyanidin-3glucoside	C21H21ClO11	2 vio	no	2.54	29.9	no	5.3	0.17
4	Diosgenin	C27H42O3	1 vio	no	1.92	96.5	yes	6.94	0.55
5	Aromadendrin	C15H12O6	No vio	no	2.16	59.07	no	3.42	0.55
6	Beta-amyrin	C30H50O	1 vio	no	2.47	93.7	no	6.04	0.55
7	Fumarprotocetraric acid	C22H16O12	1 vio	no	2.43	28.7	no	3.87	0.11
8	Alpha.-Amyrin	C30H50O	1 vio	no	2.46	94.07	no	6.17	0.55
9	Bryonolic acid	C30H48O3	1 vio	yes	2.58	98.1	no	5.91	0.85
10	Citrullonol	C30H42O2	1 vio	yes	1.94	96.0	no	6.62	0.55
11	Quassinoid analog	C27H38O7	No vio	yes	4.04	80.0	no	6.77	0.55
12	Clausenin	C14H12O5	No vio	yes	2.01	95.3	yes	2.89	0.55
13	Astragalin	C21H20O11	2 vio	no	2.54	48.05	no	5.29	0.17
14	Cynaroside	C21H20O11	2 vio	no	2.54	37.5	no	5.42	0.55
15	(−)-Epicatechingallate	C22H18O10	2 vio	no	2.55	62.09	no	4.16	0.55
16	Astilbin	C21H22O11	2 vio	no	2.58	49.0	no	5.27	0.17
17	Baicalin	C21H18O11	2 vio	no	2.63	26.2	no	5.09	0.11
18	Guajavarin	C20H18O11	No vio	no	2.58	51.8	no	5.05	0.17
19	Balsaminone A	C21H12O5	No vio	no	2.43	95.60	no	3.44	0.55
20	Canophyllic acid	C30H50O3	1 vio	yes	2.26	100	no	5.27	0.85
21	Oleanderolide	C30H48O4	1 vio	no	2.25	92.9	no	6.64	0.55
22	Oleanolic acid	C30H48O3	1 vio	yes	2.34	99.9	no	6.08	0.85
23	ProcerageninA	C30H46O4	1 vio	no	2.15	91.9	no	6.92	0.55
24	Vitexin	C21H20O10	2 vio	no	2.59	46.69	no	5.12	0.55
25	Avicularin	C20H18O11	1 vio	no	2.54	57.2	no	5.63	0.85
26	Betulic acid	C30H48O3	1 vio	yes	2.25	99.76	no	5.04	0.17
27	Ursolic acid	C30H48O3	No vio	yes	2.34	100	no	6.21	0.85
28	Hexanorcucurbitacin I	C24H32O5	3 vio	no	2.33	90.7	no	5.74	0.55
29	Subulin	C28H32O16	No vio	no	2.52	32.6	no	6.49	0.17
30	Silybin	C25H22O10	2 vio	no	2.55	61.8	no	4.92	0.55

As of now, several studies have employed virtual screening of phytocompounds against SARS-CoV-2 proteins [6,26] to discover potent compounds against Covid-19. Although Oleanderolide and Proceragenin A has slight violation for one Rule of Lipinski (mLogP> 4.15) but rest of the parameters justify the reasons for this selection.

Studies have shown that the most commonly found bonds between protein-ligand complex are – hydrophobic interaction and Hydrogen bonds [14]. In the present study, six docked complexes were analyzed in detail for molecular interaction (Fig. 1 ). Oleanderolide has −9.4 kcal/mol docking score against (6LU7) Mpro protein and it forms two hydrogen bonds with Thr199 (3.0 Å) and Leu287 (3.07 Å) residues of the protein. There are four hydrophobic interactions with Tyr237, Tyr239, Leu286 and Leu287 residues.

Fig. 1

Non-bonded interaction of phytocompounds with SARS-CoV-2 Mand RBD obtained after docking with Autodock Vina [A: Mpro-Balsaminone A; B: Mpro-Oleanderolide; C: Mpro-Proceragenin A; D: RBD-Balsaminone A; E: RBD-Oleanderolide; F: RBD-Proceragenin A]. The protein is shown in cartoon representation and the ligands are shown in stick format while interacting amino acids (shown in stick format) are labelled. The dotted lines indicate the hydrogen bonds.

Proceragenin A with a docking score of −8.6 kcal/mol with Mpro and has developed three H bonds with Arg131 (3.60 Å), Thr199 (2.89 Å) and Leu271 (3.02 Å) and created six hydrophobic interactions with Tyr237, Tyr239, Leu272, Leu286 and Leu287 residues.

Another studied complex (6lu7-Balsaminone A) has a docking score −8.1 kcal/mol which has formed two H bonds with Asp197 (3.54 Å), and Leu287 (2.87 Å), and three hydrophobic interactions with Tyr237, Leu272 and Leu286 residues of the protein. Balsaminone A also has created two salt bridges with Arg131 (4.51 Å), and Lys137 (4.34 Å). Salt bridge is proved to be the most strong interaction in comparison with H bonds and hydrophobic interactions [14].

All the three above-mentioned phytocompounds have shown good binding affinity with the same docking score −8.3 kcal/mol against spike-RBD protein (6M0J). So molecular interactions were also studied for those complexes.Oleanderolide-6M0J has formed one H bond with Ser375 (4.03 Å) one hydrophobic interaction with Val503 and one salt bridge interaction with Lys 378. Proceragenin A has created three hydrogen bonds with Arg355 (2.88 Å), Asp428 (4.01 Å) and Thr430 (3.07 Å) residues and three hydrophobic interactions with Thr430, Phe464 residues of the protein. Another phytocompound BalsaminoneA formed two H bonds with Ser371 (2.50 Å) and Ser373 (3.04 Å) residues and six hydrophobic interactions with Phe338, Phe342, Val367, Leu368, Phe374 and Trp436 residues of the protein.

Hydrogen bonds are suggested to increase the binding affinity of the protein-ligand complex by removing protein bound water into a bulk solvent [7]. Hydrogen bond interactions of selected phytocompounds with target proteins were shown in S2. The above mentioned six docked complexes were further subjected to molecular dynamic simulation to understand the stability of the docked complex.

Analysis of MD trajectories

After performing molecular docking, the binding affinity between protein and ligand has become clear in a rigid environment, but the binding properties in the solvated state can fully be analyzed only through MD simulation. A total 100 ns MD run was performed to understand structural stability of each complex.

Root mean square deviation (RMSD) of the protein backbone denotes the stability of the structure during MD run. Three protein-ligand complexes (with 6LU7) has found little fluctuations in RMSD till first 20ns, then the system became almost stable with an average RMSD-0.2 nm (Fig. 2 A). For spike protein-ligand complexes (RBD), the structures were seemed to fluctuate till the first 40ns, and after 80ns, the systems appeared to be stable. Among three protein-ligand complexes, Balsaminone A-6M0J has shown the least RMSD value for protein backbone, which indicates more structural stability (Fig. 2D). Similar to protein RMSD, ligand RMSD values were calculated (S3 C,D). The probability distribution of RMSD for each complex was shown in S3 A,B.

Fig. 2

Intrinsic dynamics stabilities of the top ranked phytocompounds (Balsaminone A, Oleanderolide and Proceragenin A complexed with SARS-CoV-2 MPro and RBD during all-atoms MD simulation of 100 ns. Colours mentioned on the bottom panel are followed in the upper images as well. (A) Dynamics stabilities of the SARS-CoV-2 MPro-phytocompound complex systems measured by the plotting the RMSD of the backbone atoms with respect to the initial structure used for production MD over the time scale of 100 ns. (B) Radius of gyration of SARS-CoV-2 MPro over the time scale of 100 ns. (C) The average Cα-RMSF profile of the each amino acid during the last 50 ns of MD in SARS-CoV-2 MPro-phytocompound complexes. (D) Backbone RMSD displaying stability of SARS-CoV-2 RBD-phytocompound complex systems with respect to the initial structure used for production MD over the time scale of 100 ns. (E) Radius of gyration of SARS-CoV-2 RBD complexes over the time scale of 100 ns. (F) The average Cα-RMSF profile of the each amino acid during the last 50 ns of MD in SARS-CoV-2 RBD-phytocompound complexes. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Radius of gyration (Rg) indicates the compactness of the system. Stability in the graph indicates conformational stability of the complex structure during MD simulation. Among three protein-ligand complexes with 6LU7 protein, ProcerageninA showed the most structural stability with an average compact of gyrate-2.25 nm and Oleanderolide was found to be conformationally stable after 40ns (Fig. 2B). But all the protein-ligand complexes with 6M0J protein showed structural stability during the whole MD simulation with an average compact radius of −1.85 nm (Fig. 2E). The steady value of Rg indicates stable folded complex structure, but variations in Rg over time indicate misfolding of the complex structure. A higher value of Rg indicates more flexibility of the conformation. Additionally, solvent-accessible surface area (SASA) indicates the contracted nature of the protein (Figure. S3 E,F). In trajectory analysis, the structural fluctuation and degree of flexibility of the protein were analyzed through (RMSF) of Cα atoms of protein-ligand complexes. A close examination of Fig. 2C,F indicates all the six complexes for both the proteins showed a similar pattern in RMSF. A close inspection of the residues in putty representation, revealed that, loops showing more fluctuations and more structural mobility (S4). Ligand binding residues with higher peaks indicate their mobility and participation in ligand binding.

Intermolecular H-bonds and PCA analysis

Intermolecular stability providing H bond dynamics for each complex were depicted in Fig. 3 (A,B). During the 100ns MD run, all of the complexes showed different intermolecular H bonding patterns, but among them, 6LU7-Oleanderolide, 6LU7-Balsaminone A, 6M0J-Oleanderolide, and 6M0J-Proceragenin A yielded a higher number of hydrogen bonds in comparison with other complexes. During the first 40 ns, an increase in the number of H bonds was observed, but later on a constant of one or more H bonds was observed for each system.

Fig. 3

Dynamics of intermolecular hydrogen bonds and principal component analysis of the SARS-CoV-2 Mand RBD complexes with Balsaminone A, Oleanderolide and Proceragenin A during all-atoms MD simulation of 100 ns. In the left panel black represents 6LU7+BalsaminoneA, red represents 6LU7+Oleanderolide, green represents 6LU7+ProcerageninA and in the right panel blue represents 6M0J + BalsaminoneA, violet represents 6M0J + Oleanderolide, and cyan blue represents 6M0J + ProcerageninA] (A) Intermolecular H-bond dynamics in SARS-CoV-2 MPro-phytocompound complexes during 100 ns MD simulation. (B) Intermolecular H-bond dynamics in SARS-CoV-2 RBD-phytocompound complexes during 100 ns MD simulation. (C) Plot displaying the eigenvalues vs. eigenvectors obtained from the main-chain atoms covariance matrix constructed from the last 50 ns MD trajectory of SARS-CoV-2 MPro-phytocompound complex systems (only the first 20 eigenvectors are considered for representation). (D) Plot displaying the eigenvalues vs. eigenvectors obtained from the main-chain atoms covariance matrix constructed from the last 50 ns MD trajectory of SARS-CoV-2 RBD-phytocompound complex systems. (E) Projection of the motion of the SARS-CoV-2 MPro-phytocompound complexes in phase space along the PC1 and PC2. (F) Projection of the motion of the SARS-CoV-2 RBD-phytocompound complexes in phase space along the PC1 and PC2. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Principle component analysis (PCA) is used to study the collective motion of each complex, based on eigenvectors and eigenvalues. To get a better knowledge about the conformational flexibility of each atom, MD trajectories were thrown to a phase space to yield a spectrum of EVs, which represents a single component of the motion, symbolic of the direction of the motion [5]. A close inspection of each docked complex showed a sharp decrease in the eigenvalues after the first two among the twenty eigenvalues. This indicates considerable flexibility of the docked protein-ligand complex structure during the initial stages of simulation, which is further reduced with simulation interval. Oleanderolide in complex with 6LU7 and Proceragenin A in complex with RBD showed a higher trace value as compared to the rest (Fig. 3C and D). High trace value indicates greater flexibility, which perfectly correlates with the higher occupancy of the phase space occupied in the scatter plot (Fig. 3E and F). This result perfectly supports RMSF analysis data.

To get a better understanding, porcupine plots were generated, by capturing eigenvectors on extreme projections on both PC1 and PC2 where the direction and length of the arrow reflects motion direction and strength respectively. The large arrows indicate increased motion of the complex structures. Here in Fig. 4 , porcupine plots were generated using both PC1 and PC2 order parameters for the protein (Mpro) in combination with phytocompounds. The Mpro-Oleanderolide complex showed more outward arrows in the middle domain for both PC1 and PC2 and the C-terminal (red) domain showed both inward and outward arrows. In comparison to Mpro protein, RBD showed less motion in porcupine plot (Fig. 5 ). However, similar to the scatter plot results, both 6LU7-Oleanderolide, and RBD-Proceragenin A showed a high degree of movement (Fig. 4, Fig. 5). A great movement was observed in the beta-strand connecting loop region in both the complexes.

Fig. 4

PCA displaying global motions occupied by the top two principal components in SARS-CoV-2 MPro-phytocompound complexes during the last 50 ns MD. (A) Porcupine plot displaying the motion of SARS-CoV-2 MPro-Balsaminone A complex obtained from PC1. (B) Porcupine plot displaying the motion of SARS-CoV-2 MPro-Oleanderolide complex obtained from PC1. (C) Porcupine plot displaying the motion of SARS-CoV-2 MPro-Proceragenin A complex obtained from PC1. (D) Porcupine plot displaying the motion of SARS-CoV-2 MPro-Balsaminone A complex obtained from PC2. (E) Porcupine plot displaying the motion of SARS-CoV-2 MPro-Oleanderolide complex obtained from PC2. (F) Porcupine plot displaying the motion of SARS-CoV-2 MPro-Proceragenin A complex obtained from PC2. obtained from PC2.

Fig. 5

PCA displaying global motions occupied by the top two principal components in SARS-CoV-2 RBD-phytocompound complexes during the last 50 ns MD. (A) Porcupine plot displaying the motion of SARS-CoV-2 RBD-Balsaminone A complex obtained from PC1. (B) Porcupine plot displaying the motion of SARS-CoV-2 RBD-Oleanderolide complex obtained from PC1. (C) Porcupine plot displaying the motion of SARS-CoV-2 RBD-Proceragenin A complex obtained from PC1. (D) Porcupine plot displaying the motion of SARS-CoV-2 RBD-Balsaminone A complex obtained from PC2. (E) Porcupine plot displaying the motion of SARS-CoV-2 RBD-Oleanderolide complex obtained from PC2. (F) Porcupine plot displaying the motion of SARS-CoV-2 RBD-Proceragenin A complex.

Free energy landscape (FEL) and clustering

The free energy landscape helps to understand conformational variability along with energy minimization. In Fig. 6 , the dark blue color represents energy minima which mean the lowest energy.The degree of dispersion in the plot indicates conformational flexibility. A larger area was observed for 6LU7-Oleanderolide which also supports the previous RMSF, PCA result in terms of flexibility. There are so no significant conformational changes been detected in the complex structures.

Fig. 6

Free energy landscape (FEL) profiles of the SARS-CoV-2 Mand RBD-phytocompound complexes as a function of PC1 and PC2. The colored scale plot shows the free energy profile (kcal/mol) and the dark blue colored profiles are deep minima. [ (A) MPro-Balsaminone A; (B) MPro-Oleanderolide; (C) MPro-Proceragenin A; (D) RBD-Balsaminone A complex; (E) RBD-Oleanderolide complex; (F) RBD-Proceragenin A. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

In Fig. 7, Fig. 8 , we have extracted the least energy conformational state from FEL analysis and compared them with the top ranked cluster which represents no significant structural variation. Thus, it indicates the structural uniformity is maintained throughout the MD run. To identify similar structures throughout the MD simulation trajectory, clustering analysis was performed using Gromacs clustering algorithm with Cα-RMSD cut off 0.2 nm. For 6LU7-Balsaminone A, 11 clusters were found with an average RMSD - 0.22 and for 6LU7-Oleanderolide, 44 clusters were found with an average RMSD -0.33. For 6LU7 - Proceragenin A complex structure, 21 clusters with an average RMSD - 0.25was found. With the same RMSD cut off for RBD protein, Balsaminone - RBD had 2 clusters with an average RMSD -0.17, Oleanderolide-6M0J had 4 clusters with an average RMSD - 0.18 and 3 clusters for ProcerageninA with an average RMSD - 0.17. In the free energy profile (Fig. 7, Fig. 8) violet color represents energy minima, so large energy minima depicts more stability, but separate peaks in violet indicate transition in the conformation of protein towards more thermodynamically stable state. In Fig. 7A-C, 6LU7-Oleanderolide and in Fig. 8A-C, RBD-Balsaminone A depicts more stability which perfectly correlates the previous FEL profile. RBD-Proceragenin A denotes transition towards thermodynamic stability.

Fig. 7

Three-dimensional FEL profile, snapshot from low energy conformational state and clustering analysis of SARS-CoV-2 M-phytocompound complexes. (A–C) FEL profiles of SARS-CoV-2 MPro-Balsaminone A, Oleanderolide complex and Proceragenin A complexes. (D–F) Snapshots extracted from the minima (low energy conformational states) of SARS-CoV-2 MPro-Balsaminone A, Oleanderolide complex and Proceragenin A complexes. (G–I) Structurally superimposed view of the clustered snapshots of SARS-CoV-2 MPro-Balsaminone A, Oleanderolide complex and Proceragenin A complexes obtained from RMSD based clustering employed in GROMACS with a cut-off of 0.2 nm.

Fig. 8

Three-dimensional FEL profile, snapshot from low energy conformational state and clustering analysis of SARS-CoV-2-RBD-phytocompound complexes. (A–C) FEL profiles of SARS-CoV-2 RBD-Balsaminone A, Oleanderolide complex and Proceragenin A complexes. (D–F) Snapshots extracted from the minima (low energy conformational states) of SARS-CoV-2 RBD-Balsaminone A, Oleanderolide complex and Proceragenin A complexes. (G–I) Structurally superimposed view of the clustered snapshots of SARS-CoV-2 RBD-Balsaminone A, Oleanderolide complex and Proceragenin A complexes obtained from RMSD based clustering employed in GROMACS with a cut-off of 0.2 nm.

Structurally superimposed view of clustered snapshots were represented in Fig. 7G-I and Fig. 8G-I. This free energy profile results clearly correlates with the previous RMSD and radius of gyrate (Rg) data for understanding stability of the system.

In computer aided drug design (CADD), understanding electrostatic surface potential of protein and binding interaction of ligands is an important tool which defines binding affinity of intermolecular complex. In Fig. 9 , it is clearly depicted that, ligands which bind on charged patches of the protein, indicates more stability in behavior throughout MD run and shows high binding affinity.

Fig. 9

Electrostatic surface potential maps of the representative conformations extracted from low-energy structural ensembles of each FEL from SARS-CoV-2 MPro(6LU7) and RBD(6M0J)-phytocompound complexes. [(A–C): 6LU7-Balsaminone A, 6LU7-Oleanderolide, 6LU7-Proceragenin A; (D–F): 6M0J-Balsaminone A, 6M0J-Oleanderolide complex, 6M0J-Proceragenin A. The image was generated using APBS-Plugin in PyMOL, where Blue, Red and White color indicates positive, negative and neutral charge surfaces in SARS-CoV-2 MPro and RBD respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Phytocompounds selected in this study were Oleanderolide from Nerium oleander, proceragenin A from Calotropis procera and balsaminone A was derived from Impatiens balsamina plants, which all are common Indian medicinal plants. Similar to our study, Sen et al. [31], also targeted the same two proteins and screened 10 Indian spices plants for their antiviral activity. Naltrexone, a FDA approved drug was targeted against spike RBD and has shown to bind in the cavity of RBD-ACE2 receptor [8]. Another study targets both Mpro and RDRP protein and screened phytocompounds against them. They selected remdesivir as reference drug. They found Mulberrosides and Punigluconin showed high binding affinity against both the proteins and Emblicanin A, Nimbolide against RdRp while Andrographolides, Anolignans, Chebulic acid, Mimusopic acid showed better binding affinity against Mpro as compared with the reference drug. (A [32]. Similarly, in our study we have used ritonavir and favipiravir as reference drugs. The three selected compounds has shown high binding affinity than those reference drugs. Although there were previous studies [21]; S [33]. reported that different phytocompounds can be found effective against COVID-19 proteins, but our study for the first time reporting these three phytocompounds with strong binding affinity against both proteins. The phytocompounds reported in this study also showed positive results in ADMET property analysis, which suggests that there is a high probability that they will perform well in clinical trial. Diversity of phytocompounds tells that still phytocompounds are very much unexplored and they can be opted as better alternative for treatment.

Conclusion

Till date accurate treatment strategy for ongoing COVID-19 pandemic is limited and controlling COVID-19 pandemic through vaccination is challenging due to the coverage and rapid emergence of variants. To cope up with viral variant development speed, in the present scenario discovery of new antiviral drugs has become inevitable. Computational technique with detailed analysis of the complex structure can pave a path for drug discovery. Natural products are non-toxic to human and have immense antimicrobial potential. In the present study, three compounds named – Oleanderolide, Proceragenin A, Balsaminone A appeared as good antiviral drug candidate. Throughout the MD simulation trajectory analysis and principal component analysis, these three compounds displayed a consistency in stability in the complex structure. But among the three phytocompounds, Oleanderolide for Mpro and Proceragenin A for spike-RBD have demonstrated encouraging results. From the above good results found in in silico studies, it can be concluded that, these three novel phytocompounds (Oleanderolide, Proceragenin A, Balsaminone A) has promising future as possible drug candidates against COVID-19 and should be further evaluated in vitro and in vivo clinical trial for their activity. So to find new antiviral drugs, a good combination of computational and medicinal techniques can be effective and which can offer the mankind new antiviral drugs in near future.

Ethics approval

The study was ethically approved by the institutional human ethical committee of ICMR – Regional Medical Research Centre, Bhubaneswar.

Declaration of competing interest

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