Combination of system biology to probe the anti-viral activity of andrographolide and its derivative against COVID-19.

The present study aimed to investigate the binding affinity of andrographolide and its derivative i.e., 14-deoxy-11,12-didehydroandrographolide with targets related to COVID-19 and their probable role in regulating multiple pathways in COVID-19 infection. SMILES of both compounds were retrieved from the PubChem database and predicted for probably regulated proteins. The predicted proteins were queried in STRING to evaluate the protein-protein interaction, and modulated pathways were identified concerning the KEGG database. Drug-likeness and ADMET profile of each compound was evaluated using MolSoft and admetSAR 2.0, respectively. Molecular docking was carried using Autodock 4.0. Andrographolide and its derivative were predicted to have a high binding affinity with papain-like protease, coronavirus main proteinase, and spike protein. Molecular dynamics simulation studies were performed for each complex which suggested the strong binding affinities of both compounds with targets. Network pharmacology analysis revealed that both compounds modulated the immune system by regulating chemokine signaling, Rap1 signaling, cytokine-cytokine receptor interaction, MAPK signaling, NF-kappa B signaling, RAS signaling, p53 signaling, HIF-1 signaling, and natural killer cell-mediated cytotoxicity. The study suggests strong interaction of andrographolide and 14-deoxy-11,12-didehydroandrographolide against COVID-19 associated target proteins and exhibited different immunoregulatory pathways. This journal is © The Royal Society of Chemistry.

Introduction

In December 2019, a severe acute respiratory syndrome caused by novel severe acute respiratory syndrome novel coronavirus 2 (SARS-CoV-2)[1] emerged as a global pandemic from Wuhan city, Hubei province, China. WHO designated this nSARS-CoV-2 infection as Coronavirus disease (COVID-19). The n-SARS-CoV-2 is a highly contagious virus that can be transmitted from person to person,[2] leading to community transmission. COVID-19 became a major global threat by influencing around 212 countries, with almost half a million deaths worldwide.[3] Presently, it has majorly affected subjects with comorbidities and low immunity that are suffering from infectious and non-infectious diseases.[4] Patients with COVID-19, especially those with severe pneumonia, showed substantially lower lymphocyte counts, and severely ill patients exhibited a reduction in CD4+ T cells, CD8+ T cells, and natural killer cells.[5,6] The higher plasma concentrations of several inflammatory cytokines, such as IL-6 and tumor necrosis factor (TNF), were observed in COVID-19 patients.[7] The pathological findings in patients with COVID-19 showed that immune-mediated lung injury was involved in acute respiratory distress syndrome (ARDS).[8] This evidence suggested an immune imbalance in COVID-19, and it was contemplated that the immune modulation can provide some prophylaxis and promising benefit against COVID-19.[9] Also, there is a need to utilize the concept to identify the new therapeutic agent with immunomodulatory action and anti-viral property against the COVID-19 as reported.[10,11]

The effectiveness of treatment based on traditional medicinal plants has been reported during 2003 SARS.[12-15] Therefore, the scientific community has already started studies on medicinal plants, based on their history and traditional uses, as plausible leads in the treatment of COVID-19.[3,16-20] For thousands of years, medicinal plants have played a vital role in managing multiple infectious and non-infectious diseases.[21-23] Among them, Andrographis paniculata (Family: Acanthaceae), also called known as ‘King of bitters’ and ‘Indian Echinacea’ reserves its importance in the management of various infectious and non-infectious diseases.[9,24-26] Further, it has been studied well for its potency as a modulator of the immune system.[25,27] In Andrographis paniculata, andrographolide[28] is a major bioactive that possesses beneficial effects in multiple pathogenic conditions, including the immunity booster role.[29] Further, two important databases, i.e., ChEBI and PCIDB, also record andrographolide (Fig. 1) as chief bioactive from Andrographis paniculata. Andrographolide and its derivative(s) also exhibited decisive immunomodulatory action[25,27] and have broad-spectrum anti-viral properties.[30] Further, it was found to be effective against multiple viral infections like dengue,[31] swine flu,[32] hepatitis C,[33] chikungunya,[24] influenza,[34] Epstein–Barr virus (EBV)[35] and herpes simplex virus 1 (HSV-1)[36] in previous experimental studies. The andrographolide derivative, i.e., 14-deoxy-11,12-didehydroandrographolide is one of the major components/derivatives of A. paniculata reported for its antiviral properties.[37-39]

Fig. 1

Structures of andrographolide and its derivative.

Recently, andrographolide has been investigated as a potential inhibitor of SARS-CoV-2 main protease (3CLpro) using an in-silico approach.[40] However, its potency to act over papain-like protease (PLpro) and spike protein has not been investigated yet. Further, there are numerous reports wherein various in silico approaches[41] such as molecular docking, fast pulling of ligand (FPL), free energy perturbation (FEP),[42] density functional theory (DFT),[43] high throughput virtual screening,[44] and drug repurposing studies[45] have been exploited to investigate various target proteins of SARS-COV-2. Likewise, there are recent reports of in-silico investigations of murine natural products[46] and some diverse scaffolds of synthetic compounds designed through in silico insights.[47,48] Since the risk of getting an infection with COVID-19 is reported to be higher in the subjects with compromised immunity,[10] it is important to consider in manipulating the immune system in them.

Hence, the present study aimed to investigate the prospective potential of andrographolide and one of the major derivatives i.e. 14-deoxy-11,12-didehydroandrographolide as a potent anti-viral agent by targeting three proteins of COVID-19, i.e., 3CLpro, PLpro, and spike protein. Further, the study also evaluated the plausible pathways to be regulated in enhancing the immune system.

Materials and methods

Prediction of targets

SMILES of 14-deoxy-11,12-didehydroandrographolide and andrographolide was retrieved from the PubChem (https://pubchem.ncbi.nlm.nih.gov/) database[49] and queried for protein-based prediction in DIGEP-Pred[50] at a pharmacological activity (Pa) > pharmacological inactivity (Pi).

Enrichment analysis

The list of up-and-down-regulated proteins was queried in the STRING database.[51] The biological process, cell component, and molecular function were recorded. The modulated protein and their associated pathways were identified using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. The interaction between the compounds, their targets, and pathways was constructed using Cytoscape (https://cytoscape.org/) and was analyzed using “edge count” of respective node.

In silico molecular docking

3D structures of 14-deoxy-11,12-didehydroandrographolide and andrographolide were retrieved from the PubChem database in .sdf format and converted into .pdb format using Discovery studio 2019. The ligand of each molecule was minimized using mmff94 forcefield and converted into ‘.pdbqt’ format. Structures of 3CLpro (PDB: 6LU7) and PLpro (PDB: 4M0W) were retrieved from the RCSB database,[52] which were complexed with water molecules and hetero atoms; removed using Discovery studio 2019 and saved in ‘.pdb’ format. The spike protein of coronavirus was homology modeled target using accession number AVP78042.1 as query sequence and PDB: 6VSB as a template using SWISS-MODEL.[53] Docking was carried using Autodock 4.0.[54] After docking ten different ligand conformations were obtained in which ligand possessing minimum binding energy was chosen to visualize the ligand–protein interaction using Discovery studio 2019.

Molecular dynamics (MD) simulation

MD simulations were carried out with the AMBER18 software package.[55] The ligands 14-deoxy-11,12-didehydroandrographolide and andrographolide were parameterized with ANTECHAMBER[56] employing GAFF force field. The amino acid residues of each protein under study were parameterized with the FF14SB force field. The xLEAP program was used to prepare the protein–ligand complexes of the target proteins 3CLpro, PLpro, and spike protein with ligands' docked poses. Each protein–ligand complex was solvated in a truncated octahedron of the TIP3P box. Appropriate counter ions Na+ and Cl− were added to neutralize the system. Thus prepared, the protein–ligand complexes were subjected to 100 ns MD simulations on Nvidia V100-SXM2-16GB GPU using the PMEMD.CUDA module. Initially, the system was subjected to energy minimization, where water and the complete system were minimized in two steps. Simulated annealing optimization, and the NVT and NPT equilibration steps of 5 ns, each was performed to equilibrate the system. The production phase MD simulations were performed at 1 atm constant pressure using Monte Carlo barostat and 300 K constant temperature by using Langevin thermostat. During the simulation, a collision frequency of 2 ps−1 and the volume exchange was attempted for every 100 fs. Hydrogen bonds were constrained by using the SHAKE algorithm, and the integration step of two fs was employed. Particle Mesh Ewald (PME) method was used to compute the long-range electrostatic interactions, while the cut-off of 8 Å was used to compute the short-range interactions. The program CPPTRAJ was used to analyze the interactions at every four ps on the result from the full trajectory. The MD simulation results were analyzed in terms of RMSD and RMSF of protein–ligand complexes.

Calculation of drug-likeness and ADMET

The compound's drug-likeness was calculated based on the rule of five using Molsoft[57] by querying the SMILES of compounds. Further, absorption, distribution, metabolism, excretion and toxicity (ADMET) profile were calculated using admetSAR2.0.[58]

Results

Andrographolide was predicted to regulate 36 proteins, of which 17 were down-regulated, and 19 were upregulated. Likewise, 14-deoxy-11,12-didehydroandrographolide regulated 48 proteins in which 21 were downregulated, and 27 were upregulated. The list of regulated proteins with their Pa and Pi of both compounds is summarized in Table 1.

Regulated proteins by andrographolide and 14-deoxy-11,12-didehydroandrographolidea

Andrographolide						14-Deoxy-11,12-didehydroandrographolide
DownRegulation			UpRegulation			DownRegulation			UpRegulation
Pa	Pi	Modulated proteins	Pa	Pi	Modulated proteins	Pa	Pi	Modulated proteins	Pa	Pi	Modulated proteins
0.548	0.132	TOP2A	0.589	0.131	VDR	0.66	0.117	CHEK1	0.701	0.077	VDR
0.559	0.163	CHEK1	0.526	0.082	CD14	0.627	0.097	TOP2A	0.599	0.043	CD14
0.387	0.028	KRT16	0.336	0.079	CLU	0.562	0.041	IVL	0.47	0.045	CLU
0.364	0.038	KRT17	0.444	0.198	AR	0.451	0.016	KRT16	0.548	0.124	AR
0.331	0.026	PTH	0.417	0.175	ID1	0.445	0.025	KRT17	0.522	0.161	CD83
0.394	0.138	ESR2	0.231	0.043	RAP1A	0.379	0.015	PTH	0.459	0.131	ID1
0.297	0.076	TIMP2	0.375	0.197	RAC1	0.452	0.092	ESR2	0.483	0.178	NPPB
0.38	0.161	CCL2	0.32	0.17	GPX1	0.474	0.127	MDM2	0.429	0.152	RAC1
0.38	0.161	IVL	0.241	0.092	KLK2	0.444	0.097	CCL2	0.41	0.134	SMN2
0.306	0.142	LEP	0.394	0.257	NPPB	0.34	0.044	TIMP2	0.39	0.124	TNFRSF1A
0.364	0.212	PRKCA	0.319	0.186	TNFRSF1A	0.357	0.086	LEP	0.262	0.027	RAP1A
0.3	0.155	CCL4	0.31	0.191	KRT18	0.403	0.139	PRKCA	0.362	0.129	KRT18
0.27	0.141	IL6R	0.171	0.058	RXRA	0.328	0.107	CCL4	0.325	0.096	CTSB
0.226	0.138	GYPA	0.356	0.248	RARA	0.382	0.187	NR3C1	0.35	0.124	GPX1
0.349	0.269	MDM2	0.205	0.103	RHOB	0.299	0.105	IL6R	0.287	0.077	KLK2
0.319	0.28	NR3C1	0.178	0.095	RHOA	0.384	0.206	CASP8	0.412	0.209	PLAT
0.221	0.206	CD44	0.354	0.328	CD83	0.139	0.051	PTHLH	0.397	0.211	RARA
			0.294	0.291	SMN2	0.247	0.166	CD44	0.371	0.201	CYP3A4
			0.191	0.19	CD38	0.362	0.293	NOS2	0.322	0.162	FKBP5
						0.276	0.229	FLT1	0.224	0.065	RHOB
						0.255	0.231	PROS1	0.195	0.04	RXRA
									0.2	0.057	RHOA
									0.231	0.094	CD38
									0.304	0.283	CAT
									0.243	0.226	PGR
									0.268	0.254	PLAU
									0.123	0.12	KRT7

Pa: pharmacological activity, Pi: pharmacological inactivity.

A total of seventy-two different pathways were identified to be regulated by the andrographolide, among which, pathways in cancer were primarily regulated by modulating nine genes, i.e., AR, ESR2, IL6R, MDM2, PRKCA, RAC1, RARA, RHOA, RXRA at the false discovery rate of 4.96 × 10−5. Similarly, 14-deoxy-11,12-didehydroandrographolide was predicted to regulate seventy-seven different pathways by modulating the Estrogen signaling pathway via seven genes, i.e., ESR2, FKBP5, KRT16, KRT17, KRT18, PGR, RARA at the false discovery rate of 7.57 × 10−6. Pathways modulated by andrographolide and 14-deoxy-11,12-didehydroandrographolide with their respective genes are summarized in Tables 2 and 3, respectively.

Enrichment analysis of andrographolide regulated targets

#Term ID	Term description	Observed gene count	False discovery rate	Matching proteins in the network (labels)
hsa05200	Pathways in cancer	9	4.96 × 10−5	AR, ESR2, IL6R, MDM2, PRKCA, RAC1, RARA, RHOA, RXRA
hsa04640	Hematopoietic cell lineage	5	7.61 × 10−5	CD14, CD38, CD44, GYPA, IL6R
hsa04972	Pancreatic secretion	5	7.61 × 10−5	CD38, PRKCA, RAC1, RAP1A, RHOA
hsa05130	Pathogenic Escherichia coli infection	4	0.00014	CD14, KRT18, PRKCA, RHOA
hsa04915	Estrogen signaling pathway	5	0.00016	ESR2, KRT16, KRT17, KRT18, RARA
hsa04932	Non-alcoholic fatty liver disease (NAFLD)	5	0.00022	IL6R, LEP, RAC1, RXRA, TNFRSF1A
hsa04062	Chemokine signaling pathway	5	0.00048	CCL2, CCL4, RAC1, RAP1A, RHOA
hsa05205	Proteoglycans in cancer	5	0.00059	CD44, MDM2, PRKCA, RAC1, RHOA
hsa04015	Rap1 signaling pathway	5	0.00063	ID1, PRKCA, RAC1, RAP1A, RHOA
hsa04670	Leukocyte transendothelial migration	4	0.00091	PRKCA, RAC1, RAP1A, RHOA
hsa04071	Sphingolipid signaling pathway	4	0.00095	PRKCA, RAC1, RHOA, TNFRSF1A
hsa04060	Cytokine–cytokine receptor interaction	5	0.0013	CCL2, CCL4, IL6R, LEP, TNFRSF1A
hsa04961	Endocrine and other factor-regulated calcium reabsorption	3	0.0013	KLK2, PRKCA, VDR
hsa05014	Amyotrophic lateral sclerosis (ALS)	3	0.0013	GPX1, RAC1, TNFRSF1A
hsa05418	Fluid shear stress and atherosclerosis	4	0.0013	CCL2, RAC1, RHOA, TNFRSF1A
hsa05206	MicroRNAs in cancer	4	0.0017	CD44, MDM2, PRKCA, RHOA
hsa04010	MAPK signaling pathway	5	0.0018	CD14, PRKCA, RAC1, RAP1A, TNFRSF1A
hsa04920	Adipocytokine signaling pathway	3	0.0024	LEP, RXRA, TNFRSF1A
hsa05152	Tuberculosis	4	0.0024	CD14, RHOA, TNFRSF1A, VDR
hsa05202	Transcriptional misregulation in cancer	4	0.0024	CD14, MDM2, RARA, RXRA
hsa05203	Viral carcinogenesis	4	0.0027	CHEK1, MDM2, RAC1, RHOA
hsa04151	PI3K-Akt signaling pathway	5	0.0031	IL6R, MDM2, PRKCA, RAC1, RXRA
hsa04510	Focal adhesion	4	0.0033	PRKCA, RAC1, RAP1A, RHOA
hsa05132	Salmonella infection	3	0.0035	CCL4, CD14, RAC1
hsa04064	NF-kappa B signaling pathway	3	0.0045	CCL4, CD14, TNFRSF1A
hsa04014	Ras signaling pathway	4	0.005	PRKCA, RAC1, RAP1A, RHOA
hsa04933	AGE-RAGE signaling pathway in diabetic complications	3	0.005	CCL2, PRKCA, RAC1
hsa04620	Toll-like receptor signaling pathway	3	0.0052	CCL4, CD14, RAC1
hsa04659	Th17 cell differentiation	3	0.0052	IL6R, RARA, RXRA
hsa04722	Neurotrophin signaling pathway	3	0.0067	RAC1, RAP1A, RHOA
hsa04919	Thyroid hormone signaling pathway	3	0.0067	MDM2, PRKCA, RXRA
hsa04310	Wnt signaling pathway	3	0.0116	PRKCA, RAC1, RHOA
hsa04150	mTOR signaling pathway	3	0.0124	PRKCA, RHOA, TNFRSF1A
hsa04921	Oxytocin signaling pathway	3	0.0124	CD38, PRKCA, RHOA
hsa04530	Tight junction	3	0.0161	RAC1, RAP1A, RHOA
hsa05144	Malaria	2	0.0161	CCL2, GYPA
hsa05164	Influenza A	3	0.0161	CCL2, PRKCA, TNFRSF1A
hsa04360	Axon guidance	3	0.0166	PRKCA, RAC1, RHOA
hsa04024	cAMP signaling pathway	3	0.0221	RAC1, RAP1A, RHOA
hsa04370	VEGF signaling pathway	2	0.0221	PRKCA, RAC1
hsa05169	Epstein–Barr virus infection	3	0.0221	CD38, CD44, MDM2
hsa04720	Long-term potentiation	2	0.0233	PRKCA, RAP1A
hsa04810	Regulation of actin cytoskeleton	3	0.0233	CD14, RAC1, RHOA
hsa05131	Shigellosis	2	0.0233	CD44, RAC1
hsa01524	Platinum drug resistance	2	0.0237	MDM2, TOP2A
hsa04115	p53 signaling pathway	2	0.0237	CHEK1, MDM2
hsa04520	Adherens junction	2	0.0237	RAC1, RHOA
hsa04664	Fc epsilon RI signaling pathway	2	0.0237	PRKCA, RAC1
hsa05100	Bacterial invasion of epithelial cells	2	0.0237	RAC1, RHOA
hsa05211	Renal cell carcinoma	2	0.0237	RAC1, RAP1A
hsa05214	Glioma	2	0.0237	MDM2, PRKCA
hsa05221	Acute myeloid leukemia	2	0.0237	CD14, RARA
hsa05223	Non-small cell lung cancer	2	0.0237	PRKCA, RXRA
hsa04918	Thyroid hormone synthesis	2	0.0239	GPX1, PRKCA
hsa05133	Pertussis	2	0.024	CD14, RHOA
hsa01521	EGFR tyrosine kinase inhibitor resistance	2	0.026	IL6R, PRKCA
hsa04350	TGF-beta signaling pathway	2	0.0287	ID1, RHOA
hsa05210	Colorectal cancer	2	0.0295	RAC1, RHOA
hsa04970	Salivary secretion	2	0.0297	CD38, PRKCA
hsa04666	Fc gamma R-mediated phagocytosis	2	0.0311	PRKCA, RAC1
hsa05146	Amoebiasis	2	0.0339	CD14, PRKCA
hsa01522	Endocrine resistance	2	0.034	ESR2, MDM2
hsa05215	Prostate cancer	2	0.0348	AR, MDM2
hsa04066	HIF-1 signaling pathway	2	0.0349	IL6R, PRKCA
hsa05231	Choline metabolism in cancer	2	0.0349	PRKCA, RAC1
hsa05142	Chagas disease (American trypanosomiasis)	2	0.0358	CCL2, TNFRSF1A
hsa04668	TNF signaling pathway	2	0.0399	CCL2, TNFRSF1A
hsa04270	Vascular smooth muscle contraction	2	0.047	PRKCA, RHOA
hsa04110	Cell cycle	2	0.0492	CHEK1, MDM2
hsa04380	Osteoclast differentiation	2	0.0492	RAC1, TNFRSF1A
hsa04611	Platelet activation	2	0.0492	RAP1A, RHOA
hsa04650	Natural killer cell mediated cytotoxicity	2	0.0492	PRKCA, RAC1

Enrichment analysis of 14-deoxy-11,12-didehydroandrographolide regulated targets

#Term ID	Term description	Observed gene count	False discovery rate	Matching proteins in the network
hsa04915	Estrogen signaling pathway	7	7.57 × 10−6	ESR2, FKBP5, KRT16, KRT17, KRT18, PGR, RARA
hsa05200	Pathways in cancer	11	7.57 × 10−6	AR, CASP8, ESR2, IL6R, MDM2, NOS2, PRKCA, RAC1, RARA, RHOA, RXRA
hsa05202	Transcriptional misregulation in cancer	7	1.32 × 10−5	CD14, FLT1, MDM2, PLAT, PLAU, RARA, RXRA
hsa04932	Non-alcoholic fatty liver disease (NAFLD)	6	9.13 × 10−5	CASP8, IL6R, LEP, RAC1, RXRA, TNFRSF1A
hsa04972	Pancreatic secretion	5	0.00016	CD38, PRKCA, RAC1, RAP1A, RHOA
hsa05152	Tuberculosis	6	0.00016	CASP8, CD14, NOS2, RHOA, TNFRSF1A, VDR
hsa04015	Rap1 signaling pathway	6	0.00023	FLT1, ID1, PRKCA, RAC1, RAP1A, RHOA
hsa05014	Amyotrophic lateral sclerosis (ALS)	4	0.00023	CAT, GPX1, RAC1, TNFRSF1A
hsa05130	Pathogenic Escherichia coli infection	4	0.00023	CD14, KRT18, PRKCA, RHOA
hsa05205	Proteoglycans in cancer	6	0.00023	CD44, MDM2, PLAU, PRKCA, RAC1, RHOA
hsa05418	Fluid shear stress and atherosclerosis	5	0.00035	CCL2, PLAT, RAC1, RHOA, TNFRSF1A
hsa05206	MicroRNAs in cancer	5	0.00054	CD44, MDM2, PLAU, PRKCA, RHOA
hsa04060	Cytokine–cytokine receptor interaction	6	0.00065	CCL2, CCL4, FLT1, IL6R, LEP, TNFRSF1A
hsa04610	Complement and coagulation cascades	4	0.00065	CLU, PLAT, PLAU, PROS1
hsa05132	Salmonella infection	4	0.00075	CCL4, CD14, NOS2, RAC1
hsa04010	MAPK signaling pathway	6	0.00095	CD14, FLT1, PRKCA, RAC1, RAP1A, TNFRSF1A
hsa04062	Chemokine signaling pathway	5	0.00095	CCL2, CCL4, RAC1, RAP1A, RHOA
hsa04064	NF-kappa B signaling pathway	4	0.00095	CCL4, CD14, PLAU, TNFRSF1A
hsa04066	HIF-1 signaling pathway	4	0.00095	FLT1, IL6R, NOS2, PRKCA
hsa04640	Hematopoietic cell lineage	4	0.00095	CD14, CD38, CD44, IL6R
hsa05203	Viral carcinogenesis	5	0.00095	CASP8, CHEK1, MDM2, RAC1, RHOA
hsa05215	Prostate cancer	4	0.00095	AR, MDM2, PLAT, PLAU
hsa04510	Focal adhesion	5	0.00097	FLT1, PRKCA, RAC1, RAP1A, RHOA
hsa04620	Toll-like receptor signaling pathway	4	0.00097	CASP8, CCL4, CD14, RAC1
hsa05142	Chagas disease (American trypanosomiasis)	4	0.00097	CASP8, CCL2, NOS2, TNFRSF1A
hsa04670	Leukocyte transendothelial migration	4	0.0013	PRKCA, RAC1, RAP1A, RHOA
hsa04071	Sphingolipid signaling pathway	4	0.0014	PRKCA, RAC1, RHOA, TNFRSF1A
hsa04151	PI3K-Akt signaling pathway	6	0.0014	FLT1, IL6R, MDM2, PRKCA, RAC1, RXRA
hsa04014	Ras signaling pathway	5	0.0015	FLT1, PRKCA, RAC1, RAP1A, RHOA
hsa04961	Endocrine and other factor-regulated calcium reabsorption	3	0.0015	KLK2, PRKCA, VDR
hsa04115	p53 signaling pathway	3	0.004	CASP8, CHEK1, MDM2
hsa04920	Adipocytokine signaling pathway	3	0.004	LEP, RXRA, TNFRSF1A
hsa01524	Platinum drug resistance	3	0.0041	CASP8, MDM2, TOP2A
hsa04621	NOD-like receptor signaling pathway	4	0.0041	CASP8, CCL2, CTSB, RHOA
hsa05133	Pertussis	3	0.0045	CD14, NOS2, RHOA
hsa05146	Amoebiasis	3	0.0085	CD14, NOS2, PRKCA
hsa04933	AGE-RAGE signaling pathway in diabetic complications	3	0.0093	CCL2, PRKCA, RAC1
hsa04659	Th17 cell differentiation	3	0.0101	IL6R, RARA, RXRA
hsa04668	TNF signaling pathway	3	0.0115	CASP8, CCL2, TNFRSF1A
hsa05145	Toxoplasmosis	3	0.0115	CASP8, NOS2, TNFRSF1A
hsa04215	Apoptosis – multiple species	2	0.0127	CASP8, TNFRSF1A
hsa04722	Neurotrophin signaling pathway	3	0.0127	RAC1, RAP1A, RHOA
hsa04919	Thyroid hormone signaling pathway	3	0.0127	MDM2, PRKCA, RXRA
hsa04210	Apoptosis	3	0.0188	CASP8, CTSB, TNFRSF1A
hsa04310	Wnt signaling pathway	3	0.0215	PRKCA, RAC1, RHOA
hsa04150	mTOR signaling pathway	3	0.0231	PRKCA, RHOA, TNFRSF1A
hsa04921	Oxytocin signaling pathway	3	0.0231	CD38, PRKCA, RHOA
hsa04530	Tight junction	3	0.0306	RAC1, RAP1A, RHOA
hsa05134	Legionellosis	2	0.0306	CASP8, CD14
hsa05164	Influenza A	3	0.0306	CCL2, PRKCA, TNFRSF1A
hsa05416	Viral myocarditis	2	0.0308	CASP8, RAC1
hsa04360	Axon guidance	3	0.0311	PRKCA, RAC1, RHOA
hsa04370	VEGF signaling pathway	2	0.0327	PRKCA, RAC1
hsa04020	Calcium signaling pathway	3	0.0328	CD38, NOS2, PRKCA
hsa05168	Herpes simplex infection	3	0.0332	CASP8, CCL2, TNFRSF1A
hsa05167	Kaposi's sarcoma-associated herpesvirus infection	3	0.0335	CASP8, RAC1, TNFRSF1A
hsa05131	Shigellosis	2	0.0343	CD44, RAC1
hsa04720	Long-term potentiation	2	0.0347	PRKCA, RAP1A
hsa04024	cAMP signaling pathway	3	0.0362	RAC1, RAP1A, RHOA
hsa04664	Fc epsilon RI signaling pathway	2	0.0362	PRKCA, RAC1
hsa05169	Epstein–Barr virus infection	3	0.0362	CD38, CD44, MDM2
hsa05211	Renal cell carcinoma	2	0.0362	RAC1, RAP1A
hsa05214	Glioma	2	0.0362	MDM2, PRKCA
hsa05221	Acute myeloid leukemia	2	0.0362	CD14, RARA
hsa05223	Non-small cell lung cancer	2	0.0362	PRKCA, RXRA
hsa04520	Adherens junction	2	0.037	RAC1, RHOA
hsa04810	Regulation of actin cytoskeleton	3	0.037	CD14, RAC1, RHOA
hsa04918	Thyroid hormone synthesis	2	0.037	GPX1, PRKCA
hsa04976	Bile secretion	2	0.037	CYP3A4, RXRA
hsa05100	Bacterial invasion of epithelial cells	2	0.037	RAC1, RHOA
hsa01521	EGFR tyrosine kinase inhibitor resistance	2	0.0409	IL6R, PRKCA
hsa04146	Peroxisome	2	0.0433	CAT, NOS2
hsa04350	TGF-beta signaling pathway	2	0.0446	ID1, RHOA
hsa05323	Rheumatoid arthritis	2	0.045	CCL2, FLT1
hsa05210	Colorectal cancer	2	0.0454	RAC1, RHOA
hsa04970	Salivary secretion	2	0.0457	CD38, PRKCA
hsa04666	Fc gamma R-mediated phagocytosis	2	0.0481	PRKCA, RAC1

Similarly, the interaction of both compounds with the proteins and regulated pathways is represented in Fig. 2 and 3.

Fig. 2

Interaction of andrographolide with the proteins and regulated pathways.

Fig. 3

Interaction of 14-deoxy-11,12-didehydroandrographolide with the proteins and regulated pathways.

Further, the number of genes in multiple cellular components, biological process, and molecular function for andrographolide and 14-deoxy-11,12-didehydroandrographolide are represented in Fig. 4 and5, respectively. Similarly, network analysis of 14-deoxy-11,12-didehydroandrographolide identified prime regulation of PRKCA protein and estrogen signaling pathway. Further, andrographolide primarily modulated PRKCA protein and pathways in cancer.

Fig. 4

GO enrichment analysis for andrographolide. (a) cellular component, (b) molecular function, and (c) biological process.

Fig. 5

GO enrichment analysis for 14-deoxy-11,12-didehydroandrographolide. (a) Cellular component, (b) molecular function, and (c) biological process.

In-silico molecular docking

Among andrographolide and 14-deoxy-11,12-didehydroandrographolide, 14-deoxy-11,12-didehydroandrographolide was predicted to have the highest binding affinity with PLpro, i.e., −6.7 kcal mol−1; however, it did not have any hydrogen bond interactions. Similarly, andrographolide showed −6.5 kcal mol−1 binding energy with PLpro with 1 hydrogen bond interaction, i.e., Tyr274. Although both molecules had equal binding energy with 3CLpro (−6.8 kcal mol−1), the number of hydrogen bond interactions were more in andrographolide due to interaction with Thr190, His163, and Cys145. Further, both molecules showed a binding affinity with spike protein, i.e., 6.9 kcal mol−1; however, andrographolide showed 1 hydrogen bond interaction with Lys807 (Table 4).

Docking hits of andrographolide and 14-deoxy-11,12-didehydroandrographolide with PLpro, 3CLpro, and spike protein

Targets	Ligand	Binding affinity (kcal mol−1)	Number of hydrogen bonds	Hydrogen bond residues
PLpro (PDB 4M0W)	14-Deoxy-11,12-didehydroandrographolide	−6.7	—	—
	Andrographolide	−6.5	1	Tyr274
3CLpro (PDB 6LU7)	14-Deoxy-11,12-didehydroandrographolide	−6.8	1	Arg131
	Andrographolide	−6.8	3	Thr190, His163, Cys145
Spike protein	14-Deoxy-11,12-didehydroandrographolide	−6.9	—	—
	Andrographolide	−6.9	1	Lys807

The interaction of each compound with the respective proteins is represented in Fig. 6.

Fig. 6

Docked poses of ligands at the binding site of each protein. Hydrogen bond donor and acceptor regions are shown as the surface around the binding site. (1) Binding site of papain-like protease (PLpro) (PDB:4M0W), (2) binding site of 3CLpro (PDB:6LU7), and (3) spike protein (panel a represents 14-deoxy-11,12-didehydroandrographolide and panel b represents andrographolide).

Molecular dynamics (MD) simulations

The flexibility at the binding site and the desolvation mechanism is not considered in the rigid docking methodology. However, molecular dynamics simulations can provide deeper insights into the interaction between ligand and protein amino acid residues at the atomistic level. The integrated workflow comprising molecular docking and molecular dynamics simulations is more suited in such situations as docking provides the most favourable bioactive poses of inhibitor molecules. In contrast, MD provides the insights of interactions and energetics in a biological environment.[59,60] The extended time scale MD simulations allow exploring a vast space of conformational optimization and its stability. In the present work, 100 ns MD simulation of well-equilibrated systems was performed on the complexes of 3CLpro, PLpro, and modeled spike protein, each of which is bound to 14-deoxy-11,12-didehydroandrographolide and andrographolide, respectively. The analysis of resulting trajectories comprising of 10 000 frames provides insights into the binding modes of inhibitor molecules, the formation of hydrogen bonds, pi–pi interactions, van der Waals interactions, and the consequent stability of the system in terms of RMSD, RMSF, and ligand-RMSD.

The MD trajectories of PLpro with 14-deoxy-11,12-didehydroandrographolide and andrographolide were analyzed for the protein RMSD, ligand RMSD and per residue fluctuations as RMSF (Fig. 7a–c). The RMSD analyses of PLpro bound with 14-deoxy-11,12-didehydroandrographolide showed a slight gradual increase in RMSD with initial equilibration at around 2 Å for the first 25 ns, after which it had gradually risen to 2.75 Å between 25 to 100 ns. These RMSD values point out the system's stability and strong binding affinity between the PLpro and the 14-deoxy-11,12-didehydroandrographolide molecule. The subtle but gradual increase in RMSD could be attributed to the binding site adaptation supported by the RMSF for the binding site residues aa150 to aa200 and aa220 to aa240. Interestingly, a similar trend in the RMSF was also observed in the complex with andrographolide. These residues are present at the binding cavity, and possibly they adopt the conformation suitable for both the ligands. These RMSD fluctuations in the case of PLpro and andrographolide complex were observed, reaching a maximum RMSD of 3 Å at around 85 ns. After that, they were gradually decreasing to around 2.5 Å towards the end of the simulation. The fluctuations in RMSD may be in part due to the C11–C12 rotatable bond, which may give rise to better conformational flexibility in the andrographolide molecule. The Lig-RMSD of 14-deoxy-11,12-didehydroandrographolide remains stable at 2.5 Å for about 75 ns and after that fluctuates and sharply rises to 15 Å. However, the Lig-RMSD of andrographolide remains stable at RMSD of 2.75 Å until 50 ns and further rises to a stable RMSD of 5 Å until the end of the simulation. The MD trajectories were visually inspected to investigate the fluctuations in Lig-RMSD (Fig. 7d and e). Both the phytochemicals adopt a conformationally more stable position by binding at the shallow binding cavity. Possibly because the hydroxyl group at the 14th position in andrographolide allows it to adopt a conformationally stable form throughout the simulation; however, the lack of this hydroxyl group and restricted rotation around the C11–C12 bond in 14-deoxy-11,12-didehydroandrographolide may be responsible for larger fluctuations in Lig-RMSD after 75 ns. The residues Val164 and Tyr274 participate in hydrogen bond formation with a carbonyl oxygen atom at C16 position in both the ligands. These bonds break in 14-deoxy-11,12-didehydroandrographolide more often after 75 ns due to restricted rotation around the C11–C12 bond. The ligand superposition also shows some structural conformational changes in both the ligands. The MDS studies on the PLpro–ligand complexes suggest that these complexes are relatively stable.

Fig. 7

Trajectory analysis for PLpro (PDB: 4M0W) bound to 14-deoxy-11,12-didehydroandrographolide and andrographolide; (a) root mean square deviation (RMSD), (b) root mean square deviation for each ligand (Lig-RMSD), and (c) root mean square fluctuations per amino acid (aa) (RMSF). Interaction analysis of the PLpro bound to ligands during the molecular dynamics simulation; (d) equilibrated structure of 14-deoxy-11,12-didehydroandrographolide bound to the PLpro before MDS production phase (green) and post-MDS production phase (red); (e) equilibrated structure of andrographolide bound to the PLpro before MDS production phase (green) and post-MDS production phase (red).

The MDS trajectories of 3CLpro bound to each ligand were analyzed, and the protein RMSD, ligand RMSD and per amino acid residue fluctuations, RMSF were recorded (Fig. 8a–c). The protein RMSD for the 14-deoxy-11,12-didehydroandrographolide rise to about 3 Å during the first 25 ns and then equilibrated at around 3 Å until 75 ns, after which it fluctuates and rose to around 4.5 Å towards the end of MDS. On the other hand, the andrographolide bound to 3CLpro equilibrates quickly, and the RMSD remains in 2 to 3 Å throughout the simulation. It indicates the fair stability of 3CLpro bound to andrographolide. This finding is also clearly reflected in the RMSF and ligand RMSD. The ligand RMSD of andrographolide is fairly constant to 2 Å until 40 ns of the MDS, increasing slightly to 2.5 Å after that and equilibrated in the same conformation for the rest of the MDS. 14-deoxy-11,12-didehydroandrographolide equilibrate initially with RMSD of 2 Å till 25 ns, and there is steep conformational change resulting in RMSD of 2.5 Å till around 70 ns.

Fig. 8

Trajectory analysis for 3CLpro (PDB: 6LU7) bound to 14-deoxy-11,12-didehydroandrographolide and andrographolide; (a) root mean square deviation (RMSD), (b) root mean square deviation for each ligand (Lig-RMSD), and (c) root mean square fluctuations per amino acid (aa) (RMSF). Interaction analysis of the 3CLpro bound to ligands during the molecular dynamics simulation; (d) equilibrated structure of 14-deoxy-11,12-didehydroandrographolide bound to the 3CLpro before MDS production phase (green) and post-MDS production phase (red); (e) equilibrated structure of andrographolide bound to the 3CLpro before MDS production phase (green) and post-MDS production phase (red).

Larger fluctuations in the RMSD were observed after 70 ns with an increase in RMSD to an average of around 7.5 Å. The results of ligand RMSD indicates better conformational stability of andrographolide than 14-deoxy-11,12-didehydroandrographolide (Fig. 8b). The per residue RMSF for both the complexes has a similar pattern of fluctuating residues involvement with the fluctuations ranging between 0.5 to 4 Å; however, the RMSF values for 14-deoxy-11,12-didehydroandrographolide are slightly higher than andrographolide. The residues aa48–aa52 and aa150–aa200 clearly show larger deviations in RMSF with 14-deoxy-11,12-didehydroandrographolide (Fig. 8c). A visual analysis of the MDS trajectories was performed to ascertain these observations, as shown in Fig. 8d and e. In the initial conformation of 14-deoxy-11,12-didehydroandrographolide bound to 3CLpro before MDS, a hydrogen bond between the carbonyl oxygen at C14 and Arg131 residue of the active site was observed. However, this hydrogen bond breaks and new hydrogen bonds were formed with other residues such as Gln109 and Thr190. Probably due to conformationally restricted bond rotation around C11–C12, these hydrogen bonds are formed less frequently, which is evident in RMSF values in these residues and ligand RMSD (Fig. 8d). In the case of andrographolide bound to 3CLpro, the initial conformation has three hydrogen bonds between C16-carbonyl oxygen–His163, C14-hydroxyl group oxygen–Cys145, and C19-hydroxyl oxygen–Thr190. During MDS's progress, some of these hydrogen bonds break, and new hydrogen bonds were formed with adjacent residues such as Ala191 and His164. However, due to conformational flexibility in andrographolide around the C11–C12 bond, the ligand stabilizes and quickly gains an energetically lower conformation (Fig. 8e). These observations suggest conformationally better stabilization of the andrographolide at the binding site of 3CLpro.

The MDS trajectories of modelled spike protein bound to both the ligands were analyzed. The protein RMSD, ligand RMSD and per amino acid residue fluctuations, and RMSF were recorded (Fig. 9a–c). The RMSD in spike protein bound to each ligand shows fluctuations in the range 5 to 12 Å, and it is acceptable due to the amino acid composition of the protein comprising more than 600 residues. In spike protein bound to 14-deoxy-11,12-didehydroandrographolide, an initial increase in RMSD to 10 Å till 25 ns simulation was observed, which remained stable with minor deviations, thereafter till the end of simulation with RMSD of 10 Å. This suggests the conformational stability of 14-deoxy-11,12-didehydroandrographolide at the binding cavity, which resulted in the system stability.

Fig. 9

Trajectory analysis for modeled spike protein bound to 14-deoxy-11,12-didehydroandrographolide and andrographolide; (a) root mean square deviation (RMSD), (b) root mean square deviation for each ligand (Lig-RMSD), and (c) root mean square fluctuations per amino acid (aa) (RMSF). Interaction analysis of the 3CLpro bound to ligands during the molecular dynamics simulation; (d) equilibrated structure of 14-deoxy-11,12-didehydroandrographolide bound to the spike protein before MDS production phase (green) and post-MDS production phase (red); (e) equilibrated structure of andrographolide bound to the spike protein before MDS production phase (green) and post-MDS production phase (red).

On the other hand, spike protein bound with andrographolide showed a similar trend in RMSD initially till 25 ns, which rises to around 12.5 Å during 25 ns to 100 ns. The binding site residues undergo conformational change during this simulation period. Possibly, the conformational change in the residues is due to conformational flexibility in the andrographolide molecule. The ligand RMSD and per residue RMSF supports this observation for andrographolide (Fig. 9b and c). The ligand RMSD of andrographolide increases sharply during initial MDS to around 15 Å until 25 ns and decreases to around 5 Å until 50 ns. However, it is unable to converge to a stable RMSD after that which suggests the major conformational changes in andrographolide and, consequently, the conformational changes in the binding site.

In contrast, the RMSD fluctuations in 14-deoxy-11,12-didehydroandrographolide are very subtle, with an initial rise to around 5 Å, and remain stable at this RMSD with minor deviations through the rest of the simulation period, suggests a stable complex and strong binding between the protein and ligand. The RMSF in spike protein residues also supports these observations. The residues aa300–aa550 clearly shows larger deviations in RMSF of around 7 to 12 Å with andrographolide (Fig. 9b). Most of these residues belong to the binding cavity. The corresponding RMSF values in the case of 14-deoxy-11,12-didehydroandrographolide for these residues range from 5 to 7 Å. A visual analysis of the MDS trajectories was also performed to ascertain these observations. The initial equilibrated conformation of 14-deoxy-11,12-didehydroandrographolide bound to spike protein has a hydrogen bond between the carbonyl oxygen at C14 and Tyr585 and C19-hydroxyl group oxygen and Pro263. However, these hydrogen bond breaks and a new hydrogen bond were formed between the C14 carbonyl oxygen and Trp582. In the case of andrographolide bound to spike protein, the initial equilibrated conformation shows a hydrogen bond between C3 hydroxyl group hydrogen and Asn285 and Gln282 residues. Due to conformational flexibility in andrographolide, during the production phase of MDS, these hydrogen bonds break. However, no new hydrogen bond formation was observed towards the end of the simulation.

Drug-likeness and ADMET profiling

14-Deoxy-11,12-didehydroandrographolide scored higher drug-likeness score, i.e., −0.52 compared to andrographolide, which was computed based on molecular weight, number of hydrogen bond donor, number of hydrogen bond acceptor, and log P value (Table 5). It has influenced both compounds' pharmacokinetic characters by affecting absorption, distribution, metabolism, excretion, and toxicity (Fig. 10).

Druglikeness score of andrographolide and 14-deoxy-11,12-didehydroandrographolide

		Andrographolide	14-Deoxy-11,12-didehydroandrographolide
Molecular formula		C20H30O5	C20H28O4
Molecular weight		350.21	332.20
Number of HBA		5	4
Number of HBD		3	2
Mollog P		2.19	3.09
Mollog S	log (mol L−1)	−1.97	−2.67
	mg L−1	3791.12	702.11
MolPSA (A2)		71.27	55.16
MolVol (A3)		416.03	421.79
Number of stereocenters		6	5
Drug-likeness model score		−0.64	−0.52

Fig. 10

ADMET profile of (A) andrographolide, and (B) 14-deoxy-11,12-didehydroandrographolide.

Discussion

The present study investigated one of the active biomolecules, andrographolide, and its derivative i.e., 14-deoxy-11,12-didehydroandrographolide from Andrographis paniculate, for the regulation of the proteins/immunomodulatory pathways and also assess their binding affinity with three targets, i.e., 3CLpro, PLpro, and spike protein involved in the COVID infection. Further, we investigated the drug-likeness property of both molecules in which andrographolide scored a lower drug-likeness character compared to 14-deoxy-11,12-didehydroandrographolide. However, on looking at the binding affinity and number of hydrogen bond interactions, andrographolide showed a higher interaction towards the selected targets. It suggests fewer modifications could be made in the andrographolide moiety to enhance its drug-likeness property without altering the binding affinity towards the targeted proteins. Further, some modifications could be made in 14-deoxy-11,12-didehydroandrographolide to amplify the hydrogen bond interaction and eventually increase the higher binding affinity with targeted proteins.

Subjects with lower immunity system are more prone to infection with COVID 19 due to compromised immunity system,[61] which is well proven in subjects suffering from an infectious and non-infectious disease(s). In this case, it is crucial to enhance the subjects' immunity to minimize the probability of viral infection. In the present study, via the enrichment analysis, we identified multiple pathways involved in boosting the immune system, which is modulated by andrographolide and 14-deoxy-11,12-didehydroandrographolide.

In the present study, we identified potential modulation of few pathways, directly or indirectly linked with the modulation of the immune system, i.e., chemokine signaling pathway, Rap1 signaling pathway, cytokine–cytokine receptor interaction, MAPK signaling pathway, NF-kappa B signaling pathway, RAS signaling pathway, p53 signaling pathway, HIF-1 signaling pathway, and natural killer cell-mediated cytotoxicity. Among the above pathways, chemokine signaling pathways, Rap1 signaling pathway, and cytokine–cytokine receptor interaction are the choice of interest pathways as they are directly linked with the immune system's regulation, and they scored minimum false discovery rate compared to the rest of the pathways. The chemokine signaling pathway was modulated by andrographolide and 14-deoxy-11,12-didehydroandrographolide, which could control the migration of immune cells in tissues.[62] Further, the Rap1 signaling pathway is involved in activating three secondary messengers, i.e., cAMP, calcium, and diacylglycerol,[63] which are needed in the signaling of cell position during viral infections; modulated by andrographolide by regulating ID1, PRKCA, RAC1, RAP1A, and RHOA and by 14-deoxy-11,12-didehydroandrographolide by regulating FLT1, ID1, PRKCA, RAC1, RAP1A, and RHOA. Similarly, the KEGG database has recorded cytokine–cytokine receptor interaction as an entry (hsa04060) in various auto-immune disorders. Since COVID-19 has a more risk over the infections on the altered immune system of subjects, modulation of this pathway could be beneficial in them, modulated by andrographolide and 14-deoxy-11,12-didehydroandrographolide. Further, the MAPK signaling pathway has been identified to play an essential role in the functioning of T lymphocytes,[64] was observed to be modulated by andrographolide and 14-deoxy-11,12-didehydroandrographolide. Additionally, other pathways like NF-kappa B signaling pathway, ras signaling pathway, p53 signaling pathway, HIF-1 signaling pathway, and natural killer cell-mediated cytotoxicity are also regulated which has been well reported to be involved in the modulation of the immune system.

In COVID-19 infection, the n-CoV-2 binds to ACE-2, enters into the cell, and starts deregulating the intracellular functions by altering the normal homeostatic stimulus.[65] Hence, it is needed to control the components by binding over them or responding towards the stimulus 3CLpro, or at least to minimize its effect by controlling the intracellular cascade initiated by the viral infection. Further, gene ontology enrichment analysis identified andrographolide and 14-deoxy-11,12-didehydroandrographolide to target the intracellular components, binding capacity towards various proteins as a molecular function and responder towards stimulus which could be the possible action of these two agents over the viral infection.

A concept of modulation of multiple proteins by a single molecule is the choice of research interest in identifying the lead hit towards respective targets. Further, andrographolide has been previously reported to possess anti-viral properties.[29] Hence, based on the same concept, andrographolide and 14-deoxy-11,12-didehydroandrographolide may also possess the anti-viral efficacy over COVID-19, which kindled us evaluating the binding affinity of these bioactives over PLpro, 3CLpro, and spike protein. Although the drug-likeness score model predicted 14-deoxy-11,12-didehydroandrographolide to behave like a drug based on “Rule of Five”, the binding affinity and number of hydrogen bond interactions reflected andrographolide to act more on three proteins of COVID-19 i.e. PLpro, 3CLpro, and spike protein.

Conclusion

The present study utilized the system biology approach to investigate the andrographolide and 14-deoxy-11,12-didehydroandrographolide against COVID-19 by modulating the multiple pathways in which the chemokine signaling pathway could be a choice of interest as it is directly linked in modulating the immune response and scored the lowest false discovery rate. Further, andrographolide could possess higher importance than 14-deoxy-11,12-didehydroandrographolide as it scored higher interaction with the targeted proteins of COVID-19. However, the present findings are completely based on computer simulations and database query; the outcome may vary based on processing units and database updates; suggests the necessity to confirm the present findings via well-designed experimental protocols and is the future scope of present findings.

Author contributions

MMW and BMP: conceptualization, supervision, investigation; RVC and RBP: molecular docking and dynamics studies; PK and YDD: network pharmacology and analysis; SNG and BV: methodology, software data analysis; SSG and RBP: writing-original draft preparation; MMW and SSG: formal analysis, visualization, reviewing and editing.

Ethical statement

This study doesn't include any animal or human study.

Funding support

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of interest

The authors declare no competing interests.