Mitigate the cytokine storm due to the severe COVID-19: A computational investigation of possible allosteric inhibitory actions on IL-6R and IL-1R using selected phytochemicals European Journal of Chemistry 11 (4) (2020) 351-363 European Journal of Chemistry ISSN 2153-2249 (Print) / ISSN 2153-2257 (Online) – Copyright © 2020 The Authors – Atlanta Publishing House LLC – Printed in the USA. This work is published and licensed by Atlanta Publishing House LLC – CC BY NC – Some Rights Reserved. http://dx.doi.org/10.5155/eurjchem.11.4.351-363.2043 European Journal of Chemistry View Journal Online View Article Online Mitigate the cytokine storm due to the severe COVID-19: A computational investigation of possible allosteric inhibitory actions on IL-6R and IL-1R using selected phytochemicals Harindu Rajapaksha 1,*, Bingun Tharusha Perera 1, Jeewani Meepage 1, Ruwan Tharanga Perera 2 and Chithramala Dissanayake 3 1 Department of Chemistry, Faculty of Science, University of Kelaniya, Dalugama, 11 300, Sri Lanka rajapaks_bs14030@stu.kln.ac.lk (H.R.), binguntharushaperera@gmail.com (B.T.P.), jeew321@gmail.com (J.M.) 2 Graduate Studies Division, Gampaha Wickramarachchi Ayurveda Institute, University of Kelaniya, Yakkala, 11870, Sri Lanka 2017_perera@kln.ac.lk (R.T.P.) 3 Department of Cikitsa, Gampaha Wickramarachchi Ayurveda Institute, University of Kelaniya, Yakkala, 11870, Sri Lanka drchithramala@gmail.com (C.D.) * Corresponding author at: Department of Chemistry, Faculty of Science, University of Kelaniya, Dalugama, 11 300, Sri Lanka. e-mail: rajapaks_bs14030@stu.kln.ac.lk (H. Rajapaksha). 10.5155/eurjchem.11.4.351-363.2043 Received: 20 September 2020 Received in revised form: 09 November 2020 Accepted: 12 November 2020 Published online: 31 December 2020 Printed: 31 December 2020 The novel corona virus 2019 (COVID 19) is growing at an increasing rate with high mortality. Meanwhile, the cytokine storm is the most dangerous and potentially life-threatening event related to COVID 19. Phyto-compounds found in existing Ayurveda drugs have the ability to inhibit the Interleukin 6 (IL-6R) and Interleukin 1 (IL-1R) receptors. IL-6R and IL-1R receptors involve in cytokine storm and recognition of phytochemicals with proven safety profiles could open a pathway to the development of the most effective drugs against cytokine storm. In this study, we intend to perform an in silico investigation of effective phyto compounds, which can be isolated from selected medicinal herbs to avoid cytokine storm, inhibiting the IL-6 and IL-1 receptor binding process. An extensive literature survey followed by virtual screening was carried out to identify phytochemicals with potential anti- hyper-inflammatory action. Flexible docking was conducted for validated models of IL-1R and IL-6R-α with the most promising phytochemicals at possible allosteric sites using AutoDock Vina. Molecular dynamics (MD) studies were conducted for selected protein- ligand complexes using LARMD server and conformational changes were evaluated. According to the results, taepeenin J had Gibbs energy (ΔG) of -10.85 kcal/mol towards IL- 1R but had limited oral bioavailability. MD analysis revealed that taepeenin J can cause significant conformational movements in IL-1R. Nortaepeenin B showed a ΔG of -8.5 kcal/mol towards IL-6R-α with an excellent oral bioavailability. MD analysis predicted that it can cause significant conformational movements in IL-6R-α. Hence, the evaluated phytochemicals are potential candidates for further in vitro studies for the development of medicine against cytokine storm on behalf of SARS-COV-2 infected patients. COVID-19 Immunology Cytokine storm Receptor blocking Molecular dynamics Computational chemistry Cite this: Eur. J. Chem. 2020, 11(4), 351-363 Journal website: www.eurjchem.com 1. Introduction Severe Acute Respiratory Syndrome Corona Virus-2 (SARS CoV-2), reported towards the end of December, 2019 has emerged as a serious public health issue [1]. The Emergency Committee of the World Health Organization (WHO) declared an outbreak in China on January 30, 2020, which was considered a Public Health Emergency of International Concern (PHEIC) [2]. By November 09, 2020, 50,728,891 patients have been reported in more than 216 countries around the world with a mortality of 761779 [3]. Although several antiviral drugs are being tested, the lack of a promising drug candidate or a vaccine against SARS-CoV-2 may underscore the urgent need to identify and develop an effective drug. At the onset of the disease, the common symptoms of patients infected with novel Coronavirus 2019 were fever, cough, fatigue, and sputum production [4]. With the development of the disease, severe acute respiratory distress syndrome (ARDS) has been recognized as a major complication similar to diseases caused by SARS-CoV and MERS-CoV [5]. Many studies have shown that SARS-CoV2 can enter cells expressing Angiotensin-converting enzyme 2 (ACE2) receptor protein on their surface [6]. ACE2 molecules on the cell surface are occupied by SARS-CoV2. As a result, Angiotensin 2 (Ang ii) levels increase in the serum because of the inability to bind with ACE2 and reduced ACE2 mediated degradation. SARS-CoV2 itself activates NF-ĸB via pattern recognition receptors (PPRS) [6,7]. Accumulated Ang ii induces inflammatory cytokines including TNFα and IL-6, followed by activation of Interleukin 6 (IL-6) amplifier (IL-6AMP), which leads to an overproduction of inflammatory cytokines initiating a vicious cycle [7]. IL-6 can bind to transmembrane IL-6 receptor (mIL-6R) and soluble IL- ABSTRACT RESEARCH ARTICLE KEYWORDS http://dx.doi.org/10.5155/eurjchem.11.4.351-363.2043 http://www.eurjchem.com/ http://dx.doi.org/10.5155/eurjchem.11.4.351-363.2043 mailto:rajapaks_bs14030@stu.kln.ac.lk mailto:binguntharushaperera@gmail.com mailto:jeew321@gmail.com mailto:2017_perera@kln.ac.lk mailto:drchithramala@gmail.com mailto:rajapaks_bs14030@stu.kln.ac.lk http://www.eurjchem.com/ https://crossmark.crossref.org/dialog/?doi=10.5155/eurjchem.11.4.351-363.2043&domain=pdf&date_stamp=2020-12-31 352 Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 6 receptors (sIL-6R). The resulting complex can combine with signal-transducing component gp130 to activate inflammatory responses via many ways [8], such as activation of the Janus kinase/signal transducer and activator of transcription (JAK/STAT) and MAPK (Mitogen-activated protein kinase) cascades [9]. As a result, inducing a further increase of cytokine release facilitates cytokine storm in SARS-CoV2 patients [10,11]. The IL-6 receptor complex comprises two molecules each of Interleukin 6 (IL-6), Interleukin 6 receptor (IL-6R), and a signaling molecule (gp130) forming a hexametric receptor complex [12]. IL-6 receptors consist of an α subunit and signal- transducing molecule (gp130) known as the β subunit in the receptor complex [12]. However, some cells do not produce the transmembrane IL-6 receptor complex and produce gp130 [8,13,14]. The soluble form of IL-6R is complexed with IL-6, which can stimulate target cells lacking the IL-6 receptor transmembrane complex but expressing gp130 [15]. To stop the cytokine storm of severe COVID-19, IL-6R blockers have been proposed as effective drugs. The drug tocilizumab has been suggested as an effective antagonist to block IL-6 signal transduction pathway [16,17]. Possible effectiveness of Tocilizumab treatment has been revealed by single-cell analysis against severe COVID-19 infection [18] and researchers hypothesized preventing COVID-19 induced pneumonia with possible cytokine blockers [19-21]. Many other studies have hypothesized Interleukin 1 receptor (IL-1R) blockade could control the hyper inflam- mation-associated lung injury in SARS-CoV2. It has been noted the use of monoclonal antibodies against Interleukin 1 (IL-1) to reduce the severity of SARS-CoV2 [22]. And Interleukin 1 (IL-1) blockade by the drug Anakinra showed a significant survival benefit in patients with hyper inflammation [23]. Apart from that, the in silico study has become critical for early stage investigations on drug development, especially antiviral treatments. Binding abilities between the main protease (Mpro) of SARS-COV-2 and most famous antiviral drugs such as chloroquine, hydroxychloroquine, umifenovir, favlplravir and galidesivir have investigated combining with antitussive drug called Noscapine and results depicted hydroxychloroquine-noscpine combination has formed stable binding with Mpro with minimal conformational variations. [24]. Similarly, repurposing usage of Noscapine and natural alkaloids also have been evaluated against SARS-COV-2 in current in silico studies. According to one molecular dynamic studies, noscapine showed remarkable potential and proximal binding to Mpro of SARS-COV-2 than traditional drugs like favlplravir, ribavirin, and chloroquine [25]. Computational investigations play a significant role in antiviral drugs as well as vaccine development processors. As example, the developed multi-epitope vaccine has subjected to in silico studies and docking analysis depicted that it has strong binding to major histocompatibility complex (MHC) receptors (MHC-1 and MHC-2) and virus-specific membrane receptor TLR-2 and molecular dynamic simulations depicted stable binding of vaccine to TLR-2 with low deformation and atomic fluctuation of the complex system expressing the high potency of the multi-epitope vaccine against the global threat of COVID- 19 [26]. Furthermore, molecular dynamic stimulation has used to investigate bindings of multi-epitope vaccines with MHC receptors (MHC-1 and MHC-2) and the virus progression- specific membrane receptor TLR-2 when optimizing a vaccine for human herpes virus-5 [27]. In addition to antiviral investigations, docking studies and molecular dynamic stimulation are utilized to give precautions against bacterial and other pathogenic microorganisms as well. As an example, in silico studies of antimicrobial peptide designing [28]. The intention of this research work is to investigate the effectiveness of Phyto compounds of some selected medicinal herbs and Ayurveda Pharmacopeia to inhibit the IL-6 and IL-1 receptors involved in the cytokine storm. To achieve this objective, a series of theoretical approaches such as homology analysis, molecular docking, and molecular dynamics studies were carried out. 2. Experimental 2.1. Plant review An extensive plant literature review was carried out to identify phytochemicals of medicinal herbs and Ayurveda Pharmacopeia in Sri Lanka [29,30], which may contain anti- hyper-inflammatory effects. The 3-dimensional (3D) structures of the phytochemical compounds of Ayurveda herbs were obtained from PubChem [31]. The acquired 3D structures were submitted to the admetSAR [32] to predict the drug-like properties and toxicities. Twenty-four phytochemical compounds were selected based on adherence to the ADMET predictions and Lipinski's rule of five. 2.2. Selection of crystallographic structures The crystallographic structures of Interleukin-6 receptor α chain (IL-6R-α) (PDB ID: 1N26) and Interleukin-1 receptor (IL- 1R) (PDB ID: 1IRA) were obtained from the RCSB protein database [33] with < 2.50 Å resolution. High energetic loop regions and missing amino acids were corrected by homology modelling followed by loop modelling. 2.3. Homology modelling and loop modelling Basic Local Alignment Search Tool (BLAST) [34] was used to find suitable templates (percentage identity > 95%, E value ≈ 0) for modelling. Modeller 9.24 [35] and Easy-Modeller 4.0 [36] GUIs were used for the model building. Models with the least DOPE score [37] were selected and further refined using the ReFOLD web server [38]. Figure 1; (a) and (b) show multiple sequence alignments of the final models with their correspon- ding templates. The 3D models were validated by VERIFY3D [39], ERRAT [40], PROCHECK [41], and PROVE [42] tests offered by SAVES v5.0 server, (VERIFY3-D - the compatibility of an atomic model (3D) with amino acid sequence, ERRAT - statistics of non-bonded interactions between different atoms, PROVE - the volume of atoms in the models, ProSA [43] - evaluate the quality of protein folds. 2.4. Prediction of binding sites Possible binding sites were predicted by the Computed Atlas of Surface Topography of proteins (CASTp) server [44]. Further refinement of the predicted results was done by blind docking using AutoDock Vina 1.1.2. Two predicted pockets were selected as the suitable binding pocket for further studies. The binding sites of IL-6R-α and IL-1R had grid box absolute coordinates of X = 37.471; Y = 55.161; Z = 60.774 and X = 33.213; Y = 11.591; Z = 42.541, respectively. The binding pockets can be seen in (c) and (d) structures of Figure 1. 2.5. Preparation of molecules Receptor preparations were done using AutoDock Tool 1.5.6. Ligand geometry optimization and energy minimization were performed by Avogadro 2 using MMFF94 [45] force field and conjugate gradient algorithm. 2.6. Parameters of molecular docking Flexible Docking was performed at the above-mentioned binding sites using the selected ligands that passed the virtual screen test and the exhaustiveness parameter was set to 30. Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 353 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 (a) (b) (c) (d) Figure 1. (a) (IL-1R) and (b) (IL-6R-α) show multiple sequence alignments of the final models of with their corresponding templates. (c) and (d) are the 3D structures of IL-1R and IL-6R-α bound to taepeenin J at the chosen binding site. The binding cavities are shown with a 7.00 Å cut off from the ligand. The total number of conformers to be generated was set to 20. Residues in the binding site were selected as flexible residues and the remaining residues were treated as rigid. To ensure the accuracy, the double docking procedure was followed and the root-mean-square deviation (RMSD) was calculated. The docking procedure was repeated with an increased exhaustiveness parameter until the desired RMSD of ≤ 2.5 Å achieved. 2.7. Calculation of electron parameters of selected ligands After the completion of all docking analyzing, the most important ligands were selected and their structures were optimized using GAMESS US version 2019.R1.P1.mkl [46]. Input file preparations were done using Avogadro 2 while DFT calculations were done using B3LYP exchange and correlation hybrid functional, 6-31G basic set, singlet multiplicity, and vacuum phase. After the completion of optimization, the output files were examined using WebMO Demo Server 19.0.009e by the creation of Hückel orbitals to find the energies and shapes of HOMO, HOMO-1, LUMO, and LUMO+1 molecular orbital. 2.8. Molecular dynamic simulation Phytochemicals that showed the best bioavailability and Binding energy ≥ - 8.00 kcal/mol were selected to be examined by molecular dynamics (MD). The highest binding energy conformers were submitted to the LARMD tool (Ligand and Receptor Molecular Dynamics) server [47]. All MD simulations were performed by AMBER16 program [48], where AMBER ff14SB [49] and General AMBER Force Field (GAFF) [50] were used for amino acid residues and ligands, respectively. The protein-ligand complex was solvated in an octahedron box of TIP3P waters [51] and Na+ or Cl- ions were added to the system as counter ions. Sander module in AMBER16 was used for energy minimization [52]. 2000 steps of the steepest descent method and 3000 steps of conjugate gradient method were used in the minimization, followed by the PMEMD module [52] in the MD simulation. Initially, the system was heated from 10 to 300 K at 30 ps. The subsequent release process was identical to the minimization process. Finally, the system was relaxed at 300 K and 1 atm by adding periodic boundary conditions. The simulation time was set to 4 ns, the maximum time allowed by the LARMD server. The RMSD, Rg, and RMSF analysis were performed by CPPTRAJ module [53] in AMBER16. The PCA analysis was done by Bio3d package [54] installed in R [55]. 3. Results and discussion 3.1. Structure-based virtual screening for approved Ayurveda drugs Structure-based virtual scanning was performed for the selected phytochemical compounds of related herbs that are approved to be used in the treatment for infective diseases (Krimi roga) in Ayurveda medicine. Among them, 24 phytochemicals were selected giving priority to the Lipinski's rule of five. According to the Lipinski's rule of five, compounds that showed molecular weight > 500, AlogP > 5, more than 5 hydrogen bond donors, more than 10 acceptor groups, and more than 10 rotational bonds are considered as less efficient for oral administration. Therefore, those compounds have not been considered for further investigations. Among the ADMET predictions, priority was given to human intestinal absorption and human oral bioavailability. According to the literature, most of these phytochemicals are known to be orally administrated. In addition, a special emphasis was given to the carcinogenicity (trinary) and acute oral toxicity when selecting the phytochemical compounds for further analysis (Table 1). 3.2. Model validation VERIFY3-D analyzed the compatibility of an atomic model (3D) with its deleted amino acid sequence. 91.97% and 80.06% of residues of IL-6R and IL-1R had 3D-1D score ≥0.2. The ERRAT server was used for analyzing the statistics of nonbonded interactions between different atom types by comparing with statistics from a highly refined structure. The scores greater than 50 are normally acceptable. Models of IL- 6R-α and IL-1R had a score of 78.69 and 89.70, respectively. In the dihedral angles (Ψ and Φ) of amino acids can have values between -180 and +180°, but the same values are prohibited due to steric interactions between the backbone and side chains. To evaluate this, Ramachandran plots were examined by PROCHECK. 354 Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 Table 1. Virtual screening results of the most promising phytochemicals showing their ADMET prediction and Lipinski rule of five results. No Compound Molecular Formula and structure Lipinski rule of five ADMET Predicted Profile Properties Value Property Probability 1 Conferone * C24H28O4 O O O O Molecular weight (<500 Da) 380.48 Human intestinal absorption 0.989 LogP (<5) 5.15 Human oral bioavailability 0.5143 H bond Acceptor (<10) 4 Carcinogenicity (trinary) 0.6604 H bond donor (<5) 0 Acute oral toxicity (c) 05993(III) Rotational bond (<10) 3 Violations 1 2 Oleanolic acid * C30H48O3 OH HO O Molecular weight (<500 Da) 456.71 Human intestinal absorption 0.9853 LogP (<5) 7.23 Human oral bioavailability 0.5429 H bond Acceptor (<10) 2 Carcinogenicity (trinary) 0.5962 H bond donor (<5) 2 Acute oral toxicity (c) 0.8316(III) Rotational bond (<10) 1 Violations 1 3 Ursolic acid * C30H48O3 HO OH O Molecular weight (<500 Da) 456.71 Human intestinal absorption 0.9853 LogP (<5) 7.09 Human oral bioavailability 0.5143 H bond Acceptor (<10) 2 Carcinogenicity(trinary) 0.5962(III) H bond donor (<5) 2 Acute oral toxicity (c) 0.8316 (III) Rotational bond (<10) 1 Violations 1 4 β-Amyrin * C30H50O O O Molecular weight (<500 Da) 468.77 Human intestinal absorption 0.9906 LogP (<5) 8.74 Human oral bioavailability 0.5571 H bond Acceptor (<10) 2 Carcinogenicity (trinary) 0.471 H bond donor (<5) 0 Acute oral toxicity (c) 0.8704(III) Rotational bond (<10) 1 Violations 1 5 Nortaepeenin B * C20H26O5 OO O O HO Molecular weight (<500 Da) 377.25 Human intestinal absorption 0.9904 LogP (<5) 5.47 Human oral bioavailability 0.5571 H bond Acceptor (<10) 1 Carcinogenicity (trinary) 0.6378 H bond donor (<5) 3 Acute oral toxicity (c) 0.4428(III) Rotational bond (<10) 2 Violations 1 * Ligands that showed significant binding energy. The IL-1R model had only 1.1% of residues in the disallowed region, while the IL-1R model had 1.2% of its residues in the disallowed region. PROVE test statistically analyses the volume of atoms in the models by comparing them with high-quality models. When the percentage of outliers >5%, it is problematic. IL-1R and IL-6R-α models only had an outlier <5%. Additionally, ProSA test could evaluate the quality of protein folds and it was performed on IL-1R and IL-6R-α. They scored a z-score of -6.77 and -6.30, which indicates the proper folding and low energy of the loop regions (Figure 2). 3.3. Molecular docking studies With the aid of the molecular mechanics, structures of the ligands were optimized and energies of the ligands were minimized. Here, the molecular mechanics method was used instead of DFT methods. Because DFT method cannot accurately simulate the molecular surroundings present in the native environment [56], although it costs a higher computational time. It is better to conduct all simulations with an exhaustiveness of 400, but with the limited time available, the exhaustiveness was set to 30 as it had only a quantitative impact and no any qualitative impact on the results obtained [57]. Due to the high degree of rotational freedom of the ligands, it was found that some RMSD values exceeded the recommended 2.00 Å. However, RMSD values reported were below 2.50 Å in all cases. The large protein-protein interaction interface between cytokine receptor and cytokine poses has been identified and it is a real challenge in identifying binding sites for small molecule inhibitor development. Allosteric inhibitor sites are distinct sites on the surface of the enzyme or receptor that are independent of the substrate-binding domain and are a promising target for receptor inhibition [58]. Additionally, compounds targeting allosteric sites give major selectivity advantages over the corresponding orthosteric ligands, including subtype selectivity within receptor families, and may also grant enhanced physicochemical properties [59]. In silico evidence suggests the presence of small molecule allosteric modulator sites on IL-1R [60] and such allosteric sites were discovered in IL-4Rα [59]. Similarly, IL-6R-α may also possess such allosteric sites and they are yet to be discovered. The predicted binding sites were hypothesized as allosteric inhibition sites and ligand binding to such sites may prevent activation of the targeted receptors by causing conformational changes [61]. In order to overcome the activation energy of ligand-induced conformational changes, it would require relatively high binding energies [62], therefore the binding energy cut-off was set to -7.0 kcal/mol and only ligand-receptor interactions which had binding energy < -8.00 kcal/mol were examined by molecular dynamics. Out of all tested phytochemicals, taepeenin J showed the highest average binding affinity (ΔG = -10.85 kcal/mol) toward IL-6R-α and IL-1R during blind docking. Even though it showed limited oral bioavailability, further investigation was done due to its unusually high binding energy. Double docking results showed binding energy of -10.70 kcal/mol (RMSD =1.010 Å) and 11.00 (RMSD =1.211 Å) toward IL-6R-α and IL-1R, respectively. Taepeenin J may be a promising candidate for receptor inhibition if its low bioavailability can be improved. Taepeenin J is 1 of 12 cassane-type diterpenes readily found in seeds of Caesalpinia bonduc [63]. Caesalpinia bonduc seed has been reported to possess anti-inflammatory properties [64], and is currently used as Ayurveda medicine to treat hyper- inflammation [65], β-amyrin, nortapeenin B, teapeenin K, and teapeenin L were the only phytochemicals from Caesalpinia bonduc that obeyed the Lipinski's rule of 5 and had binding energy less than the cut off value. They also showed high oral bioavailability and low toxicity when subjected to the ADMET test. Among them, nortapeenin B had the highest affinity toward IL-6R-α (ΔG = -8.50 kcal/mol, RMSD = 1.735 Å) and β- amyrin had the highest affinity toward IL-1R (ΔG = -8.50 kcal/mol, RMSD = 0.522 Å). Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 355 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 (a) (b) (c) (d) (e) (f) (g) (h) Figure 2. The ERRAT plots of (a) IL-6R-α and (b) IL-1R models indicate their accuracy. Red, yellow, and green colors represent errors, warning, and acceptable scores. (c) and (d) VERIFY3-D plots of (a) IL-6R-α and (b) IL-1R models show that 91.97% and 80.06% of residues of IL-6R and IL-1R have a 3D-1D score ≥0.2. (e) Ramachandran plots of IL-6R-α show that only 1.1% of its residues in the disallowed region. (f) ProSA Z-score diagrams show that the model is within the range of scores typically found for native proteins of similar size. Similarly, (g) and (h) indicate the Ramachandran plot and ProSA Z-score diagram of IL-1R. Teapeenin K and teapeenin L showed ΔG < -8.00 kcal/mol. Table 2 shows the 2D dimensional representation of the binding residues of the receptor responsible for ligand binding. Ferula foetida is a plant that is extensively used in Ayurveda [66] and oligo gum resins of Ferula foetida is highly sorted for its anti-inflammatory effect [67]. Out of all the phytochemicals of Ferula foetida, only Conferone obeyed Lipinski’s rule of 5 and showed high oral bioavailability. Conferone showed ΔG of -9.3 kcal/mol (RMSD = 2.125 Å) toward IL-6R-α. Samarcandinis, which is another phytochemical of Ferula foetida [10] that was tested, and it showed a ΔG of -8.6 kcal/mol toward IL-6R but further investigations were not done due to its limited bioavailability. Both conferone and samarcandinis showed ΔG < -8.00 kcal/mol toward IL-1R. Glycyrrhize glabra stem is widely used in Ayurveda and known for its anti-inflammatory effect [68]. Out of the tested 356 Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 phytochemicals of Glycyrrhize glabra, only 3'-hydroxy-4'-O- methylglabridin [69] obeyed the Lipinski's rule of 5 and showed high oral bioavailability. This compound showed ΔG < -8.00 kcal/mol toward IL-1R and IL-6R. Apart from that, ursolic acid is a natural penta-cyclic triterpenoid carboxylic acid and is the major component of some traditional medicine herbs [70]. Virtual screen results show that it has high bioavailability and drug-likeness. Ursolic acid has shown ΔG < -8.00 kcal/mol toward IL-6R. Oleanolic acid is a natural product that has been isolated from many food and medicinal plants [71]. It has shown significant oral bioavailability and binding ability towards IL-1R and IL-6R-α with ΔG < -8.00 kcal/mol. It is interesting to note that all compounds that showed significant binding are steroids. When considering protein-ligand interactions, it is clear that almost all the ligands have formed hydrophobic interactions with the binding pocket and those interactions are the dominant interactions. Ursolic acid, oleanolic acid and β- amyrin form H-bonds with the binding site. Another way to determine the affinity is to use the inhibition constant Ki to describe the docking values. The Ki value decreases exponentially when the binding energy increases and are qualitatively similar to the IC50 value. Naturally, IC50 values have to be experimental, because they only depend on the concentrations of protein-ligand complexes [72]. A compound or a drug with lower Ki is better, as a small amount is sufficient to maintain a stable conformation at the receptor site, therefore, a better inhibition rate. Once again, this shows that teapeenin J is the best inhibitor as it has a Ki of 0.0160 µM (toward IL-6R-α) and 0.0097 µM (toward IL-1R), which is approximately 9 times better than conferone, and 56 times better than ursolic acid, respectively. However, because of the low bioavailability of teapeenin J, conferone (Ki = 0.168 µM) may be the better choice for the inhibition of IL-6R-α and ursolic acid for IL-1R (Ki = 0.543 µM). Similarly, the binding process of the selected ligands to Glycoprotein 130 (IL-6-β) (a transmembrane protein), which is the founding member of the class of all cytokine receptors [73] was also examined by using the same procedure mentioned above. It was found that nortaepeenin B, oleanolic acid, ursolic acid, and β-amyrin were binding to IL-6-β with ΔG of -8.90 kcal/mol, -8.00 kcal/mol, -8.10 kcal/mol and -8.10 kcal/mol, respectively. Conferone had a ΔG < -7.00 kcal/mol and may be a good candidate if IL-6R-α needs to be selectively inhibited. The results show that the binding pockets have limited solvent accessibility and are hydrophobic in nature (Table 2). Moreover, it can be clearly seen that most of the ligands with significant binding energies are hydrophobic in nature and lack the ability to make a wide range of non-covalent interactions like ionic interactions. As results, the hydrophobic motifs play the most significant role in ligand binding. Furthermore, the commercial hyperinflammatory drugs prednisone, dexamethasone, and methylprednisolone were docked with the predicted binding sites of IL-1R, IL-6R and IL- 6-β. The results revealed that prednisone can bind with the identified binding sites of IL-1R, IL-6R, and IL-6-β with ΔG of - 7.5, -7.5 and -6.7 kcal/mol, respectively. dexamethasone was found to be interacting with IL-1R, IL-6R, and IL-6-β with ΔG of -7.7, -7.4 and -6.1 kcal/mol accordingly. Finally, methylpred- nisolone interacted with IL-6R and IL-6-β with ΔG < -7.00 kcal/mol. However, methylprednisolone showed ΔG of -7.8 kcal/mol with IL-1R. These ΔG values depicted that the investigated phytochemicals in this study have higher potential to act as allosteric inhibitors of IL-1R, IL-6R than the analyzed commercial anti-inflammatory drugs. Performed docking simulations revealed that above-cited compounds showed significant binding affinity to possible allosteric sites of IL-6R-α and IL-1R. However, the docking score is purely dependent on thermodynamic factors. Therefore, molecular dynamics experiments and in vitro studies are required to arrive at a definite conclusion. 3.4. Analysis of electronic properties of selected ligands The Molecular Electrostatic Potential Surface (MEPS) indicates the relative electrostatic potential of the constant electron density surface and is used to predict sites of the molecule that are prone to electrophilic attack and nucleophilic attacks [74, 75]. MAPS plots express the charge distribution profile of molecules, which helps to identify the possible sites of ligand-receptor interactions. Figure 3 shows MEPS maps of the 6 selected compounds. The negative electrostatic potentials (red color) represent sites that have loosely bound or excess electrons and that are likely to react with electrophiles (protons). Positive electrostatic potentials represent a (blue color) site that is electron-deficient and that are likely to react with nucleophiles depicted in blue color. The relative order of nucleophilicity of functional groups in amino acids is R-S− > R- NH2 > R-COO− = R-O− [76]. The electrostatic potential increases in the order of red < orange < yellow < green < blue. MEPS of taepeenin J indicates a positive potential over the whole molecule except in areas where carbonyl oxygen is present. This observation is somewhat unique as other molecules did not show such positive potential over the whole molecule. This may be a reason for the high binding energy of it, as it can bind with many nucleophilic functional groups. Conferone was somewhat similar to taepeenin J as it has increased positive potential over the ring area. β-Amyrin only had a single positive potential patch around the carbonyl carbon. In the case of ursolic acid and oleanolic acid, positive potentials were only observed in acidic and alcoholic protons, indicating their ability to undergo deprotonation. Nortapeenin B showed negative potential around oxygen which forms a double bond. In characterizing molecular behavior, the energy difference between the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO) is a very important descriptor. FMO (Frontier Molecular Orbitals) attempts to forecast relative reactivity based on the reactant's properties [77]. HOMO orbitals appear to concentrate mainly close to the aromatic group, a functional group capable of performing dipole-dipole, aromatic-hydrophobic, and ion- dipole interactions, with a high likelihood of interacting with amino acids present in the receptor periphery. Energy differences between HOMO and LUMO increase in the order of taepeenin J (2.20 eV) < conferone (2.30 eV) < β-amyrin (3.15 eV) < oleanolic acid (3.18 eV) < nortapeenin B (3.19 eV) < ursolic acid (3.21 eV) (Figure 4). Lower the HOMO-LUMO gap, greater the reactivity [78]. Therefore, taepeenin J must have the highest reactivity and this may be the reason for its highest binding energy compared to other tested compounds. As indicated by the result, conferone has relatively high reactivity and agrees with the docking results. According to the results, oleanolic acid, nortapeenin B, and ursolic acid showed similar reactivities, and may be that is the reason for them having similar binding energies. All tested compounds had negative LUMO energies and it indicates that they have high electron acceptability and readily undergo nucleophilic addition (react with HOMO of nucleophile amino acid residues). When there is more than one potential attack centre, the chosen mode of reaction would be with the centre having the largest LUMO orbital [56]. Similarly, in the case of electrophilic attack, the reaction would be with the largest HOMO of the compound and the largest LUMO of electrophiles. 3.5. Molecular dynamics studies Ligands that cause conformational changes in receptors usually show high fluctuations in Root Mean Square Deviation Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 357 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 Table 2. Molecular docking results of the most promising phytochemicals indicating their binding pocket, ΔE (kcal/mol), Ki (µmol), and RMSD (Å). Compounds IL-6 Receptor IL-1 Receptor Conferone ΔE = -9.3 Ki = 0.168 RMSD = 2.1245 ΔE = -7.0 Ki = 4.81 RMSD = 2.36 Oleanolic acid ΔE = -8.3 Ki = 0.9 RMSD = 0.648 ΔE = -8.4 Ki = 0.76 RMSD = 0.760 Ursolic acid ΔE = -7.8 Ki = 2.08 RMSD = 1.562 ΔE = -8.6 Ki = 0.54 RMSD = 1.723 Β-Amyrin ΔE = -7.7 Ki = 2.46 RMSD = 0.357 ΔE = -5.8 Ki = 0.64 RMSD = 0.522 Nortaepeenin B ΔE = -8.5 Ki = 0.64 RMSD = 1.735 ΔE = -8.0 Ki = 1.49 RMSD = 1.578 Taepeenin J ΔE = -10.7 Ki = 0.016 RMSD = 0.028 ΔE = -11.0 Ki = 0.0097 RMSD = 0.015 358 Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 (a) (b) (c) (d) (e) (f) Figure 3. Electrostatic potential maps (MEPS) of (a) β-amyrin, (b) conferone, (c) nortaepeenin B, (d) oleanolic acid, (e) ursolic acid and, (f) teapeenin J. The electrostatic potential increases in the order of red < orange < yellow < green < blue. Images are not scaled proportionally. (a) (b) (c) (d) (e) (f) Figure 4. Visualization of frontier molecular orbitals (FMO) a) β-amyrin, (b) conferone, (c) nortaepeenin B, (d) oleanolic acid, (e) ursolic acid and, (f) teapeenin J. The energy of each level is shown in the left-hand corner. Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 359 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 (a) (b) (c) (d) Figure 5. Molecular dynamics results of IL-6R-α with (a) conferone, (b) nortaepeenin B, (c) oleanolic acid, and, (d) teapeenin J. Left, middle and right columns indicate RMSD, RMSF, and Rg, respectively. (RMSD), Root Mean Square Fluctuation (RMSF), and Radius of Gyration (Rg) graphs [79]. When calculating the Free energy contribution of electrostatic energy (ELE), Van der Waals (VDW), total gas-phase energy, nonpolar and polar contributions to solvation (PBSOL and GBSOL) and entropic (T∆S) states were taken into account. Even though the complex has not attained equilibrium after 4 ns of simulation time, RMSD plots of receptor-ligand complexes show high flexibility (conformational mobility) of the backbone of the receptor. When considering the RMSD data of IL-6R-α protein with different ligands (Figure 5), they differ in magnitude, indicating variations in fluctuations with different ligands. Therefore, it can be concluded that the ligands may play a critical role in the conformational motion. A similar observation can be seen in the magnitude of RMSD data of IL- 1R proteins with different ligands. Comparison of the RMSD graphs for IL-6R-α protein with different ligands show that with ligands conferone and oleanolic acid receptor backbone shows comparatively less RMSD with average value of 2.1089 Å in the presence of ligand conferone and average value of 2.0123 Å in the presence of ligand oleanolic acid. Therefore, it can be stated that IL-6R-α protein is comparatively more stable with ligands conferone and oleanolic acid, than with ligands nortaepeenin B and teapeenin J. According to the RMSD graphs of IL-1R protein with different ligands (Figure 6), it can be seen that the receptor is not well equilibrated during the simulation time (4 ns). Therefore, it is difficult to draw a conclusion depending on RMSDs of IL-1R. However, the average RMSD values of the receptor with ligands ursolic acid and oleanolic acid show comparatively less values compared to the ligands β-amyrin and teapeenin J, revealing a somewhat similar trend which can be seen in IL-6R-α protein. These results according RMSD to reveal the effects of ligands on protein structure. To get further details, RMSF calculations were carried out. RMSF data provides information about the important residues for the fluctuation in conformation. According to the results obtained, it can be seen that the residues in a site away from the ligand binding site are highly flexible, indicating that those are residues responsible for the conformational motion. It also shows a high mobility in loop regions as expected. In order to understand the inhibition mechanism of these compounds, a longer simulation time must be used. The binding free energy calculated by MM/PB(GB)SA agrees with the binding energy trend predicted by docking. In IL-6R-α, Teapeenin J had the highest ∆GPB and ∆GGB with respect to IL-1R and IL-6R-α (For teapeenin J with IL-6R-α; ∆GPB = - 33.44 kcal/mol and ∆GGB = -38.34 kcal/mol and for teapeenin J with IL-1R; ∆GPB = -20.8 kcal/mol and ∆GGB = -28.13 kcal/mol). The binding free energy trend of β-amyrin and oleanolic acid to IL-1R was reversed, but this minor deviation can be neglected as they have a negligible energy difference. To gain insight into the hydrogen bond (H-bond) formation between ligand and protein, H-bond analysis was performed. Teapeenin J was H-bonding with LYS 45, SER 72, SER 122, LEU 123, and THR 124 residues of IL-6R-α in an alternating manner. However, in the case of IL-1R, teapeenin J only showed H- bonding with LYS 60 and in both these cases teapeenin J acted as an acceptor where the average H-bond length was 3.20 Å. On the other hand, oleanolic acid forms H-bonds with GLY 14, LYS 45, GLY 73, ASN 74, HIS 88, LEU 123, and THR 124 of IL-6R-α. 360 Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 (a) (b) (c) (d) Figure 6. Molecular dynamics results of IL-1R with (a) ursolic acid, (b) β-amyrin, (c) oleanolic acid and (d) teapeenin J. Left, middle and right columns indicate RMSD, RMSF, and Rg, respectively. (a) (b) Figure 7. PCA results trajectories for complexes of IL-6R-α with (a) conferone and (b) nortaepeenin B. The color changes from blue to white to red as the complex evolves with the time. Similarly, it had formed H-bonds with IL-1R via ARG 22, PRO 23, SER 90, and TYR 124 alternating manner. The oxygen atoms in the blue color region shown in Figures 7-9 acts as donor, while the oxygen atom in the red color region acted as an acceptor and the average H-bond length was 3.18 Å. Similar to oleanolic acid, nortaepeenin B also showed both accepter and donor abilities. Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 361 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 (a) (b) Figure 8. PCA results trajectories for complexes of IL-6R-α with (a) oleanolic acid and (b) teapeenin J. The color changes from blue to white to red as the complex evolves with the time. (a) (b) (c) (d) Figure 9. PCA results trajectories for complexes of IL-1R with (a) ursolic acid, (b) β-amyrin, (c) oleanolic acid, and (d) teapeenin J. The color changes from blue to white to red as the complex evolves with the time. 362 Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 It has formed H-bonds with IL-6R-α via GLY 14, LYS 45, GLY 73, ASN 74, HIS 88, LEU 123, and THR 124 with an average bond length of 3.19 Å. In contrast, conferone only had acceptor ability and it formed H-bond (Average bond length = 3.22 Å) with LYS 45, ALA 47, SER 72, GLY 73 and LEU 89 of IL-6R-α an alternating manner. β-Amyrin formed H-bonds with ARG 22 and LEU 12 of IL-1R (Average bond length = 3.30 Å). Finally, ursolic acid formed H-bonds with ARG 22, PRO 23, and SER 90 of IL-1R (Average bond length = 3.22 Å). Normally, the protein-ligand H- bond distance ranges between 2.70 to 3.10 Å [80,81]. However, all tested ligands formed longer H-bonds. Therefore, it can be concluded that none of the above ligands were strong H bonders. This observation confirmed the docking results where no significant were not detected. Furthermore, this showed that hydrophobic interactions were the most significant type of interaction. As the computational power is limited, principal correlation analysis (PCA) was performed to further investigate the conformational motion of the complexes even though they were not at equilibrium. Evidently, the first 3 principal components (PCs) accounted for >60% variance except in the complexes of IL-6R-α with conferone and oleanolic acid. The two states of confirmation are represented by blue and red colors. The color changes from blue to white to red as the complex evolves with the time. According to PCA results, only IL-6R-α-nortaepeenin B and IL-1R-teapeenin J complexes showed 2 distinct clusters along PC1 axis (Figures 7-9). Other complexes showed thermal motion caused due to not being properly equilibrated and no distinct separations of the clusters were observed. In summary, MD results confirmed the binding free energy trend predicted by docking, and RMSD results together with PCA results indicated that nortaepeenin B and teapeenin J can cause considerable conformational mobility in IL-6R- α and IL- 1R, respectively. 4. Conclusion Owing to the augmented SARS CoV-2 viral spread among the continents in the world, it is worthwhile to discover ways to protect affected patients from the complications. To proceed with the task, investigation of effective antagonists for the IL-6 signal transduction pathway is important to prevent covid-19 induced pneumonia and other related complications. This docking survey will address the emphasized issue and based on orally administrated capability of phytochemicals such as conferone, ursolic acid, β-amyrin, and nortaepeenin B and taepeenin J, which extracted from commonly used medicinal herbs in Ayurveda medicine. The above-mentioned phyto chemicals were used in the study as their binding energies showed high partiality towards IL-1R and IL-6R receptors. In addition to that, MD studies depicted that nortaepeenin B and taepeenin J can cause considerable conformational mobility in IL-6R and IL-1R, respectively, verifying the possibility of using them as allosteric inhibitors of the concerned receptors and ΔG values depict that some investigated phytochemicals in this study have higher potential to act as allosteric inhibitors of IL-1R, IL-6R than analyzed commercial anti-inflammatory drugs such as Prednisone, Dexamethasone and Methylprednisolone. Apart from that, a follow-up study with higher simulation time may provide detailed insights of the mechanism of action of these phytochemicals. Due to the time limitation, the study was limited to some selected herbs in Sri Lanka. A follow-up study with different medicinal herbs from different countries around the world is recommended. Based on the observations made in this study, it can be concluded that there is a reasonable probability of a promising drug for hyper-inflammation using medicinal herbs to combat SARS CoV-2 viral pandemic. Acknowledgements Author thanks to Prof. Channa De Silva, Mrs. Sachini Rathnasekara, Miss. Hiruni Jayasekera and Mr. Sudesh Hemal for providing language help, writing assistance, and proof reading the article. Disclosure statement Conflict of interest: The authors declare that they have no conflict of interest. Author contributions: All authors contributed equally to this work. Ethical approval: All ethical guidelines have been adhered. Sample availability: Samples of the compounds are available from the author. ORCID Harindu Rajapaksha http://orcid.org/0000-0002-3988-8882 Bingun Tharusha Perera http://orcid.org/0000-0002-5671-8279 Jeewani Meepage http://orcid.org/0000-0001-7723-0443 Ruwan Tharanga Perera http://orcid.org/0000-0001-5695-6579 Chithramala Dissanayake http://orcid.org/0000-0002-2653-4049 References [1]. Wu, F.; Zhao, S.; Yu, B.; Chen, Y. M.; Wang, W.; Song, Z. G.; Hu, Y.; Tao, Z. W.; Tian, J. H.; Pei, Y. Y.; Yuan, M. L.; Zhang, Y. L.; Dai, F. H.; Liu, Y.; Wang, Q. M.; Zheng, J. J.; Xu, L.; Holmes, E. C.; Zhang, Y. Z. Nature 2020, 579(7798), 265-269. [2]. Chan, J. F. W.; Kok, K. H.; Zhu, Z.; Chu, H.; To, K. K. W.; Yuan, S.; Yuen, K. Y. Emerg. Microbes. Infect. 2020, 9(1), 221-236. [3]. Culp, W. C., Jr. A & A Practice 2020, 14(6), e01218. [4]. Huang, C.; Wang, Y.; Li, X.; Ren, L.; Zhao, J.; Hu, Y.; Zhang, L.; Fan, G.; Xu, J.; Gu, X.; Cheng, Z.; Yu, T.; Xia, J.; Wei, Y.; Wu, W.; Xie, X.; Yin, W.; Li, H.; Liu, M.; Xiao, Y.; Gao, H.; Guo, L.; Xie, J.; Wang, G.; Jiang, R.; Gao, Z.; Jin, Q.; Wang, J.; Cao, B. Lancet 2020, 395(10223), 497-506. [5]. de Wit, E.; van Doremalen, N.; Falzarano, D.; Munster, V. J. Nat. Rev. Microbiol. 2016, 14(8), 523-534. [6]. Zhou, P.; Yang, X. L.; Wang, X. G.; Hu, B.; Zhang, L.; Zhang, W.; Si, H. R.; Zhu, Y.; Li, B.; Huang, C. L.; Chen, H. D.; Chen, J.; Luo, Y.; Guo, H.; Jiang, R. D.; Liu, M. Q.; Chen, Y.; Shen, X. R.; Wang, X.; Zheng, X. S.; Zhao, K.; Chen, Q. J.; Deng, F.; Liu, L. L.; Yan, B.; Zhan, F. X.; Wang, Y. Y.; Xiao, G. F.; Shi, Z. L. Nature 2020, 579(7798), 270-273. [7]. Hirano, T.; Murakami, M. Immunity 2020, 52(5), 731-733. [8]. Tanaka, T.; Narazaki, M.; Kishimoto, T. Immunotherapy 2016, 8(8), 959-970. [9]. Heinrich, P. C.; Behrmann, I.; Haan, S.; Hermanns, H. M.; Müller- Newen, G.; Schaper, F. Biochem. J. 2003, 374(1), 1-20. [10]. Rajapaksa, R. M. H.; Perera, B. T.; Nisansala, M. J.; Perera, W. P. R. T.; Dissanayake, K. G. C. Global J. Eng. Sci. Res. Manag. 2020, 7, 51-61. [11]. Radbel, J.; Narayanan, N.; Bhatt, P. J. Chest 2020, 158(1), e15-e19. [12]. Varghese, J. N.; Moritz, R. L.; Lou, M.-Z.; van Donkelaar, A.; Ji, H.; Ivancic, N.; Branson, K. M.; Hall, N. E.; Simpson, R. J. Proceed. National Acad. Sci. 2002, 99(25), 15959-15964. [13]. Ward, L. D.; Howlett, G. J.; Discolo, G.; Yasukawa, K.; Hammacher, A.; Moritz, R. L.; Simpson, R. J. J. Biol. Chem. 1994, 269, 23286-23289. [14]. Scheller, J.; Chalaris, A.; Schmidt-Arras, D.; Rose-John, S. Biochim. Biophys. Acta (BBA) – Mol. Cell Res. 2011, 1813(5), 878-888. [15]. McGonagle, D.; Sharif, K.; O’Regan, A.; Bridgewood, C. Autoimmunity Rev. 2020, 19(6), 102537. [16]. Zhang, C.; Wu, Z.; Li, J. W.; Zhao, H.; Wang, G. Q. Int. J. Antimicrob. Agent. 2020, 55(5), 105954. [17]. Lima de Oliveira, M. D.; Teixeira de Oliveira, K. M. ChemRxiv 2020, Retrieved Oct 10, 2020, https://doi.org/10.26434/chemrxiv.12044538.v4 [18]. Guo, C.; Li, B.; Ma, H.; Wang, X.; Cai, P.; Yu, Q.; Zhu, L.; Jin, L.; Jiang, C.; Fang, J.; Liu, Q.; Zong, D.; Zhang, W.; Lu, Y.; Li, K.; Gao, X.; Fu, B.; Liu, L.; Ma, X.; Weng, J.; Wei, H.; Jin, T.; Lin, J.; Qu, K. Nat. Commun. 2020, 11(1), https://doi.org/10.1038/s41467-020-17834-w http://orcid.org/0000-0002-3988-8882 http://orcid.org/0000-0002-5671-8279 http://orcid.org/0000-0001-7723-0443 http://orcid.org/0000-0001-5695-6579 http://orcid.org/0000-0002-2653-4049 https://doi.org/10.26434/chemrxiv.12044538.v4 https://doi.org/10.1038/s41467-020-17834-w Rajapaksha et al. / European Journal of Chemistry 11 (4) (2020) 351-363 363 2020 – European Journal of Chemistry – CC BY NC – DOI: 10.5155/eurjchem.11.4.351-363.2043 [19]. Wang, J.; Qiao, C.; Xiao, H.; Lin, Z.; Li, Y.; Zhang, J.; Shen, B.; Fu, T.; Feng, J. Drug Des. Devel. Ther. 2016, 10, 4091-4100. [20]. Monteleone, G.; Sarzi-Puttini, P. C.; Ardizzone, S. Lancet Rheumatol. 2020, 2(5), e255-e256. [21]. Liu, B.; Li, M.; Zhou, Z.; Guan, X.; Xiang, Y. J. Autoimmunity 2020, 111, 102452. [22]. AbdelMassih, A. F.; Ramzy, D.; Nathan, L.; Aziz, S.; Ashraf, M.; Youssef, N. H.; Hafez, N.; Saeed, R.; Agha, H. Cardiovasc. Endocrinol. Metabol. 2020, 9(3), 121-124. [23]. Velavan, T. P.; Meyer, C. G. Trop. Med. Int. Health 2020, 25(3), 278- 280. [24]. Kumar, N.; Awasthi, A.; Kumari, A.; Sood, D.; Jain, P.; Singh, T.; Sharma, N.; Grover, A.; Chandra, R. J. Biomol. Struc. Dynam. 2020, 1- 16. [25]. Kumar, N.; Sood, D.; van der Spek, P. J.; Sharma, H. S.; Chandra, R. J. Proteome Res. 2020, 19(11), 4678-4689. [26]. Kumar, N.; Sood, D.; Chandra, R. RSC Adv. 2020, 10(59), 35856- 35872. [27]. Kumar, N.; Sood, D.; Tomar, R.; Chandra, R. ACS Omega 2019, 4(25), 21370-21380. [28]. Kumar, N.; Sood, D.; Chandra, R. ACS Pharmacol. Transl. Sci. 2020. https://doi.org/10.1021/acsptsci.0c00139. [29]. Dissanayake, K.G.C.; Perera, W. P. R. T. Res. J. Med. Plant 2020, 8(2), 135-139. [30]. Dissanayake, L.; Perera, P.; Attanayaka, T.; Heberle, E.; Jayawardhana, M. Plants 2020, 9(10), 1315. [31]. PubMed, U. S. National Library of Medicine, Retrieved Oct 10, 2020, from https://pubchem.ncbi.nlm.nih.gov [32]. Cheng, F.; Li, W.; Zhou, Y.; Shen, J.; Wu, Z.; Liu, G.; Lee, P. W.; Tang, Y. J. Chem. Inf. Model. 2012, 52(11), 3099-3105. [33]. Berman, H. M. Nucleic Acids Res. 2000, 28(1), 235-242. [34]. Basic Local Alignment Search Tool, National Center for Biotechnology Information, U. S. National Library of Medicine, Retrieved Oct 10, 2020, from https://blast.ncbi.nlm.nih.gov/Blast.cgi [35]. Eswar, N.; Webb, B.; Marti-Renom, M. A.; Madhusudhan, M. S.; Eramian, D.; Shen, M.; Pieper, U.; Sali, A. Curr. Protoc. Bioinform. 2006, 15(1), 5.6.1-5.6.30. [36]. Kuntal, B. K.; Aparoy, P.; Reddanna, P. BMC Res Notes 2010, 3(1), 226. [37]. Shen, M.; Sali, A. Protein Sci. 2006, 15(11), 2507-2524. [38]. Shuid, A. N.; Kempster, R.; McGuffin, L. J. Nucleic Acids Res. 2017, 45(W1), W422-W428. [39]. Eisenberg, D.; Lüthy, R.; Bowie, J. U. Meth. Enzymol. 1997, 277, 396- 404. [40]. Colovos, C.; Yeates, T. O. Protein Sci. 1993, 2(9), 1511-1519. [41]. Laskowski, R. A.; MacArthur, M. W.; Moss, D. S.; Thornton, J. M. J. Appl. Cryst. 1993, 26(2), 283-291. [42]. Torrisi, M.; Kaleel, M.; Pollastri, G. BioRxiv 2018, 289033. [43]. Sippl, M. J. Proteins 1993, 17(4), 355–362. [44]. Tian, W.; Chen, C.; Lei, X.; Zhao, J.; Liang, J. Nucleic Acids Res. 2018, 46(W1), W363-W367. [45]. Halgren, T. A. J. Comput. Chem. 1996, 17 (5-6), 490-519. [46]. Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory. In Theory and Applications of Computational Chemistry; Elsevier, 2005; pp 1167-1189. [47]. Yang, J. F.; Wang, F.; Chen, Y. Z.; Hao, G. F.; Yang, G. F. Brief. Bioinformatics 2019, bbz141, https://doi.org/10.1093/bib/bbz141. [48]. Marques, P. R. B. de O.; Yamanaka, H. Quim. Nova 2008, 31(7), 1791- 1799. [49]. Plewczynski, D.; Lazniewski, M.; Augustyniak, R.; Ginalski, K. J. Comput. Chem. 2010, 32(4), 742-755. [50]. Wenthur, C. J.; Gentry, P. R.; Mathews, T. P.; Lindsley, C. W. Annu. Rev. Pharmacol. Toxicol. 2014, 54(1), 165-184. [51]. Wold, E. A.; Chen, J.; Cunningham, K. A.; Zhou, J. J. Med. Chem. 2018, 62(1), 88-127. [52]. Yang, C. Y. PLoS ONE 2015, 10(2), e0118671. [53]. Case, D.; Betz, R.; Cerutti, D. S.; Cheatham, T.; Darden, T.; Duke, R.; Giese, T. J.; Gohlke, H.; Götz, A.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T.-S.; LeGrand, S.; Li, P.; Lin, C.; Luchko, T.; Kollman, P. Amber 2016, University of California, San Francisco, 2016. [54]. Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. J. Chem. Theory Comput. 2015, 11(8), 3696-3713. [55]. Wang, B.; Merz, K. M. J. Chem. Theory Comput. 2005, 2(1), 209-215. [56]. Price, D. J.; Brooks, C. L. J. Chem. Phys. 2004, 121(20), 10096-10103. [57]. Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26(16), 1668-1688. [58]. Roe, D. R.; Cheatham, T. E. J. Chem. Theory Comput. 2013, 9(7), 3084- 3095. [59]. Grant, B. J.; Rodrigues, A. P. C.; ElSawy, K. M.; McCammon, J. A.; Caves, L. S. D. Bioinformatics 2006, 22(21), 2695-2696. [60]. R Core Team, R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, 2018. [61]. Kayne, F. J.; Price, N. C. Biochemistry 1972, 11(23), 4415-4420. [62]. Nayak, T. K.; Vij, R.; Bruhova, I.; Shandilya, J.; Auerbach, A. J. General Physiol. 2019, 151(4), 465-477. [63]. Cheenpracha, S.; Srisuwan, R.; Karalai, C.; Ponglimanont, C.; Chantrapromma, S.; Chantrapromma, K.; Fun, H.-K.; Anjum, S.; Atta- ur-Rahman. Tetrahedron 2005, 61 (36), 8656-8662. [64]. Kannur, D.; Sonavane, L.; Khandelwal, K.; Paranjpe, M.; Dongre, P. J. Adv. Pharm. Tech. Res. 2012, 3(3), 171. [65]. Mehta, P.; McAuley, D. F.; Brown, M.; Sanchez, E.; Tattersall, R. S.; Manson, J. J. Lancet 2020, 395(10229), 1033-1034. [66]. Mahendra, P.; Bisht, S. Phcog. Rev. 2012, 6(12), 141-146. [67]. Bagheri, S. M.; Hedesh, S. T.; Mirjalili, A.; Dashti-R, M. H. J. Evid. Based Complementary Altern. Med. 2016, 21(4), 271-276. [68]. Shah, S. L.; Wahid, F.; Khan, N.; Farooq, U.; Shah, A. J.; Tareen, S.; Ahmad, F.; Khan, T. Alternative Med. 2018, 2018, 1-8. [69]. Dissanayake, K. G. C.; Weerakoon, W. M. T. D. N., Perera, W. P. R. T. Int. J. Sci. Basic Appl. Res. 2020, 51(1), 1-11. [70]. Ikeda, Y.; Murakami, A.; Ohigashi, H. Mol. Nutr. Food Res. 2008, 52(1), 26-42. [71]. Ayeleso, T.; Matumba, M. Molecules 2017, 22(11), 1915. [72]. Hulme, E. C.; Trevethick, M. A. British J. Pharm. 2010, 161(6), 1219- 1237. [73]. Karczewska, A.; Nawrocki, S.; Brborowicz, D.; Filas, V.; Mackiewicz, A. Cancer 2000, 88 (9), 2061-2071. [74]. Murray, J. S.; Politzer, P. WIREs Comput. Mol. Sci. 2011, 1(2), 153-163. [75]. Anand, K.; Khan, F. I.; Singh, T.; Elumalai, P.; Balakumar, C.; Premnath, D.; Lai, D.; Chuturgoon, A. A.; Saravanan, M. ACS Omega 2020, 5(29), 17973-17982. [76]. Bischoff, R.; Schlüter, H. J. Proteomics 2012, 75(8), 2275-2296. [77]. Zhuo, L. G.; Liao, W.; Yu, Z. X. Asian J. Org. Chem. 2012, 1(4), 336-345. [78]. Aihara, J. J. Phys. Chem. A 1999, 103(37), 7487-7495. [79]. Mukund, V.; Behera, S. K.; Alam, A.; Nagaraju, G. P. Bioinformation 2019, 15(1), 11-17. [80]. Stryer, L.; Gumport, R. I. Student companion for Stryer's biochemistry - Biochemistry, New York, N.Y., Freeman, 1995. [81]. Poznanski, J.; Poznanska, A.; Shugar, D. PLoS ONE 2014, 9(6), e99984. Copyright © 2020 by Authors. This work is published and licensed by Atlanta Publishing House LLC, Atlanta, GA, USA. The full terms of this license are available at http://www.eurjchem.com/index.php/eurjchem/pages/view/terms and incorporate the Creative Commons Attribution-Non Commercial (CC BY NC) (International, v4.0) License (http://creativecommons.org/licenses/by-nc/4.0). By accessing the work, you hereby accept the Terms. This is an open access article distributed under the terms and conditions of the CC BY NC License, which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited without any further permission from Atlanta Publishing House LLC (European Journal of Chemistry). No use, distribution or reproduction is permitted which does not comply with these terms. Permissions for commercial use of this work beyond the scope of the License (http://www.eurjchem.com/index.php/eurjchem/pages/view/terms) are administered by Atlanta Publishing House LLC (European Journal of Chemistry). https://doi.org/10.1021/acsptsci.0c00139 https://pubchem.ncbi.nlm.nih.gov/ https://blast.ncbi.nlm.nih.gov/Blast.cgi http://www.eurjchem.com/index.php/eurjchem/pages/view/terms http://creativecommons.org/licenses/by-nc/4.0 http://www.eurjchem.com/index.php/eurjchem/pages/view/terms 1. Introduction 2. Experimental 2.1. Plant review 2.2. Selection of crystallographic structures 2.3. Homology modelling and loop modelling 2.4. Prediction of binding sites 2.5. Preparation of molecules 2.6. Parameters of molecular docking 2.7. Calculation of electron parameters of selected ligands 2.8. Molecular dynamic simulation 3. Results and discussion 3.1. Structure-based virtual screening for approved Ayurveda drugs 3.2. Model validation 3.3. Molecular docking studies 3.4. Analysis of electronic properties of selected ligands 3.5. Molecular dynamics studies 4. Conclusion Acknowledgements Disclosure statement ORCID References PrintField10: PrintField11: PrintField12: PrintField13: PrintField14: PrintField15: PrintField16: PrintField17: PrintField18: PrintField19: PrintField110: PrintField111: PrintField112: PrintField20: PrintField21: PrintField22: PrintField23: PrintField24: PrintField25: PrintField26: PrintField27: PrintField28: PrintField29: PrintField210: PrintField211: PrintField212: