In silico molecular docking study of Andrographis paniculata phytochemicals against TNF-α as a potent anti-rheumatoid drug

Abstract Tumor necrosis factor-α (TNF-α) is a proinflammatory cytokine which plays a crucial role in controlling inflammatory responses. The pathway of Rheumatoid arthritis (RA) leading to TNF-alpha is activated by macrophages and quite often by natural killer cells and lymphocytes. In the inflammatory phase, it is believed to be the main mediator and to be anchored with the progression of different diseases such as ankylosing spondylitis, Crohn's disease, and Rheumatoid arthritis (RA). The major goal of this study is to use in silico docking studies to investigate the anti-inflammatory potential of a bioactive molecule from the medicinal plant Andrographis paniculata. The three-dimensional structures of different phytochemicals of A. paniculata were obtained from PubChem database, and the receptor protein was derived from PDB database. Docking analysis was executed using AutoDock vina, and the binding energies were compared. Bisandrographolide A and Andrographidine C revealed the highest score of −8.6 Kcal/mol, followed by, Neoandrographolide (−8.5 Kcal/mol). ADME and toxicity parameters were evaluated for these high scoring ligands and results showed that Andrographidine C could be a potent drug, whereas Neoandrographolide and Bisandrographolide A can be modified in in vitro and can lead to a promising drug. Further, the top scorer (Andrographidine C) and control drug (Leflunomide) were subjected to 100 ns MD Simulation. The protein complex with Andrographidine C had more stable confirmation with lower RMSD (0.28 nm) and higher binding energy (−133.927 +/− 13.866 kJ/mol). In conclusion, Andrographidine C may be a potent surrogate to the disease-modifying anti-rheumatic drugs (DMARD’s) & Non-steroidal anti-inflammatory drugs (NSAID’s) that has fewer or minor adverse effects and can aid in RA management. Graphical Abstract


Introduction
Rheumatoid arthritis (RA) is a heterogeneous chronic inflammatory disease affecting the small synovial joints and characterized by synovium hyperplasia and accumulation of inflammatory cells, resulting in progressive joint destruction that amasses over the years, resulting in severe impairment, morbidity and increased mortality among patients (Grove et al., 2001). In RA, several cytokines such as TNF-a (Tumor necrosis factor -a), TNF-b (Tumor necrosis factor -b), Interleukins (IL À 1 to 6) are involved in particular destruction. TNF-a is a 17 kDa monomeric mature molecule that forms a biologically active homotrimer non-covalently (Mancini et al., 1999). TNF-a is activated once it interacts with receptor, TNFR1 and TNFR2, which are well-known targets in the treatment of autoimmune diseases. The main function of TNFR1 is to trigger TNF-a signaling pathway which causes inflammation that is found to be associated with diseases such as arthritis and hepatitis (Mukai et al., 2010). TNF's signaling cascade is strongly related with TNFR1 interaction, thus initiating several pathways such as mitogen activated protein kinase (MAPK), extracellular signal regulated kinase (ERK), activation of c-Jun N-terminal kinase and many such pathways (Urschel & Cicha, 2015). TNF-a being a core activator of several signaling mechanisms, this protein has been greatly focused in search of small molecule inhibitors for solving immunological dysregulations. Following to this antagonist for TNF-a have been revolutionized in therapeutic management for several autoimmune diseases including Rheumatoid Arthritis. As a result, suitable inhibitors can be chosen to block TNF-alpha by stabilizing the inactive dimeric form, preventing active trimeric TNF-alpha from attaching to its receptors and shutting down the inflammatory signaling pathway which leads to RA (Choi et al., 2010;Wajant et al., 2003).
For the medicament of rheumatoid arthritis, many biologics targeting TNF-have been approved. In 1995, Suramin, a symmetrical polysulfonated urea derivative, was reported to inhibit the interaction of TNF-a with its cellular receptors by promoting the dissociation of trimeric human tumor necrosis factor-a (TNF-a) into physiologically inactive subunits (Alzani et al., 1995). The use of disease-modifying anti-rheumatoid medications (DMARDs) is the mainstay of RA treatment that comprisemethotrexate (MTX) (Hazlewood et al., 2016), sulfasalazine (SSP) (Dale et al., 2007) and leflunomide (O'Donnell et al., 2010), Non-steroidal anti-inflammatory drugs (NSAID's) and modern therapy including b-DMARDs and ts-DMARDs.
These are helpful in postponing the onset of illnesses, including joint deterioration, but they also have particular adverse effects which involve increased toxicity, bacterial and viral infections, cancer, hematologic problems, and cardiovascular diseases. (Antoni & Braun, 2002;Nell et al., 2004;Wilsdon & Hill, 2017).
In order to obtain a potent natural drug targeting the desired protein, a computer-aided drug design (CADD) or an in silico molecular docking analysis is a low-cost and rapid solution. It is widely acknowledged that medication development and manufacture are time-consuming and energyintensive operations. As a result, an increased effort is being made to adapt computing power to the combined chemical and biological domain in order to expedite the drug development process (Marshall, 1987). Computer-aided or in silico architecture is used in the biomedical arena to speed up and promote hit detection, hit-to-lead collection, optimize the profile of absorption, distribution, metabolism, excretion, and toxicity, and prevent security problems (Kapetanovic, 2008). Molecular docking examines how tiny molecules known as ligands interact with the target macromolecule's active binding pocket (receptor proteins). In addition to generating binding energies, the position of ligand in the enzyme binding site can be visualized in these docking studies. It may be useful for the development and understanding of the binding nature of potential drug candidates (Madeswaran et al., 2012).
Moreover, this ligand-based virtual screening method is also used to classify the most likely molecules with a pharmacological activity using biological interaction models to scan for compounds. Similarly, several algorithms exist for pharmacokinetics, toxicity, and drug-likeness prediction studies, making the job simpler. In the discovery of naturally derived medicines, there is much proof that proves the application of computational instruments (Kiran et al., 2020). Therefore, the purpose of the current study is to apply in-silico screening methodology to screen the bioactives from Andrographis paniculata for its anti-rheumatoid activity.
In the pipeline of the discovery of natural medicines, recent literature studies found that plant molecule contribute as a critical source of natural anti-inflammatory drugs (Arunkumar & Rajarajan, 2018). About 1300 plants have been discerned as medicinal plants and are in use as novel drugs . The plant A. paniculata, owned by the Acanthaceae family, is an essential medicinal plant in Unani and Ayurvedic medicines and is extensively used worldwide as traditional herbal therapeutics. The compounds present in the plant, i.e. Neoandrographolide, Andrographiside, Andrographolactone, Andrographidine, 14-Deoxy-11-hydroxyandrographolide, etc. accounts for its major biological activities such as anticancer (Cheung et al., 2005), antimalarial (Misra et al., 1992), antiviral (Thirumoorthy et al., 2021), antidiarrheal (Gupta et al., 1990), anti-HIV (Nanduri et al., 2004) & antihyperglycemic (Yu et al., 2003). The primary compound in the plant was reported to be Andrographolide, which has favourable effects on inflammation-related diseases, including infection with viruses, bacterial dysentery, malaria, and rheumatoid arthritis. In recent years, commercial production of these plant extract is also used (Dai et al., 2019). However, it is still essential to standardize the preparations for their better efficacy (Sanower Hossain et al., 2014).
Using molecular docking techniques, the current work used a combination of in silico methodologies to better understand the critical direct interactions between TNF-a and phytochemicals. Further, the ADME properties and druggability nature of the molecules were predicted using SWISS ADME, and pKCSM have been carried out to find the toxicity nature of the phytocompounds present in Andrographis paniculata. Moreover, Molecular dynamic simulation was carried out for the appropriate potent drug molecule along with the control drug.

Selection and preparation of ligands
The two-dimensional structure of the ligands which included the phytochemicals from the plant A. paniculata were retrieved from PubChem database in SDF format. Structural optimization which included protonation and energy minimization was performed prior to molecular docking using Avogadro software and was saved in PDB format.

Selection and preparation of receptor
TNF-a protein ( Figure 1) which is considered as receptor protein was downloaded from RCSB protein data bank PDB ID 2AZ5 that has a three-dimensional crystal structure having a small molecule inhibitor (SPD-305) that has been shown to block TNF action in biochemical and cell-based studies. The protein was prepared for molecular docking study using AutoDock MGL Tools by eliminating the inhibitor and water molecules. Further polar hydrogen bonds were added, missing atoms were checked and repaired, Kolaman and gangsteir charge were added to the protein using AutoDock tools and saved in PDBQT format whereas for swiss dock the PDBQT file was converted to PDB using PyMol.

Binding site prediction
The TNF-a (2AZ5), 3-D structure obtained from PDB has a cocrystalized ligand, which was utilized to predict the active site. Discovery studio visualizer was employed to designate the co-ordinates of the site where the ligand (small inhibitor molecule) was located (as x ¼ À18.127, y ¼ 74.486, z ¼ 35.744) that was further confirmed as active pocket of the receptor protein.

Molecular docking
Molecular docking is an approach to knit structure-based drug design for identifying the basic interactions of amino acids residues between the chosen low energy conformation protein and generated ligands. AutoDock Vina and Swiss dock was used for docking simulation of all ligand molecules. The ligand molecules were optimized with AutoDock tools and stored in PQBQT format (which is a modified protein data bank format for ligands that include atomic charges, atom type descriptors, and topological information) for Auto dock vina whereas for Swiss dock the ligands were optimized and converted to mol2 format.
Furthermore, both the tools were employed to test the relationship of the hydrogen bond and binding affinities. The grid box parameters were set at the x, y, and z (x ¼ À18.127, y ¼ 74.486, z ¼ 35.744) points based on active site prediction after the energy minimization process. The coordinates were saved in a configuration file for AutoDock vina and the same parameters were uploaded as active site coordinates onto the Swiss dock server. The docking was proceeded and further the visualization and analysis work of the molecules to identify and predict the interaction was carried out using Discovery Studio 2020 and PyMol.

Docking validation
The small inhibitor molecule co-crystallized TNF-a was used to validate the docking technique utilizing the pose selection method. Using maestro, the heteroatom was removed from the protein, and the protein and inhibitor crystal structures were stored as separate PDB files. The same procedure was followed throughout the operation, as well as the grid settings. This was done to ensure that the inhibitor correctly binds to the active site pocket and differs just slightly from the same co-crystallized complex. The re-docked complex was seen in PyMOL 2.3 by superimposing it over the reference co-crystallized TNF-alpha complex. The RMSD value was also computed, and comparable amino acid residues in the 2D picture plotted using Ligplot þ v2.2.4 were compared and highlighted.

Evaluation of drug-like properties
Using Swiss ADME website (http://www.swissadme.ch/) Adsorption, Distribution, Metabolism, Excretion properties of the selected top scorer ligands were anticipated. The phytochemical's canonical smile format was uploaded to test its pharmacological and drug similarity on the Swiss ADME platform. The tool calculates physiochemical descriptors and predicts ADME parameters, drug-like nature, and molecules that have been shown to endorse the discovery of drugs (Daina et al., 2017).

Toxicity
Toxicity assessment was done by using pkCSM website (http:// biosig.unimelb.edu.au/pkcsm/prediction). Smile format of selected ligands were uploaded and the different parameters like AMES toxicity, Max. tolerated dose for humans, hERG I & II inhibitors, oral rat acute toxicity, oral rat chronic toxicity, hepatotoxicity, skin sensitization, T. pyriformis toxicity and Minnow toxicity was observed and results were analyzed.

Molecular dynamics simulation
The MD simulation of the complex i.e. APO, Andrographidine C & Leflunomide were carried out on GROningen MAchine for Chemical Simulations (GROMACS) version 2019.4. The 'pdb2gmx' tool was used to create the protein topology, whereas the GlycoBioChem PRODRG1 server was used to create the ligand topologies for getting the force field coordinates. Using the steepest descent technique, the system was vacuum minimized for 1500 steps. The complex structures were then solvated using a water simple point charge (SPC) water model in a cubic periodic box of 0.5 nm. Following that, the complex systems were maintained, i.e. neutralized with an adequate salt concentration of 0.15 M by introducing the necessary quantities of Na þ and Cl À counter ions (Bhardwaj et al., 2021).
A steepest descent minimization approach with a maximum of 50,000 steps and a force of 10.0 kJ/mol was used to reduce the energy of the complexes. The equilibration of the system was accomplished in two stages. The first stage was to keep an NVT ensemble with a constant number of particles, volume, and temperature for 2 ns, and the second step was to keep an NPT ensemble with a constant number of particles, pressure, and temperature for 10 ns. MD position restrictions were imposed for 100ps to both ensembles (NVT and NPT). The backbone C-atoms were restricted in this phase, but all solvent molecules were free to move, ensuring that the system's solvent equilibrium was preserved. The linear constraint solver technique restricted the system's covalent bonds. The long-range electrostatic interactions were obtained using the particle mesh Eshwald method with a 1.2 nm cut-off and 1.2 nm Fourier spacing. In terms of shape and solvent orientation, additional simulation steps were undertaken on the well-equilibrated and solvated structures. The temperature of the system was regulated at 310 K using the V-rescale weak coupling approach. To equilibrate and establish the pressure (1 atm), density, and total energy of the system, the Parrinello Rahman method was employed. The well-equilibrated complexes were then placed through a production phase with no restrictions for 250 ns with a 2 fs time step, with structural coordinates recorded every 2 ps. The last production run in NPT ensemble for 100 ns simulation time was applied to each resulting structure from the NPT equilibration phase. The GROMACS simulation, suite of Protein RMSD, RMSF, RG, and H-Bond was used to perform trajectory analysis (Bhardwaj et al., 2021).

MM-PBSA
To study the binding free energy (G binding) of an inhibitor with protein throughout simulation time, a Molecular Mechanics Poisson-Boltzmann surface area (MM-PBSA) method was used. MD scripts were extracted in order to do MM-PBSA calculations. The binding free energy gives you a good idea of how protein and ligand interact biomolecularly. The binding energy is made up of potential energy, polar and non-polar solvation energies. The MM-PBSA binding free energy was calculated using GROMACS's 'g mmpbsa' script. We calculated G for the last 20 ns in using dt 1000 frames to get an accurate result. The binding free energy was estimated using the GROMACS function g mmpbsa, which is computed using the following formula: The DG binding represents the total binding energy of the complex, while the binding energy of free receptor is G receptor , and that of unbounded ligand is represented by G ligand .

Molecular docking analysis
Specifically, small molecule inhibitor bound site is examined for the molecular interaction of twenty phytochemicals of Andrographis paniculata against TNF-a protein through molecular docking approach using AutoDock vina and Swiss dock. AutoDock vina is an open-source programme that performs molecular docking analysis. Due to its accuracy, high precision programmes, user-friendliness and ease of access, it was found to be a strong competitor to other programmes (Trott & Olson, 2009). Alongside, Swissdocka free online server which generally use PDB format, it enables the uploading of CHARMM formatted data and gives option to upload protein structures as a collection of protein structure files (PSF), coordinate files (CRD), additional topology (RTF), and parameter files (PAR) (Grosdidier et al., 2011). In this study, the receptor protein TNF-a's crystal structure was obtained from the RCSB -PDB bank and the phytochemical component structures were downloaded from the PubChem database. Initial docking analysis exhibited that the bioactives have an excellent tendency with receptor protein (TNF-a) and showed lesser binding energy contrasted to the standard control drug (leflunomide). To confirm the docking expression, redocking was performed using Autodock vina and Swiss dock. Based on their docking poses the ranking log of docked ligands and their corresponding docking scores are represented in Table 1 and redock scores are mentioned in the supplementary Table  1. The dynamic nature of the compound is based on its lowest binding energy (Madeswaran et al., 2012). The docking results of the phytochemicals with reference drugs Leflunomide: À7.1 kcal/mol against TNF-a revealed that the molecules Bisandrographolide A (À8.6 kcal/mol), Neoandrographolide (À8.5 kcal/mol), and Andrographidine C (À8.6 kcal/mol) were adequately positioned in the binding pocket of the receptor protein (Prasanth et al., 2021).

Molecular interaction studies
Among the twenty compounds, top three scorers were taken for molecular visualization to identify the interaction involved with the ligand-protein complex, which quotes for its high binding energy utilizing Discovery studio and Protein-ligand interaction profiler (https://projects.biotec.tu-dresden.de/plipweb/plip/submitCalculationJob). The binding energies of the phytochemicals Andrographidine C and Bisandrographolide A are same i.e. À8.6 kcal/mol which is higher than the binding free energy of other compounds because of the most favourable interaction at the binding pocket surrounded by polar residues. To prevent the binding of the third subunit, which forms the physiologically active trimer complex, small-molecule inhibitors of TNF-a should be somewhat hydrophobic and big enough to interact with both subunits of the TNF-a dimer concurrently (Chan et al., 2010). Andrographidine C has polar residues Ser60 (2.06 Å˚), Leu120 (2.72 Å˚), Tyr151 (2.43 Å˚), Gly121 (3.50 Å˚); hydrophobic residues Leu57 (5.16 Å˚), Tyr59 (3.86 Å˚) where Tyr59 is one of the related hydrophobic residues that has been linked to trimer contact of TNF-alpha (He et al., 2005). Furthermore, Bisandrographolide A has polar residues Tyr151 (3.06 Å˚), Leu120 (3.16 Å˚); hydrophobic residues Tyr119 (3.92 Å˚), Tyr59 (4.43 Å˚). The critical interactions of the top scorer molecules and control drug against TNF-a at the binding site are shown in Table 2 (Biswas, 2019). This showed that natural compound-based medicines could be a significant alternative in the place of synthetic drugs.

Validation of docking
The cognate docking was done to see how docking technique is effective.. For the cognate docking procedure, we utilized the same docking parameters as before. The interacting amino acids in the co-crystallized structure of TNF-alpha are Leu57, Gly121, Tyr59, Leu120, Tyr119, Gly121, Tyr119, Tyr151, Ser60, Leu120. The docked complex was overlaid on the native TNF-alpha co-crystallized structure with PyMOL, and the amino acid residue similarity was visualized with LigPlot þ v.2.2.4. The complex was completely superimposed onto the original co-crystallized complex without any changes, as shown in the overlaid 2D structure in Figure 6. In both complexes, all of the amino acid atoms were superimposed with no limitations. The alignment score was 0.000, with a total of 10 amino acid residues overlaid, and re-docking of the natural ligand into the binding area of the protein structure verified the docking approach with an RMSD value of 0.007. As a consequence, this method demonstrated the docking procedure' effectiveness and validity.

Evaluation of drug-like properties
The compounds that come in the top-tier of screening were evaluated for physico-chemical properties using Swiss ADME (http://www.swissadme.ch/). All the molecules were submitted in SMILE format, and the results of the top scorer ligand's molecular properties and drug-likeliness were retrieved and represented in Table 3. Our results indicated that among those top scorer molecules, the phytochemical Andrographidine C possessed significant drug-likeliness. The physico-chemical parameters satisfied by the compound are Topological polar surface area (TPSA) À 148.05Å2, Partition coefficient between n-octanol and water (log Po/w) of the molecule -iLOGP (2.70), SILICOS-IT (1.44) and Consensus P0/W (0.92), Solubility property (soluble À2.52e-01mg/ml). Low gastrointestinal absorption (GI), no permeable bloodbrain barrier, operates as a P-gp substrate, and does not inhibit any of the P450 cytochrome molecules except CYP3A4 were identified in the anticipated pharmacokinetic data.    The drug likeliness of a compound is majorly determined by Lipinski's Rule of 5. In addition, Lipinski's rule aids in the prediction of whether the particular compound has specific biological activity that will make it a likely human orally active drug (Kiran et al., 2020). Hence the druggability factor was found to be with zero violation which obeys Lipinski's Rules and also obeys drug likeliness including Ghose, Muegge, Veber, Egan with 0.55 Bioavailability score. Similarly, Kaloni et al. (2020) successfully applied Swiss ADME prediction in their analysis to evaluate the drug likeliness property of Murraya koenigii phytochemicals.

Toxicity
In silico toxicity profiling of the selected ligands showed a very low toxicity profile. As shown in Table 4 there is an absence of AMES toxicity in all three ligands namely, Neoandrographolide, Bisandrographolide A and Andrographidine C. Maximum tolerated dose of the drug for human should be less than or equal to 0.477 (log mg/kg/day) which the ligands did not exceed as Neoandrographolide showed À0.462 (log mg/kg/day), Bisandrographolide A showed À0.271 log mg/kg/day and Andrographidine C showed 0.263 log mg/kg/day, it does not hamper hERG I & II except Andrographidine C which reduced the chances of QT syndrome that leads to fatal ventricular arrhythmia. Hepatotoxicity caused by drugs is a significant safety issue for drug production and a major cause of drug attrition. Bisandrographolide A demonstrated hepatotoxicity, and the other two molecules had the absence of same.

MD simulation
In our study, Crystal Structure of TNF-alpha APO and complex with selected ligands from docking Andrographidine C and control drug Leflunomide was subjected to molecular dynamic simulation analysis. To understand the stability of the above-mentioned protein-ligand complexes, RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuations), RG (Radius of Gyration), H-Bonds (Hydrogen bonds), and MMPSA calculations were carried out for 100 nanoseconds.

RMSD and RMSF
The stability of both complexes was assessed using the root mean square deviation (RMSD) of back bone atoms to check  whether there was any resemblance or difference. The RMSD is a well-known measure of protein stability and equilibration (Thirumoorthy et al., 2021). To corroborate our docking poses and analyze the average behaviour of the entire protein during MD simulations, we calculated the RMSD of backbone Catoms of all the selected complexes. A lower RMSD value (RMSD $0.2 -$0.3 nm) indicates more stable complex structures, and vice versa (Pandey et al., 2018). The APO, the complex containing Andrographis paniculata molecule i.e. Andrographidine C and the control drug Leflunomide had variations of less than 0.38 nm as depicted in Figure 7. Apo initially had deviation from 0.20 nm to 0.34 nm at 10 ns simulation. Later, from 20 to 40 ns it showed average of 0.25 nm deviation. Furthermore, the complex had average of 0.30 nm. Andrographidine C showed much deviation in beginning but after 50 ns, it got stable at the average of 0.28 nm. Whereas, deviation of leflunomide was stable initially at an average of 0.3 nm till 50 ns. Thereafter, it got increased and went above 0.34 nm.
Whereas Figure 8 illustrates the stability of the ligandreceptor combination. To evaluate the system's residue-byresidue fluctuation, the RMSF (root mean square fluctuation) of C atoms is employed. Because of the restricted motions during simulations, residues with a low RMSF value are more stable. By identifying the protein area that fluctuates during the simulation, RMSF (root means square fluctuation) leads to a better understanding of how ligand binding impacts protein flexibility (Prasanth et al., 2021). When compared to Leflunomide, Andrographidine C shows less variation in majority of locations. Andrographidine C has a lower potential RMSD energy, which allows it to bind to TNF-alpha with greater stability. Andrographidine C variation was decreased in all peaks except one at 100-120 amino acid residues.

Rg
The radius of gyration (Rg) was utilized to measure any changes in the amount of compaction in the Apo protein, Andrographidine C, and control drug Leflunomide complexes during the simulation and to quantify the overall dimensions of the structures. The Rg is the mass-weighted root mean square distance between a collection of atoms and their common centre of mass (Kumar et al., 2014;Maiti et al., 2010). A stable Rg value indicates a securely folded structure, whereas a movement in Rg value throughout the course of the simulation indicates an unfolded structure. When the Rg value is highest, loose packing of amino acids is indicated, and vice versa when the Rg value is lowest (Liao et al., 2014;Lobanov et al., 2008). From Figure 9, stability is observed in the complex as the fluctuation is not major. Andrographidine C initially drops from 2.06 nm to 2 nm and then shows stable deviation whereas Leflunomide shows fluctuations in between range of 2 nm to 2.05 nm. Rg helps identify any structural changes inside a structure and therefore by the figure it is seen that Andrographidine C is significantly more stable than Leflunomide.

H-bond analysis
Hydrogen bonds and their respective strength in a water environment aid protein-ligand interaction . Hydrogen bonds are one of the most significant molecular interactions in biological systems. Hydrogen bonds serve as the cornerstone for molecular identification and selectivity by providing directionality and explicitness to molecular interactions. Changes in secondary structures, which were regulated by hydrogen bonds, guided interactions between proteins and ligands. According to MD simulations, several conformations of a protein may be observed in real biological situations. A protein's form is designed to interact with the ligand in a certain way (Bhardwaj et al., 2021). Figure 10 shows the total number of hydrogen bonds formed for selected complexes during the course of MD simulations. Andrographolide C has a majority of 6 to 7 Hydrogen Bonds, with the maximum being 8 in a few locations. Leflunomide, on the other hand, has the largest number of hydrogen bonds (7).

MM-PBSA
The MmPbSaStat.py python script from the g mmpbsa package was used to compute the average free binding energy of the protein-ligand combination. Table 5 displays van der Waal energy,     electrostatic energy, Polar solvation energy, SASA energy, SAV energy, and WCA energy. By putting all of the energy components together, the ligand's binding energy with the protein is calculated. The binding energies of Andrographidine C and the control drug Leflunomide were À133.927 þ/À 13.866 kJ/mol and À126.106 þ/À 8.671 kJ/mol, respectively. The influence of vander-Walls interactions on total interaction energy is higher than that of electrostatic interactions.

Conclusion
In the present study, we uncovered the finest interacting flavonoids to inhibit the TNF-a inflammation protein. The possible inhibitory potential of the phytochemicals from Andrographis paniculata was analyzed via molecular docking. We selected 20 phytochemicals, which have disclosed anti-inflammatory activity against several inflammatory proteins, were docked against TNF-a. The docking score for twenty phytochemicals was in a range from À8.6 to À7.1 kcal/mol. The final selection was made based on the hydrophilic and hydrophobic amino acid interactions, docking, ADME properties, Toxicity and Molecular dynamics simulations. The results suggested that the compounds Neoandrographolide, Bisandrographolide A and Andrographidine C can be further explored and Andrographidine C could be used as progenitor molecule to institute new druggable compounds targeting destability of trimeric TNF-a against inflammation in Rheumatoid arthritis.

Acknowledgments
Authors are grateful to CHRIST (Deemed to be University) for all the facilities and support provided in this research. We are also thankful to Dr. Usha Talambedu, Maharani Lakshmi Ammanni College for Women for helping us in Molecular Dynamic simulation studies.

Disclosure statement
On behalf of all authors, the corresponding author states that there is no conflict of interest.

Funding
The author(s) reported there is no funding associated with the work featured in this article.