Phytochemical constituents of Inula britannica as potential inhibitors of dihydrofolate reductase: A strategic approach against shigellosis

Abstract Shigella dysenteriae type 1 is considered as an epidemic in different developing countries, which is responsible for the most severe form of bacterial dysentery. It habitually can develop to the most severe form of dysentery with deadly complications. Development of drugs against this disease is still ongoing. Therefore, we used in silico studies to screen the Inula britannica phytocompounds that are used in traditional Chinese and Kampo Medicines and have activities against different diseases. Spinacetin, eupatin, chrysoeriol and diosmetin were successfully passed through the docking-based screening and absorption, distribution, metabolism, excretion and toxicity (ADMET) filtration. The estimated docking affinities of eupatin, diosmetin, chrysoeriol and spinacetin with Dihydrofolate reductase type 1 (DHFR-1), were −6.5, −6.5, −6.3 and −6.1 kcal/mol, respectively. Which were selected for further investigations based on their favorable ADME/Tox characteristics. Then, the 100 ns molecular dynamics (MD) simulations of apo DHFR, spinacetin-DHFR, eupatin-DHFR, chrysoeriol-DHFR and diosmetin-DHFR complexes were carried out. The RMSD fluctuations of the spinacetin, eupatin, chrysoeriol and diosmetin inside the binding site were explored. Subsequently, the effect of binding Spinacetin, eupatin, chrysoeriol and diosmetin upon the dynamic stability of protein was assessed. Additionally, Principal Component Analysis (PCA) and Hydrogen bond analysis was performed for the apo protein and the protein ligand complexes. The results revealed that chrysoeriol and eupatin has good inhibitory effects against DHFR-1 as treatment for Shigella dysenteriae type when compared to other compounds under study. Hence this study implies that eupatin and chrysoeriol are a significantly potential drug like molecule for the treatment of Shigellosis and must undergo validation through in vivo and in vitro experiments. Communicated by Ramaswamy H. Sarma


Introduction
Shigella are a group of Gram-negative, nonmotile, facultatively anaerobic, non-spore-forming rods bacteria (Paradh, 2015). Shigella can be distinguished from the closely associated Escherichia coli on the basis of pathogenicity, composition, physiology (inability to ferment lactose or decarboxylate lysine) and serology (Hale & Keusch, 1996;Jacewicz et al., 1986;Van den Beld & Reubsaet, 2012). The genus is classified into four serogroups with numerous serotypes: A (S. dysenteriae, 12 serotypes); B (S. flexneri, 6 serotypes); C (S. boydii, 18 serotypes); and D (S. sonnei, 1 serotype). It is a bacterial type that can cause severe diarrhea, generally in children (Agtini et al., 2005;Riddle et al., 2006;Zafar et al., 2005). S. dysenteriae type 1 is responsible for the most severe form of bacterial dysentery, which is considered as epidemics in various developing countries (Faruque et al., 2003;Pazhani et al., 2008;Taneja & Mewara, 2016;Tuttle et al., 1995;von Seidlein et al., 2006). S. dysenteriae type 1 can habitually develop to the most severe form of dysentery with deadly complications. Usually, the death is caused by severe colitis and are instantaneously accompanied with septicemia and pneumonia (Bennish et al., 1990;von Seidlein et al., 2006). The major source of this epidemic spread is found to be Shigella contaminated food and drink. A recent study also specified that Shigella can endure and persist in surface waters (Jatt et al., 2018). Possible symptoms of shigellosis include abdominal pain, tenesmus, watery diarrhea, and/or dysentery (multiple scanty, bloody, mucoid stools). Others may comprise abdominal tenderness, fever, vomiting, dehydration, and convulsions (Baveja, 2014;Niyogi, 2005). Infection usually began with ingestion of Shigella (usually via fecal-oral route) followed by an early symptom, diarrhea (possibly elicited by enterotoxins and/or cytotoxin), may arise as the organisms move through the small intestine (Desta, 2010;Girma, 2015). The trademarks of shigellosis are bacterial invasion of the colonic epithelium and inflammatory colitis. Both are considered to be interdependent processes which are amplified by resident release of cytokines and by the infiltration of inflammatory substances. Shigellosis can be appropriately identified in utmost patients on the presence of fresh blood in the stool, additionally, neutrophils in fecal smears. The most effective control strategy is the prevention of fecal-oral route of transmission (Keusch, 2009;Mehta et al., 1986;Speelman et al., 1987). More severe dysentery is treated with ampicillin, trimethoprim-sulfamethoxazole, or in patients over 17 years old, a 4-fluoroquinolone such as ciprofloxacin. Development of promising candidates against this disease is still ongoing.
Protein inhibition is being used for the treatment of several pathogenic infections. DHFR, is an enzyme that is found in S. dysenteriae which reduces dihydrofolic acid to tetrahydrofolic acid, using nicotinamide adenine dinucleotide phosphate (NADPH) as electron donor, where DHFR-1 is characterized specifically in Shigella. Indeed, NADPH can be converted to different varieties of tetrahydrofolate cofactors that are mandatory for DNA synthesis (S. K. Bhattacharya & Sur, 2003;Nelson et al., 1976). And, tetrahydrofuran (THF) is a vital cofactor involved in the transmission of methyl, methylene, and formyl groups from one molecule to another throughout the formation of nucleotides and several amino acids (Hawser et al., 2006). A high-level resistance to trimethoprim in S. dysenteriae is associated with plasmidencoded DHFR-1 (Moore, 1975;M. Navia et al., 2005). Tetrahydrofolate (H 4 -folate) is the key C 1 carrier in the synthesis of purines, thymidine, glycine, methionine, and pantothenate in bacteria and eukaryotes. In bacteria, H 4 -folate is also essential for the synthesis of formyl methionyl tRNA fMet (Moore, 1975). Dihydrofolate (H 2 -folate) consists of dihydropterin connected to p-aminobenzoate and to one or more glutamate residues that are connected to the p-aminobenzoate moiety (Chanarin et al., 1992;Heikkil€ a et al., 1990;Huovinen et al., 1995;Maden, 2000). In most bacteria and eukaryotes, the reduction of H 2 -folate to H 4 -folate is performed by the enzyme DHFR (Brown, 1971).
The conventional long duration antibiotic therapies available for Shigellosis possess various off-target effects and toxicity concerns. Widely used medications like Azithromycin and sulfonamides are found ineffective in several cases against Shigella universally. From 2006 to 2011, nearly half of the patients recovered in the Bay of Bengal which were reported as a resistant to Azithromycin (Orsomando et al., 2006). The sulfonamides resistance in Shigella affected patients in Japan were found to be increased to 90% in a duration of five years (D. D. Bhattacharya et al., 2014;Chen et al., 1994;S. Mitsuhashi et al., 1960). Henceforth there is a consistently expanding demand to investigate new compounds given the underlying emergency of emerging antibiotic resistance and increase in adverse effects. Exploring natural compounds, specifically phytocompounds aids by offering least toxic effects and novel target binding mechanisms. Drug repurposing against pathogenic organisms is solely significant as these organisms are resistant against almost every possible drug available in the market nowadays. Phytochemicals are screened with a major goal of inhibiting the pathogenic protein of the organism with minimal side effects. Previous studies suggest that Inula britannica which is a Eurasian species of plant in the genus Inula is found to comprise numerous phytocompounds that can be used to design lead molecules for a drug (Mitsuhashi, 1971) (Kim et al., 2002;Marberg et al., 1958;Song et al., 2000).
Inula britannica, commonly called British elecampane, British yellowhead or meadow fleabane, is an erect, rhizomatous, sunflower-like, herbaceous biennial or perennial in the composite family. This plant is innate to Europe and Asia where it naturally flourishes in a variety of moist environments including marshes, meadows, ditches, stream banks, wet grassland, wet woods, roadsides and wastelands (M. M. Rafi et al., 2005). Inula britannica L., family Asteraceae, is used in traditional Chinese and Kampo Medicines for various diseases. Flowers or the aerial parts are a rich source of secondary metabolites. These consist mainly of terpenoids (sesquiterpene lactones and dimmers, diterpenes and triterpenoids) and flavonoids, listed with their properties in Table  2. Such Inula derived compounds have shown diverse biological activities: anticancer, antioxidant, anti-inflammatory, neuroprotective and hepatoprotective activities (Das et al., 2020;Fischedick et al., 2013). Using in silico approaches such as molecular docking, ADMET and MD simulations (Figure 1), to study the interactions between the phytomolecules and the protein target, the enzyme, dihydrofolate reductase and predict the draggability of selected compounds 2. Computational methods

Sequence retrieval and physicochemical characterization
The FASTA sequence DHFR type I (DHFR-1) was retrieved from NCBI and the sequence was saved. Template search was performed using SWISS MODEL "Template Search" option which was based on basic local alignment search tool (BLAST) and HHblits which was an accurate and sensitive process (Waterhouse et al., 2018). The template was selected by evaluating the generated templates based on the sequence identity and global model quality estimate (GMQE). Physicochemical characterization of the retrieved sequences was performed using in silico ExPASy-ProtParam tool (http:// web.expasy.org/protparam/) that firms the physical and chemical parameters of the protein sequences. Different parameters, viz., number of amino acids, theoretical isoelectric point (pI), grand average of hydropathicity (GRAVY), molecular weight, aliphatic index and instability index, were computed (Ramsby & Makowski, 2005).    Secondary Structure Prediction by GOR IV Secondary Structure Prediction Method was used to predict the secondary structure of the amino acid sequence which is a Garnier-Osguthorpe-Robson method (Garnier et al., 1996) (https:// npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/ npsa_gor4.html). The three-dimensional structure of the target protein was determined using SWISS-MODEL (https://swissmodel.expasy.org/assess). SWISS-MODEL is a fully automated server that works based on the principles of homology modeling. Open Structure computational structural biology framework and the ProMod3 modelling is exploited by SWISS MODEL to perform this process accompanied by quality estimation (Waterhouse et al., 2018). The modelled protein is downloaded and saved in the .pdb format. The modeled structure was energy minimized using UCSF Chimera, an application that reduces the energy of the protein which helps in fixing the bond angles, and distorted protein geometrical arrangements (Pettersen et al., 2004).

Structure refinement and validation
GalaxyRefine2 is used to refine the protein using stronger harmonic restraints and the modelled and refined structures are validated and analyzed using Swiss structure assessment (https://swissmodel.expasy.org/assess). The Ramachandran plots of the modelled and refined structure are plotted. In regions where information is missing, structures predicted based on template structures or related sequences need to be improved further. It is critical to refine a predicted protein structure with particularly insufficient information on certain regions because these regions may be linked to functional specificity that is not conserved amongst associated proteins. It is critical to refine a predicted protein three-dimensional conformation with insufficient data on certain regions because these regions may be linked to functional specificity that is not conserved among related proteins. The GalaxyRefine2 (http://galaxy.seoklab.org/refine2) is an improved version of the GalaxyRefine protein assembly refinement server that incorporates latest advancements that have been effectively tested using CASP prediction experiments.

Active site prediction
Binding pockets of the protein was predicted using Computed Atlas of Surface Topography of proteins (CASTp) in which pockets are empty concavities on a protein surface into which solvent (probe sphere 1.4 A) can gain access, i.e. these concavities have mouth openings connecting their interior with the outside bulk solution (Tian et al., 2018). The CASTp web server objects to provide a complete and exhaustive computable characterization of protein interior voids and pockets that are present on its surface, which are noticeable concave areas of proteins that are often related with binding events. CASTp is built on the alpha shape and pocket process, both of which were developed in computational geometry.

Ligand selection
The twenty-nine different compounds like flavonoids, terpenoids, sesquiterpene and triterpenoids were selected based on their biological activities, molecular weights and the three-dimensional structures were downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/).

Docking studies
We utilised Autodock Vina to carry out the docking studies with the receptor and the ligand molecules and prepared them using MGL tools; when compared to the molecular docking software AutoDock 4, Autodock Vina utilises multithreading on multi-core machines and achieves approximately two orders of magnitude speed-up and a superior accuracy of the predicted binding conformations (Trott & Olson, 2010). Using the refined protein structure, the receptor was prepared by the deletion of water molecules, addition of polar hydrogen and Kollman's charges. The ligands .sdf structure was converted to .pdb structure using PyMol  (Lill & Danielson, 2011) and were then set for docking using MGL tools and there. pdbqt files were saved.
Upon saving the receptor.pdbqt file, the grid box was set for blind docking and the grid file was saved. The grid box was set at x ¼ 21.717, y ¼ 34.267 and z ¼ 22.920 and the output.pdbqt files were saved and analysed for the proteinligand interactions. The output files were visualised at first using PyMol, converted to .pdb files and visualised using Discovery studio (https://discover.3ds.com/discovery-studiovisualizer-download) for binding interactions.  . Ramachandran plot of protein structure before (a) and after refinement (b) process reducing the residues in disallowed region to zero and increase the residues in most favored regions to 92.6% from 86.1%.

ADMET filtering
SWISS-ADME (http://www.swissadme.ch/) was used to examine and filter the molecules based on their adsorption, distribution, metabolism and excretion (ADME) (M. O. Rafi et al., 2020) using the smiles obtained from PubChem. These parameters were used to identify their drug likeness and toxic nature. The compounds were strictly screened based on Lipinski's , Ghose's (Ghose et al., 1999), Veber's (Veber et al., 2002), Egan's (Egan et al., 2000) and Muegge's rules (Muegge, 2002), Bioavailability Score (Martin, 2005), PAINS (Baell & Holloway, 2010) and Lead Likeness (Teague et al., 1999). Any compound that did not satisfy these set of rules were exempt from performing further experimental studies due to their limitations in being an ideal drug-like molecule. And toxicity (mutagenicity and carcinogenicity) of nominated compounds were evaluated via using Discovery studio (Biovia, 2013).

Molecular dynamics (MD) simulation
MD simulation is now considered an authoritative step in the computational investigations of drug discovery at atomic level. Many mysterious biological functions in proteins and their profound dynamic mechanisms can be revealed by studying their internal motions (Al-Khafaji & Taskin Tok, 2021a). It exhibits the dynamical changes with the time scale where we can evaluate whether the protein-ligand complex is stable or not (Aldahham et al., 2020). In the present study, we subjected 4 ligands obtained from molecular docking to 100-ns timescale MD simulations of ligand-free and with ligand-bound forms. We utilized GROMACS 2018.1 package (Abraham et al., 2015) to run MD simulations. The protein-ligand complexes parameterized by using CHARMM27 force field for all atoms (Brooks et al., 1983). The three-point transferable intermolecular potential (TIP3P) was chosen as solvent (Al-Khafaji & Taskin Tok, 2020). And the charge was adjusted by adding Na þ or Cl À ions by imitating the environment of physiological systems. These systems were then energy minimized by using the steepest descent algorithm at the tolerance value of 1000 kJ/mol.nm (Al-Khafaji & Taskin Tok, 2021b). In the next step, NVT and NPT ensembles counterpoised the complexes with position restraint on the protein molecules for 0.1 ns. The Particle Mesh Ewald (PME) (Essmann et al., 1995) was appointed for dealing with the long range of Coulomb interactions through a Fourier grid spacing of 0.12 nm. While the short-range van der Waals interactions were calibrated through the Lennard-Jones potential with a cut-off distance of 1 nm. Now the MD simulations were carried out without any restraint on the protein molecules or ligand to determine the stability for 100 ns. Finally, some of Gromacs modules were used to analyze the MD trajectories such as gmx rms, gmx rmsf, gmx gyrate, gmx hbond. And the PCA estimations were calculated via gmx anaeig and gmx covar (Rafi et al., 2020).

Sequence retrieval and physicochemical characterization of the target
The protein sequence of interest was retrieved from NCBI which is a 136-residue length sequence taken from S. dysenteriae, strain1617 (https://www.ncbi.nlm.nih.gov/protein/ AHA64897.1), GenBank accession code AHA64897.1, named DHFR. DHFR is an enzyme responsible for reduction of dihydrofolic acid to tetrahydrofolic acid using NADPH which is inter-convertible to varieties of molecules that are essential for DNA synthesis like tetrahydrofolate cofactors.
Physicochemical characteristics like molecular weight, grand average of hydropathicity (GRAVY), theoretical pI, instability index, Coarse Packing Quality Control and aliphatic index were calculated. The total number of negatively charged residues (Asp þ Glu) and positively charged residues (Arg þ Lys) equal 16 (Supplementary Table S1). Protein instability index and aliphatic index gives confirmation regarding the solidity of the structure. The refined and energetically minimized structure had an instability index less than 40 while 93.31 aliphatic index gives relatively a high volume occupied alanine, valine, isoleucine, and leucine which might contribute to the thermostability of the protein.
A marginally negative value of GRAVY, the ratio of sum of hydropathy values and number of residues in the given sequence indicates moderate hydrophilic nature of the complex. Table 1 represents the values calculated for the physicochemical characteristics for the protein structure.

Protein structure prediction
With the template with highest sequence identity was selected to build the model of the protein, 46.27 sequence identity, the short-chain dehydrogenase from Pseudomonas syringae ((PDB ID:3GEM_A) was selected as the template for model building. The GOR (Garnier-Osguthorpe-Robson method) predicted Alpha helix (Hh) region of 65 residues (47.79%), 14 residues forming extended strand (Ee) (10.29%) and 57 residues forming random coil (Cc) (41.91%) as represented in Figure 2.

Protein structure refinement and validation
The modelled protein energy minimization and refinement was performed and the structures were comparatively validated. Figure 3(a) and (b) discusses the validation parameters and it was found that the refinement process helped to reduce the residues in disallowed regions to zero and increase the residues in most favored regions to 92.6% from 86.1%. MolProbity score was reduced from 2.34 to 1.76 for the refined protein while clash score was markedly reduced from 40.52 to 9.0. This helped to make the complexes Figure 6. Interaction diagrams generated using Discovery studio Visualizer of the protein DHFR-1 enzyme and ligand, spinacetin (a) is exhibiting the formation of two hydrogen bonds between the ligand and Phe33 and Asn13 residues of the protein and Pro9, Ilu30 and Leu12 forming pi-alkyl bonds, eupatin (b) is forming a hydrogen bond with Phe33 and five pi-alkyl bonds with residues Leu12, Leu 16, Leu 20, Ile30 and Pro9, chrysoeriol (c) is forming two hydrogen bonds with Phe33 and His5 residues and five pi-alkyl interactions with Pro9, Leu12 and Ile30, diosmetin (d) is establishing two bonds with Asn13 and Lys52 residues, interacting with Pro9 and Leu12 via pi-alkyl residues and is contacting to Phe33 via pi-pi T-shaped. energetically stable with minimum steric clashes. Comparing the modelled protein and the refined protein based on the Qualitative Model Energy ANalysis (QMEAN) score which explains better structural models, refinement decreases the QMEAN4 normalized score which should ideally fall into [0, 1] which was strictly followed by both the structures. Figure 4 explains detailed analysis of quality estimate based on QMEAN Z score, local quality estimate and comparison with non-redundant set of PDB structures.

Active site prediction
Active site prediction was performed to find the voids and pockets present on the surface of the protein. The protein 3 D structure analysis using CASTp showed that the modelled protein structure has 17 binding pockets. Each binding pocket has amino acid residues that can interact with the ligands. From the docking analysis, we came to the conclusion that among the 17 binding pockets, we have chosen the 1 st ranked binding site based on the area and volume of the predicted active sites.
Hence, the major binding pocket under study is marked in the green pocket as shown in Figure 5. One of the striking features of these pockets is that they all fall into the surface of the protein which implies that the ligands, importantly bind to the surface of the receptor.

Molecular docking, ADME analysis and visualization
Compounds were selected based on their previously identified properties and this consisted of 29 compounds of the plant Inula Britannica. Consisting of seven sesquiterpenes, five triterpenoids, and seventeen flavonoids. Molecular docking is used to evaluate the potential of several selected compounds to treat Shigella dysenteriae type 1 via targeting DHFR-1. Using Autodock Vine, the binding mode and free energy binding mode of the selected compound with DHFR-1 were estimated and compared with natural substrate of DHFR-1 (dihydrofolic acid) ( Table 2). Docking studies showed that the docking score of dihydrofolic acid (i.e. natural substrate) was À5.9 kcal/mol while the docking scores of 29 compounds have various binding affinities toward DHFR protein which are ranging from À8.8 kcal/mol (B-Amyrin) to À5.5 kcal/ mol (6b-O-(2-Methylbutyryl) Britannilactone). Depending on the docking scores' cpmparsion, 1,6-O, O-Diacetylbritannilactone, O-Acetylbritaanilactone, Neobritannilactone B and 6b-O-(2-Methylbutyryl) Britannilactone were filtered out due to having equal or lower docking scores than dihydrofolic acid's docking score. The rest 25 compounds were subjected to AMDET studies.
Evaluating the efficacy or risks of these small molecules, the properties of ADMET have become one of the most important aspects of biology (Singh et al., 2017). Therefore, as part of drug discovering, ADMET estimation is a crucial. Drug metabolism and pharmacokinetics (DMPK) studies, commonly referred to as ADMET is carrying out and it is the most important and challenging step in drug discovery and development (Khan et al., 2010). The high biological activity and low toxicity of the drug or drug-like compound are not enough to make the compound a good drug candidate. All of optimal pharmacokinetic characteristics are essential for the process of discovering the drugs. Therefore, it is very important to evaluate the ADMET profile to avoid wasting time/resources. Therefore, we have used the swiss-ADME online software to predict the ADME characteristics and Tox of the selected compounds was predicted by Discovery studio via evaluating, Ames toxicity and carcinogenicity.
The Lipinski violations are one of the number one standard to test the drug likeness of any compound if it had been to be administered orally. The five rules recommend that a molecule with the best traits for use as a drug should have the subsequent regulations: following regulations: molecular weight should be less than or same to 500 DA, number of hydrogen bond donors should be less or equal to 5, hydrogen bond acceptors should be less or equal to 10, and XlogP3 should be less or equal to 5. Thus, it is imperative that a compound considered for drug development, should have least violations to this rule (Lipinski et al., 2001). The Ghose filter determines the drug likeliness of a molecule through the subsequent constraints, a log p value should range between À0.4 and 5.6, the molecular weight ought to be among a 160 and 480, a molecular refractivity (MR) among 40 and ve130, and the whole number of atoms ought to be among 20 to 70 (Ghose et al., 1999). Veber's rule (Veber et al., 2002) decrees that vital standards ought to be met for a molecule to be suitable for oral bioavailability. First, the range of hydrogen bonds withinside the molecule ought to now no longer exceed ten, and second, the polar floor place ought to now no longer exceed 140Å 2 , which in Figure 8. PCA. (a) projection of the motion for apo DHFR-1 (black) and spinacetin-DHFR-1 (red) in phase space, (b) projection of the motion for apo DHFR-1 (black) and eupatin-DHFR-1 (green) in phase space, (c) projection of the motion for apo DHFR-1 (black) and chrysoeriol-DHFR-1 (blue) in phase space and (d) projection of the motion for apo DHFR-1 (black) and diosmetin-DHFR-1 (yellow) in phase space. (e) Eigenvalues for the first fifteen eigenvectors apo DHFR-1 (black), spinacetin-DHFR-1 (red), eupatin-DHFR-1 (green), chrysoeriol-DHFR-1 (blue) and diosmetin-DHFR-1 (yellow).
the end corresponds to the molecule having much less than twelve hydrogen donors and acceptors. Egan's rule states that a molecule has proper oral bioavailability in the event that they fulfil the log p value withinside the variety of À1.0 and 5.8 and a topological polar surface area (TPSA) value much less than or equal to 130Å 2 (Egan et al., 2000). Muegge's rule (Muegge, 2002) edicts for a molecule's draggability are that the molecular weight need to be withinside the variety of 200 to 600, a lipophilicity profile (xlopP3) between À2 and 5, the TPSA much less than or same to 150Å 2 . Additionally, the variety of rings need to be much less than or same to seven, the variety of carbon atoms need to be greater than or same to four; the variety of hetero atoms greater than one. The variety of rotatable bonds need to be lesser than or same to fifteen. Finally, the variety of hydrogen bond acceptors need to be lesser than or same to ten, and hydrogen bond donors need to be lesser than or same to five. The Brenk and pan assay interference compounds (PAINS) structural indicators had been utilized in medicinal chemistry to expect unstable, reactive, poisonous fragments gift withinside the structure. have 0 indicators in Brenk and PAINS descriptors, offering different promising signs to be drug candidates (Singh et al., 2017).
Pharmacokinetic profile is very important in drug discovery. Gastrointestinal absorption and entry into the brain are two pharmacokinetic behaviours. Therefore, four drugs have successfully passed the drug likelihood test, spinacetin, eupatin, chrysoeriol and diosmetin were well absorbed from the gastrointestinal tract, but cannot be transferred through the BBB or be a P-glycoprotein substrate. From these accumulated data, it can be concluded that the route of administration includes oral and sublingual (the drug dissolves under the tongue) (Metwally et al., 2019).
Toxicity profiling of all docked compounds was also performed by implementing toxicity prediction extensible protocol of discovery studio 3.5. Toxicity profile involves AMES mutagenicity and carcinogenicity. Based on the aforementioned characteristics, we obtained only diosmetin, chrysoeriol, spinacetin and eupatin as good and safe drugs for Shigellosis.

Conformational of protein
Binding of a compound in the binding site of the protein can lead to observable conformational changes in the dynamics of targeted protein. Root mean square deviation (RMSD) is one of the most important fundamental properties for establishing whether the protein is stable and close to the experimental structure (Al-Khafaji & Taskin Tok, 2020). The RMSD plot suggests that the binding of eupatin and chrysoeriol kept the dynamics of DHFR-1 while the binding of spinacetin and diosmetin led to more structural deviations from its native conformation (Figure 7(a)). In the case of the apo DHFR-1 and chrysoeriol-DHFR-1 complex, the nature of dynamics of them were having the same style during the 100 ns of MD simulation and in less degree the complex of eupatin-DHFR-1 (Figure 7(a)). Moreover, the dynamics of understudied drugs inside the active site were compared and presented in Figure 7(b). It can be observed that all of these ligands have nearly the same nature of movements inside the active site.
To investigate the dynamics of the protein's backbone residues in the protein-ligand complexes compared to the apo form of protein, the root mean square fluctuations (RMSF) of the backbone atoms of the protein were depicted in Figure 7(c). The figure reveals that the RMSF values of eupatin and chrysoeriol when complexed with protein reduced to be lower than the RMSF of apo form of protein.
While spinacetin and diosmetin make the residual fluctuation of DHFR-1 more elevated than apo form. Thus, RMSF variation shows that eupatin and chrysoeriol stabilize the protein, which could indicate more efficiency in blocking.

Radius of gyration (Rg)
Radius of gyration (Rg) factor refers to the compactness of protein during the time scale of molecular simulation. It is intelligible that it is used to measure the distance between the center of mass of the protein atoms and its terminal in a given time step. Generally, a stable protein tends to preserve a relatively less fluctuation in Rg value which is a determinant factor in evaluation of dynamic stability. In current work, the Rg values are depicted in Figure 7(d). The radius of gyration averages was 1.724897329, 1.791844007, 1.731573548, 1.682601293, and 1.758484603 for apo DHFR-1, spinacetin-DHFR-1, eupatin-DHFR-1, chrysoeriol-DHFR-1 and diosmetin-DHFR-1, respectively. This suggests that compactness of chrysoeriol-DHFR-1 is lower than apo form of apo and eupatin-DHFR-1 is comparable to the apo DHFR-1.

PCA
The plot in Figure 8(a,b,c and d) shows the 2 D projection of the trajectories for two major principal components PC1 and PC3 for apo DHFR-1 compared with spinacetin-DHFR-1, eupatin-DHFR-1, chrysoeriol-DHFR-1 and diosmetin-DHFR-1. It is evident from the 2 D plot that the protein-eupatin complex and chrysoeriol-DHFR-1 showed higher stability and occupies lesser phase space along with well-defined clusters compared to spinacetin-DHFR-, diosmetin-DHFR-1 and apo DHFR-1, indicating a more stable complex formation, showing the least dynamics and stable conformational space.

Hydrogen bond analysis
Hydrogen bond number and distribution in the protein-spinacetin and protein-eupatin complexes were studied to determine the stability of spinacetin and eupatin inside the binding site during the 100 ns simulation period. Hydrogen bond numbers results show that spinacetin and eupatin formed stable interactions with protein Figure 9(a). The average of hydrogen bonding between DHFR-1 with spinacetin, eupatin, chrysoberyl and diosmetin were 0.414758524, 0.879612039, 0.872412759 and 0.369063094, respectively.
Besides, the results of the distribution of hydrogen bond lengths (Figure 9b) also indicates that the spinacetin, eupatin, chrysoberyl and diosmetin have the ability to establish hydrogen bonds with different lengths. The hydrogen bond results help in understanding the functionality and ability of eupatin and chrysoberyl to efficiently hinder.

Discussion
Inula britannica belongs to the family Asteraceae and is widely used by the traditional Chinese and Kampo Medicines for various diseases (Khan et al., 2010). Secondary metabolites are vastly found in the flowers or the aerial parts (Singh et al., 2017). Secondary metabolites mostly consist of terpenoids like sesquiterpene lactones and dimmers, diterpenes, triterpenoids and flavonoids (Khan et al., 2010). Various biological activities such as anticancer, antioxidant, anti-inflammatory, neuroprotective and hepatoprotective activities were observed from the isolated compounds (Song et al., 2002). This relates that the compounds of this plant are promising for the designing of lead molecules against S. dysenteriae. Reduction of prevalence and death rates of the disease can be overcome by the Antibiotic treatment of Shigella infections (Kim et al., 2002). The following antibiotics Ampicillin, amoxicillin, Nalidixic acid, ciprofloxacin, norfloxacin, Azithromycin and Trimethoprim are extensively used to treat Shigella dysentery (Bai et al., 2006). In the recent studies, S. dysenteriae shows resistance to the number of antibiotics available in the market. This is due to the worldwide rise of broad-spectrum resistance against many antibiotics (Kotloff et al., 2018). Drug resistance of Shigella spp. can be a result of many mechanisms, like decrease in cellular permeability, extrusion of drugs with the active efflux pumps, and overexpression of drug-modifying and the inactivating enzymes or target modification by mutation (Murphy et al., 1993). There is an increasing need for the identification and evolution of alternative therapeutic strategies. The most common mechanism of trimethoprim (TMP)-resistance is the acquisition of dihydrofolate reductase enzyme resistant to this drug (Blair et al., 2015). Therefore, DHFR-1 is considered as an important target to design new anti-resistant and antibacterial drugs. DHFR-1 is a very important enzyme because it produces cofactors that are necessary for DNA synthesis. Specifically, DHFR-1 catalyses the reduction of folate and 7, 8-dihydrofolate (DHF) to 5,6,7,8-tetrahydrofolate (THF) (Fern andez & Hancock, 2013). The transfer of methyl, methylene, and formyl groups from one molecule to another, during the production of nucleotides and several amino acids is done by an essential cofactor THF. Utilization of carbon units from a THF cofactor by thymidylate synthase to make thymidine nucleotides is an example of the activity of THF (M. Navia et al., 2005). Therefore, by targeting DHFR-1 protein, viral replication can be prevented (M. M. Navia et al., 2003). Phytomedicines can become a great alternative to chemical compound-based drugs. With the advancement of bioinformatics by using computer aided drug discovery, the identification and virtual screening of compounds can be done easily. 29 different compounds like flavonoids, terpenoids, sesquiterpene from Inula britannica were taken into consideration. The phytocompounds that successfully docked well into the binding site of receptor protein and passed from ADMET filtration are: eupatin (-6.5 kcal/mol), diosmetin (-6.5 kcal/mol), chrysoeriol (-6.3 kcal/mol) and spinacetin (-6.1 kcal/mol) with DHFR-1. It is evident that four flavonoid compounds namely eupatin, diosmetin, chrysoeriol and spinacetin are effective compounds for treating shigellosis due to they have diverse of biological activities (Choi et al., 2005;Chougouo et al., 2016;Miller & Ruiz-Larrea, 2002;Santiago-Saenz et al., 2019). Flavonoids are well known as antibacterial agents against a wide range of pathogenic microorganisms (Liao et al., 2003). Similar In silco study was conducted with flavonoids of Inula Britannica which weakened the toxicity induced by glutamate in primary cultures of rat cortical cells. Another recent study was also conducted to study the antiinflammatory and anti-cancer effect of 1 b hydroxy alantolactone, a derivative from Inula Britannica.
With increasing prevalence of untreatable infections induced by antibiotic resistance bacteria, flavonoids have attracted much interest because of the potential to be substitutes for antibiotics. Inhibition of nucleic acid synthesis, inhibition of cytoplasmic membrane function, inhibition of energy metabolism, inhibition of the attachment and biofilm formation, inhibition of the porin on the cell membrane, alteration of the membrane permeability, and attenuation of the pathogenicity are the determined antibacterial mechanisms of flavonoids (Miranda et al., 2016). Spinacetin, eupatin, chrysoeriol and diosmetin compounds were selected for analysing their binding energy from Autodock simulations followed by their screening based on ADMET properties.
In silico docking and toxicity analysis suggests that the compounds that successfully interacted with the target protein with good binding affinity and favourable ADMET properties are safe for their extensive use as a therapeutic agent and also molecular dynamics (MD) simulations are also done which played a foremost role in these attempts to advance docking procedures (Xie et al., 2015) but due to lockdown we were unable to carry out their in silico and in vivo validations and effectiveness of these compounds. However, we understand that this finding may set forth an idea to use these phytocompounds for the researchers who have been working to produce anti-resistant antibiotics against resistant pathogenic bacteria.

Conclusions
Inula britannica phytocompounds have significant biological implications against different diseases. Moreover, inhibition of DHFR-1 plays a key role in reducing the danger of S. dysenteriae type 1. Phytocompounds screened at two levels of molecular docking followed by ADME studies and molecular dynamics simulation to investigate the stability of the successful candidate with target protein. Eupatin and Spinacetin showed good docking affinities and good pharmacokinetic properties. While Molecular Dynamic simulation showed promising results about eupatin and this is in agreement with docking results as discussed in detail. Hence, it is concluded that eupatin and chrysoeriol would be a crucial biomolecule for the treatment of shigellosis caused by S. dysenteriae type 1 after appropriate in vitro and in vivo studies.

Disclosure statement
No potential conflict of interest was reported by the authors.