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

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 and “cytokine storm” syndrome, which is sustained by an aberrant systemic inammatory response and elevated pro-inammatory cytokines. Thus, phytocompounds with broad-spectrum anti-inammatory activity that target multiple SARS-CoV-2 proteins will enhance the development of effective drugs against the disease. In this study, an in-house library of 106 steriodal plant-derived pregnanes (PDPs) was docked in the active regions of human glucocorticoid receptors (hGRs) in a comparative molecular docking analysis. Based on the minimal binding energy and a comparative dexamethason binding mode analysis, a list of top twenty ranked PDPs docked in the agonist conformation of hGR, with binding energies ranging between -9.8 and -11.2 Kcal/mol, was obtained and analyzed for interactions with the human Janus kinases 1 and Interleukins-6 and SARS-CoV-2 3-chymotrypsin-like protease, Papain-like protease and RNA-dependent RNA polymerase. For each target protein, the top three ranked PDPs were selected. Eight PDPs (bregenin, hirundigenin, anhydroholantogenin, atratogenin A, atratogenin B, glaucogenin A, glaucogenin C and glaucogenin D) with high binding tendencies to the catalytic residues of multiple targets were identied. A high degree of structural stability was observed from the 100 ns molecular dynamics simulation analyses of glaucogenin C and hirundigenin complexes of hGR. The selected top-eight ranked PDPs demonstrated favourable druggable and in silico ADMET properties. Thus, the therapeutic potentials of glaucogenin C and hirundigenin can be explored for further in vitro and in vivo studies. 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 lters. In the same vein the 8 PDPs displayed favorable 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-inammatory 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.


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 those having older age and 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 rst 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 in ammatory 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, broblast growth factor, macrophage in ammatory protein 1, compared to healthy individuals [8]; circulating levels of IL-6, IL-10, and TNF-α also correlated with illness severity as they were signi cantly 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-in ammatory 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 in ammatory 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-in ammatory 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 the urgent need for the development of natural compounds with higher anti-in ammatory activity compared to standard GCs, alongside antiviral potential and low toxicity. Janus kinases (JAKs) mediate the signaling of numerous cytokines and growth factors involved in the regulation of in ammation, immunity and hematopoiesis [13]. Among the JAK family members, the JAK1 has the broadest cytokine signaling pro le,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 [14]. Inhibitors of the JAK-STAT pathway, such as baricitinib and Ruxolitnib, are used for suppressing proin ammatory cytokine production and systemic in ammation. 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 [15].
Along with structural proteins, the SARS viral genomes encode non-structural proteins, including 3chymotrypsin-like protease (3CL pro ), papain-like protease (PLpro), helicase and RNA-dependent RNA polymerase (RdRp) which are important target for the development of therapeutics [16]. The proteolytic processing of the polyproteins is performed by the viral cysteine proteases to yield 16 non-structural proteins; 3CL pro cleaves and modi es the viral polyproteins at 11 sites while PL pro cleaves the rst three sites at the N-terminus [4,17]. 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 [18][19][20].
Recently, dexamethasone, a potent synthetic anti-in ammatory glucocorticoid was declared as the world's rst treatment proven effective in reducing the risks of death through cytokine storm among severely ill COVID-19 patients based on clinical trial results [21,22]. Through computational and biological comparison, few plant-derived steroidal compounds have been suggested as modulators of in ammation through interactions with GR (Dean et al., 2017; Morsy et al., 2019). Such plant-derived anti-in ammatory steroids like glycyrrhetinic acid [23], guggulsterone [24], boswellic acid [25], withaferin A [26] and diosgenin [27] have a common cyclopentanoperhydrophenanthrene steroid ring structure. Pregnanes are naturally occurring C21 steroidal compounds that have been documented with wide range of bioactivities including anti-in ammatory activity [20,[28][29][30]. 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.1
Retrieval of protein structure The three-dimensional (3D) structure of human glucocorticoid receptors in the agonist conformation

Protein preparation
The crystal structures of the of proteins were processed by removing existing ligands and water molecules while missing hydrogen atoms were added according to the amino acid protonation state at pH 7.0 utilizing Autodock version 4.2 program (Scripps Research Institute, La Jolla, CA). Thereafter, nonpolar hydrogens were merged while polar hydrogens were added to each protein. The process was repeated for each protein and subsequently saved into a dockable pdbqt format for molecular docking.

Ligand preparation
PDPs ( [31,32] was employed for a structure based identi cation of agonist of the hGR protein. The approach combined separate molecular docking models for hGR in the agonist and antagonist conformations. True agonists and antagonists that were native ligands to the cocrystalized structures were extracted and rst used to evaluate the ability of this approach to differentiate agonists and antagonists. An initial docking analysis of the PDPs (103) and reference compounds (positive control: dexamethasone and negative control: mifepristone) to the hGRag (4UDC) was conducted. Ranking based on minimum binding energies and interactions with catalytic residues was employed to generate a list of the top twenty PDPs hits with the the highest binding tendencies. (Table  S1). A competitive docking analysis of the top twenty PDPs to another hGRagt (1NHZ) in the antagonist conformation (reference compounds: mifepristone as positive control and dexamethasone as negative control) was further conducted. Both hGRagand hGRagtwere co-crystallized with dexamethasone and mifepristone (the positive controls) respectively. The PDPs, reference compounds and the hGR proteins were loaded into PyRx (Python prescription) 0.8 [33] with the incorporation of Autodock vina [34]. For each of the docking steps, the ligands were imported via the OpenBabel [35], a plug-in tool in PyRx 0.8. The Universal Force Field (UFF) was used as the energy minimization parameter and conjugate gradient descent as the optimization algorithm. The binding site coordinates of the active site regions of the hGRag as de ned by the grid boxes were used for docking analysis (table 1a). All the other parameters were kept as default. After the completion of the docking process, the binding a nities of the protein for the compounds for the selected clusters were recorded. The compounds were then ranked by their binding energies.

2.4.1
Active site targeted molecular docking to other proteins targets Using the same protocol above, the top twenty PDPs with the lowest binding energies to the hGRag in the agonist conformation and the reference inhibitors were docked to the active region of ve proteins: human interleukin-6, human janus kinases, SARS-CoV-2 3-chymotrypsin-like protease, SARS-CoV-2 papain-like protease and SARS-CoV-2 RNA-dependent RNA polymerase as de ned by the grid boxes (table  1b). The molecular interactions of the top three PDPs with the highest binding a nities to each of the proteins and the reference compounds were viewed with Discovery Studio Visualizer version 16.

Molecular Dynamics Simulation
Molecular Dynamics Simulation (MDS) was performed on the hGR in the agonist conformation (apo) protein and top-two PDPs complexed with hGRag protein using NAMD software version 2.13 [36]. Necessary les for MDS were generated using CHARMM-GUI webserver [37,38]. For each complex or apo protein, the system was minimized for 10000 steps in constant number of atoms, constant volume and constant temperature (NVT) ensemble then a production run for 100 ns in NVT ensemble was performed. Temperature was set to be 310 K and salt concentration was set to be the physiological concentration 0.154 M NaCl. Afterwards, calculations of Backbone-Root Mean Square Deviation (RMSD), Per residue Root Mean Square Fluctuations (RMSF), Radius of Gyration (RoG), Surface Accessible Surface Area (SASA) were performed using VMD TK console scripts [39]. TTClust version 4.7.2 were used to cluster the trajectory automatically according to the elbow method, a representative structure for each cluster was produce [40]. These representative conformations were analyzed using Protein Ligand Interaction Pro ler (PLIP) for pregnane atom-amino acid residue interactive analysis [41]. The images were created using PyMol V2.2.2 [42] 2.6 ADMET study Eight compounds which are the top 3 compounds to the 6 protein targets were selected for evaluation of the drug-likeness and ADMET ltering analysis. The drug-likeness analysis which includes Lipinski, Veber, Ghose, Egan and Muegge were performed on the SwissADME (http://www.swissadme.ch/index.php) webserver. [43], while the predicted Absorption, Distribution, Metabolism, Excretion and toxicity (ADME/tox) study was analysed using the SuperPred webserver (http://lmmd.ecust.edu.cn/admetsar1/predict/) [44]. The SDF le and canonical SMILES of the compounds were downloaded from PubChem Database or copied from ChemDraw to calculate ADMET properties using default parameters.

Molecular docking
The binding a nities from the docking analysis of the proteins for the PDPs (103) and the reference compounds are shown in Table S1 (supplementary material). Based on the minimum binding energies and interactions with catalytic residues, the top twenty PDPs with binding energies ranging from -9.8 to -11.2 Kcal/mol for hGRag was compared to the binding energies of the controls (positive controldexamethasone = -12.2 Kcal/mol and negative control -mifepristone = -6.0 Kcal/mol; Figure 1). Using competitive docking approach,these results were compared with the results obtained from the docking analysis of the top twenty PDPs with the antagonist conformation of the same protein (hGRagt). The binding energies of the positive control (mifepristone: -11. 5 Kcal/mol), negative control (dexamethasone: -8.7 Kcal/mol)and the top twenty PDPs (between -7.7 and -8.8 Kcal/mol) are presented in Table S1 (supplementary material) and Figure 1. It was also observed that the binding a nities of hGRag for glaucogenin C, hirundigenin and bregenin (-11.2, -10.8 and -10.6 Kcal/mol respectively), the top three PDPs, were higher compared to those of hGRagt for them (-8.8, -8.7 and -8.7 Kcal/mol respectively).
From the interaction of the top twenty ranked PDPs with hGRag, hIL-6, hJAK1,s3CL pro , sPL pro and sRdRp, top three PDPs with the lowest binding energies for each of the proteins were obtained, yielding a combined list of eight pregnanes: bregenin, hirundigenin, anhydroholantogenin, atratogenin A, atratogenin B, glaucogenin A, glaucogenin C and glaucogenin D. From this list, glaucogenin C, hirundigenin, glaucogenin A and glaucogenin D were part of the top three ranked PDPs with least binding energies for at least two proteins, while glaucogenin C, with the least binding energy for hGRag, was listed among the top three ranked PDPs for hIL-6, hJAK1, sPL pro and sRdRp, thereby exhibiting multiplicity of binding properties. It was also observed that apart from GR and JAK, the three top pregnanes for each protein had binding energies for other proteins that were lower than those of the reference compounds.

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 Figures 2 to 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  (Figure 3a). ARG 30 and GLN 175 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 33 ( Figure   3b).Alkyl interactions were formed between the 4β-methyl moiety and ARG 30 and LEU 30 while Pi-alkyl interactions were formed by the B and furan rings with LEU 178 and LUE 33 of hIL-6 respectively. Glaucogenin C was docked into the same binding site and interacted with some of the amino acid residues as methotrexate (Figure 3c). A conventional hydrogen bond and carbon-hydrogen bond were formed with ARG 182 and GLN 175 of hIL-6 respectively while most of the alkyl interactions were formed with 5-and 19-methyl groups. ASP 34 of hIL-6 donated a hydrogen atom to form hydrogen bond with 7hydroxyl group, while the alkyl interactions were formed by the four rings of anhydroholantogenin ( Figure   3d). The amino group of the pyrimidine ring of roxolitinib (reference inhibitor) contributed two hydrogen atoms to form hydrogen bonds with GLU 957 and GLY 1020 of hJAK1, while that of pyraxole ring formed two hydrogen bonds with ARG 1007 and ASN 1008 of hJAK1. Two Pi-sigma bonds were formed between the pyrrole ring of roxolitinib and LEU 881 and LEU 1010 . An alkyl interaction was observed between the cyclopentane ring of roxolitinib and ARG 1007 (Figure 4a).

Results for Molecular dynamics
The stability, structural/conformational uctuations that occurred in the docked PDPs-hGRag systems were monitored in a simulated dynamic environment. The apo form and two complexes of hGRag with glaucogenin C and hirundigenin were used in a MDS study for 100 ns in NVT ensemble. The results were analyzed using VMD Tk console scripts to calculate RMSD, SASA, RoG, and RMSF. The mean RMSD of the apo form, hGRag -glaucogenin C, and hGRag -hirundigenin complexes are 1.58Å, 1.  Table 4 shows the results of the number of clusters that were generated and the interactions at different clusters, using a representative conformer. The most common interactions in both complexes are hydrophobic with few hydrogen bonds. The most amino acid that commonly participates in the interactions of hGRag -glaucogenin C complex is LEU 563 while in hGRag -Hirundigenin complex there are two amino acids, which are LEU 563 and LEU 566 . Figures 13 and 14 show the rst and last cluster of the hGRag -glaucogenin C complex while gures 15 and 16 show the rst and last clusters of protein-Hirundigenin complex.

Results for In Silico Drug-likeness and ADMET properties of top docked plant derived pregnanes to selected human and SARS-CoV-2 proteins
From the docking analysis, eight plant pregnanes (bregenin, hirundigenin, anhydroholantogenin, atratogenin A, atratogenin B, glaucogenin A, glaucogenin C and glaucogenin D) with high binding tendencies to hGRag with corresponding high binding tendencies to hIL-6, hJAK1, s3CL pro , sPL pro , and sRdRp were subjected to the predictive drug-likeness and ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) ltering analyses. The results for the predictive ltering analyses are presented in Table 5.

Discussion
Parallel advances in protein crystallography and various virtual-screening software for the modeling of ligand-receptor interactions have enhanced computer-aided drug design [45]. 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 analysed for multiplicity of inhibitory tendencies against the hIL-6 and hJAK1(used as anti-proin ammatory targets), and s3CL pro , sPL pro and sRdRp (used as SARS-CoV-2 therapeutic targets). The docking of the PDPs to the hGRs identi ed 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 [46,47]. The polar residues on the LBP interacted with the dexamethasone and the top ranked PDPs via several hydrogen bonds [47]. 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 [48]. Among the several amino acids involved in the interactions, Cys-736 has been implicated in the interactions with heat shock proteins [48], Tyr-735, has been shown to be important for transactivation [49], while Gln-642 have been reported to play a unique role in steroid recognition [48]. 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-in ammatory potency and glucocorticoid receptor a nity [50]. The identi ed 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-in ammatory 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 in ammatory disease and it is a promising target, as well as its signaling pathways; however, orally available small-molecule drugs speci c for IL-6 have not been developed [51]. 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. [52]. From the Xray 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) [53], in which ARG 179 is known to be the key residue [53]. AB loop and helices A and D is important in receptor binding and signal transduction [54]. Compounds that interact strongly with residue ARG 179 may interfere with the binding of the receptor to its ligands [55], thus these PDPs may proffer antiin ammatory 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 [56]. Activation of STATs by the JAK kinases promotes the transcriptional activation of target genes controlling cell proliferation and survival, angiogenesis, and immune function [57]. Some JAK family inhibitors such as tofacitinib [58] and ruxolitinib [59] 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 whichhJAK1 had the highest a nities, interacted with the hinge residues LEU 959 , GLU 957 and the side chain of ASN 1008 and the backbone carbonyl oxygen of ARG 1007 of the catalytic residues. These residues are involved in the inhibitory activities of selected compounds in both in silico and in vitro analyses [60][61][62].
The catalytic dyad (His 41 and Cys 145 ) of 3CL pro is domiciled between its domain I (residues 8-101) and domain II (residues 102-184) [63]. A long loop (residues 185-200) that connects domain II and domain III (residues 201-303) completes the 3CLpro monomer [64]. In the same binding pattern as ritonavir, the top-three ranked pregnane (glaucogenin D, hirundigenin and anhydroholantogenin) were docked into the dyad, interacting with various catalytic residues in respective domains I and II. The strong interactions with these amino acid residues may cause some conformational changes at the active site which may inhibit the catalytic activity of 3CL pro . The strong hydrogen bonding and hydrophobic interaction may suggest them as potential 3CL pro inhibitors. The catalytic triad of SARS-CoV-2 PL pro is formed by CYS 111 , HIS 272 and ASP 286 [65,66] , while TRP 106 , GLY 256 , and LYS 274 are catalytic residues . LEU 162 , GLU 167 , ASP 164 and TYR 264 have been reported to be crucial for deubiquitinating activity of PLpro [67]. The host innate immune system is critical to controlling SARS-CoV-2 infection. Reverse posttranslational modi cations 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. [65,68], 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 pro and 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 [69]. 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 759 , ASP 760 , and ASP 761 in the β-turn structure [69]. 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 les of the top two ranked PDP -hGRag complexes revealed a high degree of stability throughout the period of the MDS run as compared to the apo receptor. A further evaluation of the MDS trajectories through clustering analysis showed that for each of the representative conformers from 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-in ammatory drugs, they still pose signi cant side effects [70], Hence the top-docked PDPs to various proteins were subjected to in silico physiochemical and ADMET analysis. The eight top-ranked PDPs ful lled the all the requirements for the ve physicochemical ltering analysis as reported by Lipinski [71] Ghose [72], Veber [73], Egan [74] and Muegge [75] 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 [76]. 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 [77], 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 hours) and clearance ratefall within the moderate range. The three phytochemicals presented a tolerable LD 50 between (51~500 mg/kg). The hERG channel plays a vital role in the repolarization and termination stages of action potential in cardiac cells [78]. Compounds that block the hERG channel may cause cardiotoxicity [79]. The top-eight PDPs did not exhibit the potential of being hERG channel blockers, suggesting that they may not cause hERG channel-related cardiotoxicity. 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 analysed 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.

Conclusion
In this study we employed a competitive docking approach for the screening of 103 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 proin ammatory targets (hIL-6, hJAK1) and three SARS-CoV-2 therapeutic targets (s3CL pro , sPL pro and sRdRp). The 8 PDPs ful lled the 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 exibility in a 100 ns simulated dynamics environment. These promising hGRag agonists with anti-in ammatory and SARS-CoV-2 replication inhibitory potential is recommended for further in vitro and in vivo experiments. The authors con rm that the data supporting the ndings of this study are available within the article [and/or] its supplementary materials.