Dual targeting of cytokine storm and viral replication in COVID-19 by plant-derived steroidal pregnanes: An in silico perspective

Abbreviations

COVID-19

The Coronavirus disease 2019

SARS-CoV-2

Severe Acute Respiratory Syndrome Corona Virus 2

sRdRp

SARS-CoV-2RNA-dependent RNA polymerase

GR

Glucocorticoid receptors

hGRag

Human glucocorticoid receptors in agonist conformation

hGRagt

Human glucocorticoid receptors in atagonist conformation

JKs

Janus kinases

3CLpro

3-chymotrypsin-like protease

PLpro

Papain-like protease

ARDS

Acute Respiratory Distress Syndrome

PDPs

Plant derived pregnanes

IL

Interleukins

IFN

Interferons

1. Introduction

Coronavirus disease 2019 (COVID-19) is a clinical syndrome, caused by Severe Acute Respiratory Syndrome Corona Virus 2 (SARS-CoV-2) [1]. The clinical presentation of SARS-CoV-2 infections ranges from asymptomatic condition or mild symptoms (such as fever, cough, and generalized malaise) in the majority of the cases to severe respiratory failure. The early stage of infection, progresses to interstitial pneumonia and acute respiratory distress syndrome (ARDS) in nearly 10–20% of the cases, especially in the elderlies and people with co-morbidities [2]. The pathophysiology of SARS-CoV-2 infection is a complex mechanism that is known to mobilize several biomolecules of the immune and hematologic systems [3].

Cytokines are a group of polypeptide signaling molecules responsible for regulating a large number of biological processes via cell surface receptors [4]. The term “cytokine storm”, a condition characterized by an exaggerated activation of the immune system was first associated with onset of the graft-versus-host disease [5] and later known to be involved in several viral infections [6]. The exaggerated cytokine release in response to viral infection, has emerged as one of the mechanisms leading to acute respiratory distress syndrome and multiple-organ failure in COVID- 19 [7]. In this regard, recent studies have shown that patients with COVID-19 have higher levels of inflammatory cytokines, such as interleukin (IL)-1β, IL-2, IL-6 IL-7, IL-8, IL-9, IL-10, IL-18, tumor necrosis factor (TNF)-α, granulocyte colony-stimulating factor (G-CSF), granulocyte-macrophage colony-stimulating factor, fibroblast growth factor, macrophage inflammatory protein 1, compared to healthy individuals [8]; circulating levels of IL-6, IL-10, and TNF-α also correlated with illness severity as they were significantly higher in intensive care unit (ICU) patients compared to mild/moderate cases. At this point, anti-viral treatment alone is not enough and should be combined with appropriate anti-inflammatory treatment. Anti-rheumatic drugs, which are tried for managing cytokine storm of SARS-CoV-2 infection include: corticosteroids, JAK inhibitors, IL-6 inhibitors, IL-1 inhibitors, anti-TNF-α agents, hydroxychloroquine, intravenous immunoglobulin (IVIG), and colchicines [9].

The interaction of glucocorticoids receptors (GR) and its ligands, glucocorticoids (GCs), have been explored for the modulation of cytokines in acute and chronic inflammatory diseases [10]. Glucocorticoids modulate cytokine expression by several genomic mechanisms. The activated GR complex: (i) binds to the promotor responsive elements, thereby inactivating key pro-inflammatory transcription factors (e.g. AP-1, NF kappa B); (ii) upregulates the expression of cytokine inhibitory proteins, e.g. I kappa B, which inactivates the transcription factor NF kappa B, thereby suppressing the secondary expression of a series of cytokines; and lastly, (iii) reduces the half-lives and utility of cytokine mRNAs [11]. Unfortunately, the use of GCs is limited by unwanted severe side effects, such as osteoporosis, disorders of glucose and lipid metabolism, and hypertension [12]. Therefore, there is an urgent need for the development of natural compounds with higher anti-inflammatory activity compared to standard GCs, alongside antiviral potential and with accompanied no or low toxicity [13,14]. Janus kinases (JAKs) mediate the signaling of numerous cytokines and growth factors involved in the regulation of inflammation, immunity and hematopoiesis [15]. Among the JAK family members, the JAK1 has the broadest cytokine signaling profile,being the only isoform that pairs with the other three JAKs. The pairing of JAK1with JAK3 regulates the signaling of the gamma common (γc) cytokines. The pairing of JAK1with JAK2 regulates the signaling of type I interferons (IFNα, IFNβ), type II interferon (IFNγ) and the IL-10 family of cytokines [16]. Inhibitors of the JAK-STAT pathway, such as baricitinib and Ruxolitnib, are used for suppressing proinflammatory cytokine production and systemic inflammation. Interleukin-6 (IL-6) is a pleiotropic cytokine. In general, IL-6 inhibitors prevent human IL-6 from binding to IL-6 receptors, thus impeding the formation of immune signaling complexes on cell surfaces [17].

Along with structural proteins, the SARS viral genomes encode non-structural proteins, including 3-chymotrypsin-like protease (3CL^pro), papain-like protease (PLpro), helicase and RNA-dependent RNA polymerase (RdRp) which are important target for the development of therapeutics [18]. The proteolytic processing of the polyproteins is performed by the viral cysteine proteases to yield 16 non-structural proteins; 3CL^pro cleaves and modifies the viral polyproteins at 11 sites while PL^pro cleaves the first three sites at the N-terminus [4,19]. The RNA-dependent RNA polymerase (RdRp), is a central component of coronaviral replication/transcription machinery that catalysis RNA-template dependent formation of phosphodiester bonds between ribonucleotides In our recent work, we have demonstrated the potential of some natural compounds as inhibitors of these proteins [[20], [21], [22]].

Recently, dexamethasone, a potent synthetic anti-inflammatory glucocorticoid was declared as the world's first treatment proven effective in reducing the risks of death through cytokine storm among severely ill COVID-19 patients based on clinical trial results [23,24]. Through computational and biological comparison, few plant-derived steroidal compounds have been suggested as modulators of inflammation through interactions with GR (Dean et al., 2017; Morsy et al., 2019). Such plant-derived anti-inflammatory steroids like glycyrrhetinic acid [25], guggulsterone [26], boswellic acid [27], withaferin A [28] and diosgenin [29] have a common cyclopentanoperhydrophenanthrene steroid ring structure. Pregnanes are naturally occurring C₂₁ steroidal compounds that have been documented with wide range of bioactivities including anti-inflammatory activity [22,[30], [31], [32]]. Due to the present COVID-19 pandemic, there is urgent need for such plant-derived steroids that may possess dual interference with cytokine storm and viral replicase/transcriptase complex but with fewer side effects. Thus, the aim of this study was to screen an in-house library of plant-derived steroidal pregnanes for hGR agonist using a comparative molecular docking approach.

2. Materials and methods

3. Results

3.2. Amino acid interaction of selected pregnanes with target proteins

The amino acid interactions of hGR in the agonist conformation, hIL-6, hJAK,s3CL^pro, sPL^pro and sRdRp with reference inhibitors and the top three ranked PDPs that demonstrated the highest binding tendencies are represented in Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7 . The interacting residues of the proteins with respective ligand groups were majorly through H-bond, hydrophobic interactions and few other bonds (Table 2 ). The revalidation of the docking pattern of the native ligand (dexamethasone) co-crystalized with hGRag showed that dexamethasone was docked into to ligand-binding domain (LBD) of hGRag. The A-ring of dexamethasone was positioned adjacent to the β-strands 1 and 2 while the D ring was close to helix 12 of hGRag. The 3′-carbonyl oxygen of the A ring formed a hydrogen bond to the guanidinium group of ARG⁶¹¹ of hGRag. On the C ring, the 11α-hydroxyl group formed a hydrogen bond with the carbonyl group of LEU⁵⁶³ while on the D ring, the 20-hyroxyl and 21-carbonyl groups formed a hydrogen bond with Thr739 and ASN⁵⁶⁴ of hGRag respectively. The 16α formed 2 alkyl bonds to TYR735 and LEU732, while the 18 and 19-methyl groups displayed an alkyl interaction with, CYS736 and MET604 of hGRag respectively (Fig. 2a). A Pi-alkyl interaction was observed between the D ring and the remaining amino acid residues of hGRag. Glaucogenin C, the topmost ranked PDP for hGRag was also docked into the LBD of hGRag. The 3-hydroxyl and 15, 20α-diepoxy groups of glaucogenin C interacted via hydrogen bonds with GLN⁵⁷⁰ and THR⁷³⁹, while the double bond between carbons 13 and 18 formed 2 hydrogen bonds with ASN⁵⁶⁴ and MET⁵⁶⁰. An alkyl interaction was observed between TRY⁷³⁵ and 19-methyl moiety, while several Pi-alkyl interactions were observed between the A, B and C rings and the remaining amino acid residues of hGRag (Fig. 2b). In a similar manner to the 3′-carbonyl oxygen of dexamethasone, the 2-hydroxyl group of hirundigenin formed a hydrogen bond with ARG⁶¹¹, while the other hydrogen bond was formed between ASN⁵⁶⁴ of hGRag and 20-oxahexacyclo group. Numerous alkyl interactions were formed by 5- and 19-methyl groups while the pi-alkyl interactions were formed by the rings (Fig. 2c). Bregenin was also docked into the LBD of hGRag, interacting with the amino acid residues of the active site. The 3-, 16- and 17-hydroxyl groups of bregenin formed 3 hydrogen bonds with GLN⁵⁷⁰, GLN⁶⁴² and LEU⁷³² of hGRag. The 10- and 13-methyl groups formed alkyl interactions with MET⁶⁰¹, CYS⁷³⁶ and MET⁶⁰⁴ of hGRag. The remaining residues interacted via Pi-alkyl interaction with the B and C rings of bregenin (Fig. 2d). l(+)-tartaric acid, the reference inhibitor, and the native molecule bound to the crystallographic structure of hIL-6 were docked into the “site 1” binding site. Five hydrogen bonds were observed between tartarate and IL-6. Direct hydrogen bonds to which ARG¹⁷⁹ and ARG¹⁸² of hIL-6 served as the donors of four pairs of hydrogen atoms were formed with α-carboxyl moiety, while the α-hydroxyl group of the tartarate donated the hydrogen atom for the hydrogen bond with GLN¹⁷⁵ (Fig. 3a). ARG³⁰ and GLN¹⁷⁵ of hIL-6 donated the hydrogen atoms for the hydrogen bonds formed with the carbonyl group of atratogenin A, while a carbon-hydrogen bond was formed between the furan ring and LEU³³(Fig. 3b).Alkyl interactions were formed between the 4β-methyl moiety and ARG³⁰ and LEU³⁰ while Pi-alkyl interactions were formed by the B and furan rings with LEU¹⁷⁸ and LUE³³ of hIL-6 respectively. Glaucogenin C was docked into the same binding site and interacted with some of the amino acid residues as methotrexate (Fig. 3c). A conventional hydrogen bond and carbon-hydrogen bond were formed with ARG¹⁸² and GLN¹⁷⁵ of hIL-6 respectively while most of the alkyl interactions were formed with 5- and 19- methyl groups. ASP³⁴ of hIL-6 donated a hydrogen atom to form hydrogen bond with 7-hydroxyl group, while the alkyl interactions were formed by the four rings of anhydroholantogenin (Fig. 3d). The amino group of the pyrimidine ring of roxolitinib (reference inhibitor) contributed two hydrogen atoms to form hydrogen bonds with GLU⁹⁵⁷ and GLY¹⁰²⁰ of hJAK1, while that of pyraxole ring formed two hydrogen bonds with ARG¹⁰⁰⁷ and ASN¹⁰⁰⁸ of hJAK1. Two Pi-sigma bonds were formed between the pyrrole ring of roxolitinib and LEU⁸⁸¹ and LEU¹⁰¹⁰. An alkyl interaction was observed between the cyclopentane ring of roxolitinib and ARG¹⁰⁰⁷ (Fig. 4a). The 6-hydroxyl group of atratogenin B donated the hydrogen atom for the only hydrogen bond formed with ARG¹⁰⁰³, while the 2- and 4-methyl groups and the methyl group attached to the furan ring interacted via alkyl interactions with VAL⁸⁸⁹, ALA⁹⁰⁶ and LEU¹⁰¹⁰ of hJAK1 respectively (Fig. 4b). Hirundigenin was docked into the same active site as roxolitinib; 8- and 16-hydroxyl groups and 20-oxahexacyclo ring of hirundigenin formed hydrogen bonds with LEU⁹⁵⁹ and ASN¹⁰⁰⁸, while the alkyl interactions were formed between the A and B rings and VAL⁸⁸⁹, LEU¹⁰¹⁰ and LEU⁸⁸¹ of hJAK1(Fig. 4c). For glaucogenin C, the 21-carbonyl group formed two hydrogen bonds and the 8-hydroxyl group formed a hydrogen bond with ARG¹⁰⁰⁷, SER⁹⁶³ and GLY¹⁰²⁰ of hJAK1 respectively. The alkyl interactions were majorly contributed by 5-methyl group of glaucogenin C (Fig. 4d). Ritonavir was docked into the receptor-binding site and interacted with amino acid residues that form the catalytic dyad (Cys-145 and His-41) of s3CL^pro via a conventional hydrogen bond to LEU¹⁴¹ while the remaining interactions with HIS¹⁶⁴, THR²⁴ and CYS¹⁴⁵ were via carbon-hydrogen bonds. It further interacted via Pi-alkyl, Pi-Pi T-Shaped and Pi-sulphur with LEU²⁷, HIS⁴¹ and MET⁴⁹ of s3CL^pro respectively (Fig. 5a). The three top ranked PDPs for s3CL^pro were docked in the same binding site as the reference compound (ritonavir). Glaucogenin D interacted via conventional hydrogen and carbon-hydrogen bonds with GLY¹⁴³ and HIS⁴¹ of s3CL^pro, while it interacted with MET¹⁶⁵ and THR²⁴ via alkyl interactions (Fig. 5b). A conventional hydrogen bond was formed between the 8- and 16-hydroxyl groups of hirundigenin and the catalytic residues (THR²⁴ and CYS¹⁴⁵) of s3CL^pro, while the 19-methyl group formed an alkyl interaction with MET¹⁶⁵ (Fig. 5c). The 7-hydroxyl group of anhydroholantogenin formed hydrogen bond with THR²⁴ of s3CL^pro while the remaining hydrophobic interactions were formed by 10- and 17-methyl groups and the rings (Fig. 5d). Disulfiram, a known inhibitor of PL^pro, was docked into the binding cavity of SARS-C0V-2 PL^pro. It interacted with the amino acids HIS²⁷² and TRP¹⁰⁶ via a pi-sulphur interaction in the binding cavity of sPL^pro (Fig. 6a). In the same vein, glaucogenin D, glaucogenin A and glaucogenin C were docked into the same binding site. The carbonyl and 19-methyl groups interacted via conventional hydrogen and Pi-alkyl interaction with HIS²⁷², while an alkyl interaction was observed with TRP¹⁰⁶ of sPL^pro (Fig. 6b–d). The 7- and 8-hydroxyl groups of glaucogenin A formed two hydrogen bonds with ASP²⁸⁶, 8-hydroxyl group interacted via carbon-hydrogen bond with HIS²⁷², 19-methyl group interacted via alkyl interaction with CYS²⁷⁰, while the pentacyclo ring formed multiple Pi-Sigma bonds with TRP¹⁰⁶ (Fig. 6c). Two hydrogen bonds were formed between the carbonyl group of glaucogenin C and HIS²⁷² and TRP¹⁰⁶, while multiple Pi-alkyl interactions were formed with the same amino acid residues of sPL^pro (Fig. 6d). The 4-aminopyrrolo[2,1-f] [1,2,4]triazin-7-yl, 5-cyno and carbonyl groups of sRdRp served as hydrogen donors for all of the conventional hydrogen bonds with the catalytic residues. The alkyl end of the 2-ethylbutyl moiety of sRdRp interacted via alkyl interaction with CYS⁶²², while an electrostatic force was formed between the phosphoryl group and APS⁶¹⁸ (Fig. 7a). In case of glaucogenin A, 7- and 8-hydroxyl groups of sRdRp formed two hydrogen bonds with GLUS⁸¹¹ and TRP⁶¹⁷ respectively, while a conventional hydrogen bond and an alkyl interaction were formed by the double bonds at positions 21 and 10 with ASP⁷⁶⁰ and LYS⁷⁹⁸ respectively (Fig. 7b). The 12-hyroxyl group of Glaucogenin D formed the only two hydrogen bonds with APS⁷⁶⁰ and TRY⁶¹⁵ of sRdRp while Glaucogenin D exhibited similar binding pattern with that of Glaucogenin A (Fig. 7 c & d) (see Table 3).

Fig. 2.

Open in a new tab

Fig. 3.

Open in a new tab

Fig. 4.

Open in a new tab

Fig. 5.

Open in a new tab

Fig. 6.

Open in a new tab

Fig. 7.

Open in a new tab

Table 2.

Combined list of the top three ranked plant-derived pregnanes with the lowest binding energies for each of the protein targets in human and SARS-COV 2^a.

S/No	Plant Pregnanes		Plant species
1	Bregenin		Sarcostemma brevistigma
2	Hirundigenin		Vincetoxicum officinale
3	Anhydroholantogenin		Holarrhena antidysenterica
4	Atratogenin A		Cyanchum Atratum
5	Atratogenin B		Cyanchum Atratum
6	Glaucogenin A		Cynanchum glaucesens hand-maazz
7	Glaucogenin C		Cynanchum glaucesens hand-maazz
8	Glaucogenin D		Cynanchum glaucesens hand-maazz

Open in a new tab

^a

Protein Targets: Human glucocorticoid receptors in the agonist conformation (hGRag); human Interleukin-6 (hIL-6), human Janus kinases (hJAK1); SARS-CoV-2 3-chymotrypsin-like protease (s3CL^pro); SARS-CoV-2 papain-like protease (sPL^pro) and SARS-CoV-2 RNA-dependent RNA polymerase (sRdRp).

Table 3.

Interaction of reference inhibitors and plant derived pregnanes with amino acid residue of various targets^a.

Compounds	Protein	Numbers	Hydrogen bonds (Bond distance Å)	Numbers	Hydrophobic Interaction (Bond distance Å)	Numbers	Other interactions (Bond distance Å)
			Interacting residues		Interacting residues		Interacting residues
Dexamethasone	hGRag	2	ASN564 (2.14) THR739 (3.07) LEU563 (3.07) ARG611 (2.21)	8	TYR735 (4.46) MET601 (5.44) LEU732 (4.68) CYS736 (4.22) MET646 (5.04) MET604 (5.14)		none
Glaucogenin C		4	THR739 (3.07) MET560 (3.07) ASN564 (3.07) GLN570 (3.07)	7	TYR735 (3.07) PHE623 (3.07) LEU732 (3.07) LEU608 (3.07) MET601 (3.07) MET646 (3.07) LEU563 (3.07)		none
Hirundigenin		2	ASN564 (1.64) ARG611 (3.17)	7	TYR735 (5.30) PHE623 (4.72) LEU6085.41) MET604 (5.23) MET646 (5.33) LEU563 (4.21)MET601 (5.23)		none
Bregenin		4	GLN570 (2.03) GLN642 (2.24) LEU732 (2.57)	10	CYS736 (3.71) PHE623 (4.72, 4.81) LEU6085.42) MET604 (5.23, 5.20) MET646 (5.35) LEU563 (4.88, 3.73) MET601 (5.26)	2	LYS336 (2) (4.08, 3.82)
Methotrexate	hIL-6	5	ARG179 (2.76, 3.04) ARG182 (2.78, 2.74) GLN175 (2.80)				none
Atratogenin A		5	GLN175 (2.01) ARG30 (2.79) LUE33 (3.68) LUE178 (2.89)	4	LUE178 (4.67) ARG30 (4.20) LUE33 (4.28, 4.87)		none
Glaucogenin C		8	GLN175 (3.29) ARG182 (2.44)	5	ARG30 (3.92) LUE33 (4.05, 3.82) LUE178 (5.27) ARG179 (4.40)	2
Anhydroholantogenin		1	GLN175 (2.78)	3	LUE178 (5.22, 4.50) ARG179 (5.49)
Ruxolitnib	hJAK1	4	ARG1007 (2.93) ASN1008 (2.16) GLU957 (2.69) GLY1020 (2.23)	4	LEU881 (3.77) LEU1010 (3.82) ALA906 (4.80) ARG1007 (4.50)
Atratogenin B		1	ARG1003 (2.85)	9	LEU881 (4.63, 4.63) LEU1010 (3.46, 4.81, 4.81) ALA906 (4.38) ARG1007 (4.50) VAL889 (3.88, 5.15)
Hirundigenin		3	ARG1007 (1.95, 3.77) LEU959 (2.65) ASN1008 (2.37)	5	LEU881 (5.30, 4.09) LEU1010 (3.46, 4.28) VAL889 (4.38)
Glaucogenin C		4	ARG1007 (3.02) SER963 (2.81) GLY882 (3.66) GLY1020 (2.17)	5	LEU881 (5.36, 4.88) LEU1010 (4.46) VAL889 (4.32, 3.17)
Ritonavir	s3CLpro	4	CYS145 (3.84) HIS164 (3.07) LEU141 (2.88) THR24 (3.67)	3	LEU27 (5.41) HIS41 (5.24) MET49 (4.82)
Glaucogenin D		2	GLY143 (2.74) HIS41 (3.36)	2	MET165 (5.07) THR24 (3.67)
Hirundigenin		2	CYS145(3.57) THR24 (2.50)	1	MET165 (4.26)
Anhydroholantogenin		2	ASN142(3.67) THR24 (2.70)	7	CYS145(5.47) HIS41 (3.36, 5.36, 4.14) MET49 (4.43) MET165 (4.33,5.46)
Disulfiram	sPLpro			3	HIS272 (3.84, 4.01) TRP106 (1.06)
Glaucogenin D		1	HIS272 (2.82)	5	TRP106 (4.21, 5.05, 4.30) HIS272 (3.54,5.09)
Glaucogenin A		3	HIS272 (3.75) ASP286 (2.20,2.35)	7	TRP106 (3.86, 3.81, 4.44, 4.60, 4.95) TRP270 (4.27) HIS272 (5.31)
Glaucogenin C		2	HIS272 (2.76) TRP106 (2.72)	6	HIS272 (3.54, 5.08) TRP106 (4.23, 5.28, 4.30, 4.70)
Remdesivir	sRdRp	5	TYR619 (2.11) ASP761 (2.72) CYS813 (2.79) SER814 (2.13) ASP760 (2.17)	2	CYS622 (3.94) ASP613 (5.27)
Glaucogenin A		3	GLU811 (2.48) TRP617 (3.45) ASP760 (2.17)	1	LYS798 (4.95)
Glaucogenin D		2	TRP619 (3.45) ASP760 (2.17)		noon
Glaucogenin C		2	GLU811 (2.18) ASP760 (3.36)	1	LYS798 (3.98)

Open in a new tab

^a

Targets: Human glucocorticoid receptors in the agonist conformation (hGRag); human Interleukin-6 (hIL-6), human Janus kinases (hJAK1); SARS-CoV-2 3-chymotrypsin-like protease (s3CL^pro); SARS-CoV-2 papain-like protease (sPL^pro) and SARS-CoV-2 RNA-dependent RNA polymerase (sRdRp).

4. Discussion

Parallel advances in protein crystallography and various virtual-screening software for the modeling of ligand–receptor interactions have enhanced computer-aided drug design [49]. In this study a structure based virtual-screening of PDPs was employed via competitive docking approach for hGR agonist with a dual inhibitory potential against cytokine storm syndrome and viral replication in COVID-19. The top potential agonists were further analyzed for multiplicity of inhibitory tendencies against the hIL-6 and hJAK1(used as anti-proinflammatory targets), and s3CL^pro, sPL^pro and sRdRp (used as SARS-CoV-2 therapeutic targets). The docking of the PDPs to the hGRs identified the top ranked twenty PDPs with dexamethasone binding mode to hGR agonist. A total of eight PDPs were selected. From these eight PDPs, top three ranked PDPs (glaucogenin C, hirundigenin and bregenin) for hGR agonist were competitively and selectively docked to the hGRag. They were docked into the hydrophobic ligand binding pocket (LBP) which is located in the bottom half of the GR ligand binding domain, LBD [50,51]. The polar residues on the LBP interacted with the dexamethasone and the top ranked PDPs via several hydrogen bonds [51]. The binding of the amino acid residues on helix 12 and the loop preceding helix 12 have been earlier hypnotized to stabilize the helix in the active conformation that could serve as the molecular basis for the ligand-dependent activation of GR [52]. Among the several amino acids involved in the interactions, Cys-736 has been implicated in the interactions with heat shock proteins [52], Tyr-735, has been shown to be important for transactivation [53], while Gln-642 have been reported to play a unique role in steroid recognition [52]. In addition to the steroid structures of glucocorticoids, the 3′-carbonyl oxygen, 2′-carbonyl oxygen, double bonds between C4 and C5, 17ʹ hydroxyl group and 21ʹ hydroxyl group, are critical for anti-inflammatory potency and glucocorticoid receptor affinity [54]. The identified PDPs contained similar and analogous functional groups that interacted with GR; thus, the binding of these plant steroidal pregnanes may initiate the ligand-dependent activation of GR since they share similar binding patterns with dexamethasone. The activated glucocorticoid‐receptor complex can: (i) bind the promotor responsive elements (RE) of key pro-inflammatory transcription factors (e.g. AP‐1, NF kappa B) to inactivate them; (ii) upregulate the expression of cytokine inhibitory proteins, e.g. I kappa B, via glucocorticoid RE; and (iii) reduce the half‐life time and usefulness of cytokine mRNAs [11]. IL-6 is a major causative factor of inflammatory disease and it is a promising target, as well as its signaling pathways; however, orally available small-molecule drugs specific for IL-6 have not been developed [55]. From the PDPs with high binding tendencies to hGRag, three PDPs (atratogeninA, glaucogenin C and anhydroholantogenin) exhibited the lowest binding energy poses for hIL-6 in the same binding site as observed for the co-crystallized ligands (tartaric acid) of hIL-6. In a similar study, furosemide exhibiting the same binding mode as tartaric acid was further found to inhibit hIL-6 activity in vitro [56]. From the X-ray crystal diffraction of hIL-6 structure, it was shown to contain four alpha helices (helices A, B, C, and D), which were linked with loops. The receptor-binding domain is located at the C-terminus (residues 175–181) [57], in which ARG¹⁷⁹ is known to be the key residue [57]. AB loop and helices A and D is important in receptor binding and signal transduction [58]. Compounds that interact strongly with residue ARG¹⁷⁹ may interfere with the binding of the receptor to its ligands [59], thus these PDPs may proffer anti-inflammatory activity via hIL-6 inhibition. The Janus kinases (JAK) family of proteins function as critical mediators of cytokine signaling from membrane receptors to various signal transducers and activators of transcription (STAT) family of proteins [60]. Activation of STATs by the JAK kinases promotes the transcriptional activation of target genes controlling cell proliferation and survival, angiogenesis, and immune function [61]. Some JAK family inhibitors such as tofacitinib [62] and ruxolitinib [63] have progressed into clinical trials, and FDA approvals. In comparison with the reference inhibitor, ruxolitnib, the top-three ranked PDPs (atratogenin B, hirundigenin, glaucogenin C) with the best binding modes, for which hJAK1 had the highest affinities, interacted with the hinge residues LEU⁹⁵⁹, GLU⁹⁵⁷ and the side chain of ASN¹⁰⁰⁸ and the backbone carbonyl oxygen of ARG¹⁰⁰⁷ of the catalytic residues. These residues are involved in the inhibitory activities of selected compounds in both in silico and in vitro analyses [[64], [65], [66]].

The catalytic dyad (His⁴¹ and Cys¹⁴⁵) of SARS-CoV-2 3CL^pro is domiciled between its domain I (residues 8–101) and domain II (residues 102–184) [67]. A long loop (residues 185–200) that connects domain II and domain III (residues 201–303) completes the 3CLpro monomer [68]. The enzymatic activity of 3CLpro resides in its catalytic dyad of Cys¹⁴⁵ and His⁴¹ [69]. The catalytic dyad and the residue Glu¹⁶⁶ have been reported to be involved in the protein dimerization and substrate cleaving through the catalytic activities present in the cleft between domains Ⅰ and Ⅱ. The substrate-binding pocket lies in the cleft consisting of residue 140–145 and 163, 166 of domain II [70,71]. The substrate-binding pocket is divided into a series of subsites (S1–S6), each accommodating a single but consecutive amino acid residue in the substrate. The key residues in the substrate binding pockets of 3CLpro are His41, Tyr⁵⁴, Met⁴⁹, Phe¹⁴⁰-Cys¹⁴⁵, His¹⁶³−Pro¹⁶⁸, His¹⁷², and Asp¹⁸⁷− Gln¹⁹² regions [72]. In the same binding pattern as ritonavir (a known inhibitor of the protease), the top-three ranked pregnane (glaucogenin D, hirundigenin and anhydroholantogenin) were docked into the substrate binding pocket, interacting with various catalytic residues listed above. Considering our results and references to existing literature, the strong interaction of these PDPs to the critical residues (most especially HIS⁴¹ Gly¹⁴³, Cys¹⁴⁵ and MET¹⁶⁵) will greatly impair the dimerization and substrate binding of the SARS CoV-2 3CLpro.

The catalytic triad of SARS-CoV-2 PL^prois formed by CYS¹¹¹, HIS²⁷² and ASP²⁸⁶ [73,74], while TRP¹⁰⁶, GLY²⁵⁶, and LYS²⁷⁴ are catalytic residues (Li et al., 2020). LEU¹⁶², GLU¹⁶⁷, ASP¹⁶⁴ and TYR²⁶⁴ have been reported to be crucial for deubiquitinating activity of PLpro [75]. The host innate immune system is critical to controlling SARS-CoV-2 infection. Reverse post-translational modifications of immune proteins, such as interferon factor 3 and NF-κB via ubiquitination and the suppression of interferon-stimulated gene product 15 (ISG15), have also been implicated in the activities of PL^pro of SARS-CoV-2. [73,76], these, in turn, assist SARS-CoV-2 to escape the host innate immune responses. Pregnanes (glaucogenin D, hirundigenin and anhydroholantogenin) interacted with the catalytic triad and residues that are involved in deubiquitination. Such interactions may alter the catalytic conformation of PL^proand inhibit its ability to reverse ubiquitination. SARS-CoV-2 RdRp plays a central role in coronaviral replication/transcription machinery; it is, therefore, accepted as an excellent target for new therapeutics for which lead inhibitors, such as remdesivir, have been approved by the FDA [77]. Glaucogenin A, glaucogenin D and glaucogenin C were docked into the Motif C of the enzyme, exhibiting the same binding pattern as remdesivir. Motif C, the region comprising amino acid residues 753 to 767, contains the catalytic residues SER⁷⁵⁹, ASP⁷⁶⁰, and ASP⁷⁶¹ in the β-turn structure [77]. The stability of the complexes formed by the pregnanes with the enzyme stemmed from the vast number of interactions with the catalytic residues in the Motif C of the active site of the enzyme.

The several thermodynamics parameters (RMSD, RMSF, SASA and RoG) that were analyzed from the 100 ns full atomistic MDS trajectory files of the top two ranked PDPs-hGRag complexes revealed a high degree of stability throughout the period of the MDS run as compared to the apo receptor.

The RMSD plots showed that the binding of glaucogenin C and hirundigenin in the same manners as Dex to the active region of hGRag still preserved the structural integrity of the protein [78]. The.

The RMSF plots indicates the flexibility of different regions of a protein and the amino acid residue along the trajectory, which can be related to crystallographic B factors [78]. The RMSF plots of the PDPs-hGRag complexes shows similar plot pattern with the Dex-hGRag complex. Though a lower amount of fluctuation occurred at with the interacting residues, it has been established that greater amounts of structural fluctuations usually occur in regions known to be involved in ligand binding and catalysis, notably the catalytic loop regions [79] The RoG and SASA was assessed to evaluate the compactness and the accessibility of hGRag upon the binding of the PDPs and Dex. The PDPs and Dex maintained a reasonably steady RoG and SASA over the simulation time, indicating a highly compacted hGRag - PDPs and Dex complexes and well folded protein structure with intact intermolecular bonds [80]. The approximately close H-bonds between the top two ranked PDPs and Dex - hGRag complexes as compared to the unbound hGRag protein further indicated that the structural integrity of the protein was preserved in each of the system. At a quantitative level, simulation-based methods provide substantially more accurate estimates of ligand binding affinities (free energies) [81]. These results are calculated based on the total binding free energy of the complex. In these calculations, the binding free energy (ΔG_bind) measures the affinity of a ligand to its target protein. The free energy difference between the ligand-bound state (complex) and the corresponding unbound states of proteins and ligands are also employed in the calculations. Thus, the ΔG_bind calculations are important to gain in-depth knowledge about the binding modes of the hits in drug design [82]. The result from the MMPBSA calculation further corroborated the docking studies. Though Dex a known inhibitor to the hGRag protein presented the highest ΔG_bind, hirundigenin the top ranked PDP presented a very high and close ΔG_bind to Dex. A further evaluation of the MDS trajectories through clustering analysis showed that for each of the representative conformers from the several clusters, the interactions (H-bonds and hydrophobic interaction) were preserved at different time frames, indicating that the interactions can be maintained in a dynamic environment, thus can be well adapted for experimental procedures.

Despite the various efforts to improve current glucocorticoids and anti-inflammatory drugs, they still pose significant side effects [83], Hence the top-docked PDPs to various proteins were subjected to in silico physiochemical and ADMET analysis. The eight top-ranked PDPs fulfilled the all the requirements for the five physicochemical filtering analysis as reported by Lipinski [84] Ghose [85], Veber [86], Egan [87] and Muegge [88] thereby suggesting favourable physicochemical/druggable properties. The top-eight ranked PDPs expressed positive and high probability of human intestinal absorption and non-substrate but inhibitor of the permeability-glycoprotein (P-gp). These PDPs are predicted to be well absorbed into the blood stream subverting the capability of P-gp to pump them back into the intestinal lumen, bile ducts, urine-conducting ducts and capillaries respectively [89]. The blood brain barrier (BBB) penetration descriptor, predicts the ability of the PDPs to penetrate the blood brain barrier. The top-eight PDPs displayed the properties that suggested their ability to cross the BBB. SARS-CoV-2 has been reported to infect the brain, thus indicating its ability to cross the BBB [90], these PDPs may cross the BBB for to exert an overall viral clearance.

The top-eight PDPs displayed a probability of at least 65% ability to be bond to the plasma protein, suggesting their ability to be transported by these proteins. The estimated half-life time (less than 2 h) and clearance ratefall within the moderate range. The three phytochemicals presented a tolerable LD₅₀ between (51–500 mg/kg). The hERG channel plays a vital role in the repolarization and termination stages of action potential in cardiac cells [91]. Compounds that block the hERG channel may cause cardiotoxicity [92]. The top-eight PDPs exhibit low probability of being a potential hERG channel blockers, suggesting that they may not cause hERG channel-related cardiotoxicity [92]. Using the mutagenicity and skin sensitization descriptors, the top-eight PDPs did not display the properties to be mutagenic in silico, thereby suggesting that they may not cause genetic mutations, which do initiate the pathophysiology of other diseases. The impact of the PDPs on the liver phase I drug metabolism was also analyzed using the various cytochrome P450 descriptors. The top-eight PDPs demonstrated no inhibitory potential for the various cytochrome P450, thus may not adversely affect phase I drug metabolism in the liver. ADME/tox analysis indicated high aqueous solubility, ability to pass the high human intestinal absorption, low acute oral toxicity with a good bioavailability score. Therefore this natural plant pregnane may be considered to be non toxic with druggable potential.

5. Conclusion

In this study we employed a competitive docking approach to screen 117 plant derived pregnanes (PDPs) for hGR agonist, with a dual inhibitory potential against SARS-CoV-2 infection and the accompanied cytokine storm syndrome. Eight PDPs (bregenin, hirundigenin, anhydroholantogenin, atratogenin A, atratogenin B, glaucogenin A, glaucogenin C and glaucogenin D) with high agonist binding tendencies to the hGRag displayed different levels of multiplicity of inhibitory potentials to other pro-inflammatory targets (hIL-6, hJAK1) and three SARS-CoV-2 therapeutic targets (s3CL^pro, sPL^pro and sRdRp). The 8 PDPs fulfilled the requirements for various physicochemical and ADMET descriptors thereby suggesting favourable druggable properties. The top two ranked PDPs (glaucogenin C and hirundigenin) complexed to the hGRag demonstrated a high degree of structural stability and flexibility in a 100 ns simulated dynamics environment. These promising hGRag agonists with anti-inflammatory and SARS-CoV-2 replication inhibitory potential is recommended for further in vitro and in vivo experiments.

Declaration of competing interest

The authors declare no conflicting interest.

Biographies

Summary

The high morbidity and mortality rate of Severe Acute Respiratory Syndrome CoronaVirus 2 (SARS-CoV-2) infection arises majorly from the Acute Respiratory Distress Syndrome (ARDS) and “cytokine storm” syndrome, which is sustained by an aberrant systemic inflammatory response and elevated pro-inflammatory cytokines. Thus, the identification of compounds which target multiple proteins in the virus and exhibit anti-inflammatory activity will enhance the development of effective drugs against the disease. In this study, we carried out an in silico evaluation of some plant-derived pregnanes for their activities against selected human pro-inflammatory and SARS-CoV-2 replication proteins targets. This was carried out by a virtual screening of an in-house library of steroidal plant-derived pregnanes (PDPs). One hundred and six (106) PDPs were docked into the active regions of human glucocorticoid receptors (hGRs) in the agonist (hGRag) and antagonist (hGRagt) conformation, in a competitive molecular docking approach. Based on the minimal binding energy and a comparative dexamethason binding mode analysis, a hit-list of the top twenty ranked PDPs that were docked in the agonist conformation of hGR, with binding energies ranging between −9.8 and −11.2 kcal/mol, was defined. The top twenty ranked PDPS were further analyzed for interactions with the human Janus kinases 1 andInterleukins-6 (hJAK1 and hIL-6 respectively), and SARS-CoV-2 3-chymotrypsin-like protease, Papain-like protease and RNA-dependent RNA polymerase (3CLpro, PLpro and RdRP respectively). For each of the 6 targeted proteins (3 humans and 3 SARS CoV-2), the top three ranked PDPs were selected, to give a sum of eight PDPs (bregenin, hirundigenin, anhydroholantogenin, atratogenin A, atratogenin B, glaucogenin A, glaucogenin C and glaucogenin D) with multiplicity of high binding tendencies to the catalytic residues of different targets. From this eight PDPs, glaucogenin C and hirundigenin having the highest agonist tendencies to the hGR were further subjected to a 100 ns atomistic molecular dynamics simulation. A high degree of structural stability was observed from molecular dynamics simulation analyses of glaucogenin C and hirundigenin complexes of hGRag. A further clustering of the MDS trajectories of the complexes of glaucogenin C and hirundigenin with the hGRag shows that the interactions of these PDPS with the active site residues of hGRag were preserved in different representative structures of the clusters. The selected top-eight ranked PDPs demonstrated favourable druggable properties over the Lipinski, Veber, Ghose, Egan and Muegge predictive filters. In the same vein the 8 PDPs displayed favourable in silico ADMET properties over a wide range of predictive molecular descriptors, such as, ability to pass the blood brain barriers, high intestinal absorption, non-substrate to the permeability glycoprotein, non hERG blockers, non inhibitors of the cytochrome p450 etc. Thus, these promising hGRag agonists, especially glaucogenin C and hirundigenin, with potential anti-inflammatory and SARS-CoV-2 replication inhibitory activity is recommended for lead optimization for drug candidate and further evaluation in an in vitro and in vivo experiment.

Funding

The research did not receive any funding or grants.

Ethical

Approval: Not applicable.

Consent to participate

Not applicable.

Consent to publish

Not applicable.

Authors contributions

G. A. Gyebi Conceived and designed the analysis.

O. M. Ogunyemi Performed molecular docking analysis.

I. M. Ibrahim Performed molecular simulations.

J. O. Adebayo Editing and review of manuscript

S. O. Afolabi Editing and review of manuscript.

Availability of data and materials

The authors confirm that the data supporting the findings of this study are available within the article [and/or] its supplementary materials.

Appendix A. Supplementary data

The following is the Supplementary data to this article: