In silico analysis of Triphala-derived polyphenols as inhibitors of TIR–TIR homodimerization in the inflammatory pathway

TL;DR

This study explores how polyphenols from Triphala may inhibit inflammation by interacting with a key protein domain in the NF-κB pathway.

Contribution

The study provides computational evidence that Triphala polyphenols may inhibit MyD88 TIR domain homodimerization, offering a novel anti-inflammatory mechanism.

Findings

Punicalagin and chebulagic acid showed strong binding and stability in the TIR domain of MyD88.

MD simulations confirmed that these polyphenols enhance protein stability and reduce inflammatory signaling.

The findings suggest Triphala polyphenols could serve as safer anti-inflammatory agents.

Abstract

Downstream signaling of the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) pathway is mediated by the adaptor protein myeloid differentiation primary response gene 88 (MyD88). The TIR domain present in MyD88 plays a pivotal role in regulating the expression of pro-inflammatory cytokines. Although synthetic drugs, including M20 and TJ-M2010-5, have been studied to mitigate the overexpression of MyD88, their prolonged usage is known to cause adverse side effects, highlighting the need for a safer, risk-free alternative. An Ayurvedic formulation named Triphala, which is rich in polyphenols and traditionally used to treat various ailments, was selected for this investigation. Although polyphenols are gaining attention as anti-inflammatory agents, their precise mode of action remains insufficiently understood. Previous studies have explored the anti-inflammatory properties of Triphala in a broad spectrum, but this study notably focuses on the interactions of Triphala-derived polyphenols with the TIR domain of the MyD88 adaptor protein in the NF-κB signaling pathway. This study employs computational docking and a molecular dynamics (MD) simulation to study the interaction and stability of the polyphenols with the target protein. The polyphenols were virtually docked to the TIR domain of MyD88 using AutoDock tools 1.5.7. Among them, the top three protein–polyphenol complexes with the highest binding affinities were selected and subjected to MD simulation for 200 ns to evaluate their interaction properties in detail. The findings of the MD simulation corroborated the docking results, showing that two complexes (protein–punicalagin and protein–chebulagic acid) demonstrated better interaction patterns. The MD trajectory revealed that polyphenol binding enhanced the stability of the target protein, as indicated by lower root-mean-square deviation (RMSD) (∼0.25 nm), solvent accessible surface area (SASA) (∼96.848–100.666 nm2), and stabilized radius of gyration (Rg) (∼1.50–1.53 nm) values for punicalagin and chebulagic acid complexes compared to the reference complex. Our findings have supported the hypothesis that Triphala polyphenols may interact with the TIR domain of MyD88, thereby inhibiting the production of inflammatory cytokines. This study provides a combination of computational validation of specific molecular targets and mechanistic insights into the anti-inflammatory potential of Triphala-derived polyphenols.

Figures

Structure of Triphala ligands with the top three binding affinities and the reference compound. (A) Punicalagin, (B) chebulagic acid, (C) sennoside C, and (D) TJ-M2010-5.

Graphical representation of the binding affinities of Triphala-derived ligands against the target protein. The graph illustrates the docking scores of various ligands, wherein the maximum negative score represents the highest binding affinity for the target protein.

The top three protein–ligand complexes with hydrogen bond interactions (studied using LigPlot + v.2.2). (A) Interaction of punicalagin with the target protein. (B) Interaction of chebulagic acid with the target protein. (C) Interaction of sennoside C with the target protein. (D) Interaction of TJ-M2010-5 with the target protein.

Graphical plots for the top three complexes, the reference complex, and the target protein alone. (A) RMSD shows the interaction, stability, and conformational changes of the molecular structures of the different complexes (p-value = 0.0227). (B) RMSF shows the positional fluctuation of individual residues of the protein in different complexes (p-value = 0.0118). (C) SASA shows the degree of the protein surface exposed to the solvent present in the medium (p-value = 0.0130). (D) Radius of gyration (Rg) measures the compactness, folding, and unfolding of the protein structure in different complexes (p-value = 0.0465). (E) H-bond shows the number of H-bonds formed in the complexes during dynamic conditions (p-value = 0.0118). All of the above-mentioned p-values are <0.05, indicating statistical significance. Note: one-way ANOVA with a Brown–Forsythe test was used for comparing each parameter individually in duplicate for significance.

MM-PBSA analysis of the 200°ns trajectory to estimate the binding free energies of the top three protein–ligand complexes, along with the reference complex.

Tables

TABLE 1. Binding affinity or docking scores of different Triphala polyphenols and the reference compound with the TIR domain of the MyD88 protein.

S.No	Ligand	Docking score/binding affinity (kcal/mol)
1	Arjunolic acid	−7.1
2	Ascorbic acid	−4.9
3	Betasitosterol	−6.5
4	Brevifolin	−5.1
5	Chebulagic acid	−8.5
6	Chebulic acid	−6.8
7	Chebulinic acid	−7.8
8	Corilagin	−7.2
9	Ellagic acid	−7.1
10	Emblicanin A	−7.8
11	Emblicanin B	−8
12	Epigallocatechin	−6.7
13	Gallocatechin	−6.8
14	Kaempferol	−7.3
15	Luteolin	−7.9
16	Maslinic acid	−7.4
17	Melissic acid	−4.4
18	Myristic acid	−4.5
19	Palmitic acid	−3.9
20	Pentagalloylglucose	−8.1
21	Phyllaemblicin A	−7.6
22	Phyllaemblicin B	−7.2
23	Phyllaemblicin C	−7.3
24	Phyllanthin	−5.1
25	Progallin A	−6
26	Punicafolin	−7.6
27	Punicalagin	−9.1
28	Quercetin	−7.6
29	Quinic acid	−5.1
30	Sennoside C	−8.4
31	Sennoside E	−8.2
32	Syringic acid	−5.1
33	Tannic acid	−6.3
34	TJ-M2010-5	−7.4

TABLE 2. Formation of hydrogen bonds between the ligands and the target protein with their interacting residues.

Ligands	No. of H bonds formed by the ligands with the protein	Interacting amino acid residues in the protein
Punicalagin	7	Arg133, Asn123, Lys101, Phe109, and Glu108
Chebulagic acid	7	Asp40, Gln81, Ser75, Cys13, Glu77, and Lys76
Sennoside C	11	Trp50, Glu77, Ser15, Gln81, Ser39, Val49, Val43, Leu44, and Thr47
TJ-M2010-5	0	Nil

TABLE 3. MM-PBSA of the top three complexes along with the reference complex.

Energy	Complex
	Punicalagin	Chebulagic acid	Sennoside C	TJ-M2010-5
Complex (ΔGcomplex)	−1974.13	−2,234.98	−2092.39	−2,374.91
Target protein (ΔGprotein)	−2,320.87	−2,353.13	−2,316.08	−2,334.55
Ligands (ΔGligand)	369.52	145.02	23.51	−36.12
Total (ΔGbind)	−22.78	−26.88	−15.83	−4.24

TABLE 4. Drug likeness and pharmacokinetic properties of the top three Triphala-derived polyphenols along with the reference molecule.

Drug likeness and pharmacokinetic properties	Ligands
	Punicalagin	Chebulagic acid	Sennoside C	TJ-M2010-5
Lipinski’s rule of five	3 Violations	3 Violations	3 Violations	0 Violations
Bioavailability score	0.17	0.11	0.11	0.55
Water solubility (log mol/L)	−2.892	−2.892	−2.892	−3.812
Intestinal absorption (%)	100	93.835	0	91.865
BBB permeability (log BB)*	−4.981	−3.889	−2.44	0.287
CYP450 substrate	No	No	No	Yes (CYP2D6 substrate, CYP3A4 substrate)
CYP450 inhibitors	No	No	No	Yes (CYP2D6 inhibitor, CYP3A4 inhibitor)
Total clearance (log mL/min/kg)	−0.246	−2,632	−1.06	0.813
Renal OCT2 substrate	No	No	No	Yes
AMES toxicity	No	No	No	No
Max. Tolerated dose (log mg/kg/day)	0.438	0.438	0.458	0.505
Oral rat acute toxicity (LD50) (mol/kg)	2.482	2.482	2.482	2.66
Oral rat chronic toxicity (LOAEL) (log mg/kg_bw/day)	12.221	10.72	7.099	1.112
Hepatotoxicity	No	No	No	Yes

Full text

Introduction

Myeloid differentiation primary response gene 88 (MyD88) is a cytosolic critical adaptor protein that plays a crucial role in the downstream signaling of the innate immune response (Di Padova et al., 2018). Although there is a vast array of mechanisms and proteins that contribute to the activation of the innate immune response, MyD88 acts as a central mediator in triggering the process (Olson et al., 2015). It facilitates the downstream signaling of the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) pathway and releases pro-inflammatory cytokines. Imbalance of MyD88 triggers a wide spectrum of inflammatory diseases, highlighting its significance in the homeostasis of innate immunity. MyD88 is predominantly activated by Toll-like receptors (TLRs) (except TLR3). TLRs are the major class of pathogen recognition receptors (PRRs) that can recognize the pathogen-associated molecular patterns (PAMPs) present in the microbe (Wang et al., 2014). PAMPs such as lipopolysaccharides (LPS) can bind to the leucine-rich repeat (LRR) motifs of TLRs. This interaction promotes the dimerization of the Toll/interleukin-1 receptor (TIR) domains, subsequently recruiting the TIR domain-containing adaptor proteins, including MyD88 (Clabbers et al., 2021; Saikh, 2021). MyD88 consists of three domains, namely, the N-terminal death domain (DD), the C-terminal Toll/interleukin-1 receptor (TIR) domain, and the intermediate domain (INT). The TIR domain (155–296 aa) facilitates the interaction of MyD88 with the TIR domain of TLR or IL-1R receptors to form a higher-order complex, thereby recruiting additional TIR domains in the vicinity and propagating the downstream signaling cascade. This strongly suggests that TIR–TIR homodimerization may be a viable drug target to inhibit inflammation. Previous studies have focused on developing mimetic small-molecule inhibitors that target the BB-loop of the TIR domain to elucidate the inhibitory mechanism of MyD88 homodimerization (Loiarro et al., 2005; Saikh, 2021).

Several studies have demonstrated that synthetic analogs such as ST2825 effectively inhibit the dimerization of MyD88 by binding to the BB-loop within the TIR domain, thereby disrupting NF-κB-mediated signal transduction in pro-inflammatory cytokine production (Loiarro et al., 2007). Furthermore, ST2825 has been studied to reduce oxidative stress by attenuating the production of reactive oxygen species (ROS) and suppressing NLRP3 inflammasome activation, ultimately inhibiting caspase-1 inflammatory pathways (Zhang et al., 2022). Additionally, TJ-M2010-5 is a mimetic drug designed to interact with the TIR domain of MyD88 to inhibit homodimerization. It binds to the αE, βD, βC, αA, DD, and EE loops but not to the BB-loop of the TIR domain (Xie et al., 2016). However, studies show that this drug has some key limitations and is more effective when administered prophylactically rather than therapeutically (Zou et al., 2023). M20 is another small-molecule inhibitor designed to reduce MyD88 homodimerization by binding to the hydrophobic (αC′–βD–αD) pocket of the TIR domain (Song et al., 2021). These synthetic drugs need to be investigated further for their limitations, mainly regarding their toxicity associated with prolonged usage. In contrast, ancient medicinal systems such as Siddha and Ayurveda provide alternative therapeutic approaches that play a rescuing role in such cases. These natural medicines have potentially lower toxicity profiles and are rich in bioactive compounds that exhibit strong anti-inflammatory properties. Triphala is a well-known ancient polyherbal formulation made up of three fruits, namely, Terminalia chebula, Terminalia bellerica, and Phyllanthus emblica, in a ratio of 1:1:1. This formulation is rich in bioactive compounds such as polyphenols, flavonoids, and tannins. They have been found to play a positive role in the treatment of various inflammatory diseases such as arthritis, cancer, diabetes, and gastrointestinal tract disorders. (Shanmuganathan and Angayarkanni, 2018). Studies have demonstrated that Triphala can hamper the production of pro-inflammatory mediators such as interleukin-6 (IL-6), interleukin-1 beta (IL-1β), and tumor necrosis factor-alpha (TNF-α) by inhibiting the NF-κB signaling pathway (Kalaiselvan and Rasool, 2016). All the above-mentioned studies have explored the general immunomodulatory effects of Triphala. However, no in silico or experimental evidence is currently available for Triphala-derived polyphenol-mediated inhibition of MyD88. Our study provides insights into the direct molecular interactions between Triphala polyphenols and the MyD88 adaptor protein, unveiling the underlying mechanism behind the immunomodulatory properties of these polyphenols through computational screening approaches.

Materials and methods

Database and software

The databases and software applications used in this study include RCSB PDB (Protein Data Bank) (https://www.rcsb.org), PubChem (https://pubchem.ncbi.nlm.nih.gov/), ChemSpider (http://www.chemspider.com), Open Babel version 3.1.1 (https://github.com/openbabel/openbabel/releases), pkCSM webserver (https://biosig.lab.uq.edu.au/pkcsm/prediction), PyMOL 3.0.3 (https://www.pymol.org), AutoDock tools 1.5.7 (https://ccsb.scripps.edu/mgltools/downloads/), Swiss-PDB Viewer version 4.1.0 (https://spdbv.unil.ch/download/binaries/SPDBV_4.10_PC.zip), AutoDock Vina version 1.1.2 (https://vina.scripps.edu/wp-content/uploads/sites/55/2020/12/autodock_vina_1_1_2_linux_x86.tgz), LigPlot + v.2.2 (https://www.ebi.ac.uk/thornton-srv/software/LigPlus/), and GROMACS version 2023.1 (https://zenodo.org/records/7852175/files/gromacs-2023.1.tar.gz?download=1).

Preparation of the protein

The crystal structure of the TIR domain of the adaptor protein MyD88 (PDB ID: 4DOM) was retrieved from the RCSB PDB database, which contains 151 amino acids (Clabbers et al., 2021). The structure was determined by X-ray diffraction to a resolution of 2.30 Å. Heteroatoms, water molecules, and undesired ligand molecules were removed using PyMOL. Missing atoms were fixed using Swiss-PDB Viewer version 4.1.0. Using AutoDock tools 1.5.7, polar hydrogens and Kollman charges were added to stabilize the protein structure and were saved in the PDBQT file format for molecular docking (Hu et al., 2024).

Ligand preparation

A total of 33 Triphala-derived bioactive compounds with strong interaction properties toward the target protein were selected based on a literature survey. For instance, all the selected ligands have been frequently explored for their strong anti-inflammatory, anticancer, and antioxidant properties. The 3D structures of these ligands were retrieved from the PubChem and ChemSpider databases (Tarasiuk et al., 2018; Prasad and Srivastava, 2020; Rudrapal et al., 2022; Wang et al., 2023). Energy minimization of all the ligands was performed by applying Merck Molecular Force Field 94 (MMFF94) to obtain stable conformations, as MMFF94 was specifically developed to provide accurate intermolecular geometries for small organic molecules. All the ligands were then converted to the PDBQT format for molecular docking (Murugesan et al., 2021; Olaokun et al., 2022).

Molecular docking

The processed ligands were subjected to docking against the TIR domain of MyD88 using AutoDock Vina version 1.1.2. In this study, the ligands that were energy minimized were given as the input file, and the prepared protein was given as the target file. To perform high-throughput molecular docking, a Perl script was used containing multiple ligand files to execute AutoDock Vina (Eberhardt et al., 2021). The blind docking covers the whole protein, and the center of the grid box was precisely positioned at the coordinates X: 17.67, Y: 1.496, and Z: 12.318 with a size range of X: 58 Å, Y: 54 Å, and Z: 46 Å. This configuration facilitates the sampling of the entire protein surface by the ligands to explore all the potential binding pockets. A maximum of 10 docking conformers were set, from which the most favorable binding energies were selected, and the root-mean-square deviation (RMSD) values were saved (Aliye et al., 2021; Eberhardt et al., 2021). Complexes of the saved confirmations of ligands and target proteins were written using PyMOL 3.0.3. Non-covalent interactions of the complexes were studied using LigPlot+ (Laskowski and Swindells, 2011). Furthermore, an MD simulation along with calculations of binding free energies was carried out for these complexes to understand the dynamic properties involved in the interaction and stability of the protein–ligand complexes (Rudrapal et al., 2022).

Molecular dynamics simulation

An MD simulation was performed for the complexes with the top three binding affinities, the reference drug (TJ-M2010-5), and the MyD88 protein alone using GROningen MAchine for Chemical Simulations (GROMACS) version 2023.1 (Abraham et al., 2015; Rudrapal et al., 2022). The CHARMM (Chemistry at HARvard Macromolecular Mechanics) force field was employed for the MD simulation of all the complexes due to its robust parameterization and reliability in modeling biomolecular systems, including small molecules and phytochemicals. The topology file for the target protein was prepared using GROMACS, and the target ligands were prepared using CGenFF (Vanommeslaeghe et al., 2010). The periodic boundary was set using a rectangular box, and the water box with solute molecules surrounded by the solvent was extended to 10 Å. Bond constraints were handled using the LINCS algorithm. All the complexes were solvated using the TIP3P water model, followed by sodium and chloride ions being added for neutralization (Ghosh et al., 2021; Murugesan et al., 2021). Energy minimization was carried out for all the complexes using the steepest descent algorithm. Equilibration was carried out under isothermal–isochoric (NVT) and isothermal–isobaric (NPT) conditions for 0.3 ns each (Rudrapal et al., 2022). The systems were then subjected to an MD simulation under a pressure of 1 atm and a temperature of 300 K for 200 ns (duplicates were run to ensure reproducibility). Here, the leapfrog algorithm was used at 2 fs for integrating the motion equation of atoms. Finally, the average distance between the atoms of the overlying protein–ligand complex was studied by RMSD. Deviations in the positions of the atoms were studied by root-mean-square fluctuation (RMSF). The overall size and compactness were studied using the radius of gyration (Rg), and the accessibility of the surface area to the solvent molecules was quantified by solvent accessible surface area (SASA). All the above four parameters corroborated the stability of the target protein in the presence of ligands, suggesting the sustainability of the interaction. The non-covalent associations (H-bonds) were calculated to evaluate the stability and specificity of the molecular interactions using GROMACS scripts (Ghosh et al., 2021; Rudrapal et al., 2022).

Molecular mechanics Poisson–Boltzmann surface area

Binding free energy is the energy that tends to be released when a protein and a ligand form a complex, and it can be calculated using the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method. MM-PBSA helps in understanding the thermodynamics and kinetics involved in the processes of molecular binding. Dielectric constants are the key parameters to calculate MM-PBSA, which enables the accurate estimation of the electrostatic contribution to the reaction medium. An internal dielectric constant (indi) of 4.0 was adopted to approximate the electric polarizability of the solute (protein), whereas an external dielectric constant (exdi) of 80.0 was set to simulate the environment of bulk water at room temperature (Hou et al., 2011). Snapshots were recorded at an interval of 0.1 ns (100 picoseconds) throughout the 200 ns simulation. The frames of the trajectory used for the MM-PBSA calculation were extracted using the parameter nstxtcout = 50,000, corresponding to one frame every 0.1 ns. A total of 2,000 frames were taken to cover the whole trajectory of 200 ns rather than focusing on a segment of the trajectory. The parameters, such as changes in bond energies, angles, electrostatics, and van der Waals and solvation-free energy, which is the sum of polar and nonpolar solvation energies, were calculated. The overall binding free energy was calculated using the following formula (Kumari et al., 2014; Murugesan et al., 2021; Rudrapal et al., 2022):

$$[eqn]$$

where ΔG_bind_ is the binding free energy between the ligand and protein, ΔG_complex_ is the binding free energy released by the protein–ligand complex, ΔG_protein_ is the binding free energy of the unbound protein, and ΔG_ligand_ is the binding free energy of the unbound ligand.

Drug-likeness and pharmacokinetic profile

Drug-likeness is a qualitative property that indicates whether a drug molecule can act as a potential drug. The rule of five (Ro5) or Lipinski’s rule of five provides a set of guidelines to evaluate the drug-likeness of oral drugs (Lipinski, 2004). These rules provide criteria for certain molecular properties that play a major role in pharmacokinetic profiles, such as absorption, distribution, metabolism, and excretion (ADME). In this study, SwissADME and pkCSM were used to study drug-likeness, physicochemical properties, and pharmacokinetics (Pires et al., 2015; Daina et al., 2017; Ghosh et al., 2021).

Results

Molecular docking

A list of 33 polyphenols derived from the Ayurvedic formulation Triphala were taken and computationally docked with the crystal structure of the TIR domain of the MyD88 protein to evaluate their interaction and binding affinities (Tarasiuk et al., 2018; Prasad and Srivastava, 2020; Wang et al., 2023). TJ-M2010-5, a synthetic anti-inflammatory drug, was taken as a reference molecule (Xie et al., 2016). Binding affinities were assessed based on the docking scores, expressed in kcal/mol (Song et al., 2021; Hu et al., 2024). All the polyphenols exhibited favorable interaction patterns, and few polyphenols, such as chebulagic acid, chebulinic acid, emblicanin B, emblicanin A, luteolin, malinic acid, maslinic acid, pentagalloylglucose, punicalagin, punicafolin, quercetin, sennoside C, and sennoside E, exhibited good binding affinities toward the target protein compared to the reference molecule. In particular, punicalagin (−9.1 kcal/mol), sennoside C (−8.4 kcal/mol), and chebulagic acid (−8.5 kcal/mol) showed the highest binding affinities compared to the reference drug TJ-M2010-5 (−7.4 kcal/mol), suggesting that the polyphenols have better interaction properties with the target protein than the reference drug (Table 1) (Figures 1, 2). However, the practical application of polyphenols as inhibitors of MyD88 depends on several additional factors such as molecular interaction, stability, and binding free energy, which were also elaborately investigated in the present study. The non-covalent interactions between the protein–ligand complexes were examined in detail using LigPlot + V.2.2. In this study, punicalagin was observed to interact with Arg133, Asn123, Lys101, Phe109, and Glu108 amino acid residues present in the TIR domain of MyD88 through seven hydrogen bonds (reversible non-covalent bonds). Similarly, chebulagic acid was observed to have established seven hydrogen bonds with Asp40, Gln81, Ser75, Cys13, Glu77, and Lys76 amino acid residues of the target protein. Notably, sennoside C established 11 hydrogen bonds with Trp50, Glu77, Ser15, Gln81, Ser39, Val49, Val43, Leu44, and Thr47 acid residues of the target protein. In contrast, the reference drug TJ-M2010-5 did not form any hydrogen bonds with the target protein, which shows that the interaction profile of the reference drug was weaker than that of Triphala-derived polyphenols (Table 2) (Figure 3).

Structure of Triphala ligands with the top three binding affinities and the reference compound. (A) Punicalagin, (B) chebulagic acid, (C) sennoside C, and (D) TJ-M2010-5.

Graphical representation of the binding affinities of Triphala-derived ligands against the target protein. The graph illustrates the docking scores of various ligands, wherein the maximum negative score represents the highest binding affinity for the target protein.

The top three protein–ligand complexes with hydrogen bond interactions (studied using LigPlot + v.2.2). (A) Interaction of punicalagin with the target protein. (B) Interaction of chebulagic acid with the target protein. (C) Interaction of sennoside C with the target protein. (D) Interaction of TJ-M2010-5 with the target protein.

Molecular dynamics simulation

The top three complexes with higher binding affinities, the reference complex, and the target protein were studied for their stability, interaction, and confirmation in detail by performing an MD simulation for 200 ns (Ghosh et al., 2021; Aljarba et al., 2022). The stability of the complexes was assessed by studying key parameters, including RMSD, RMSF, SASA, Rg, and the hydrogen bonds formed. Graphical representations for these parameters were generated using the Xmgrace 2D plotter (Ghosh et al., 2021; Rudrapal et al., 2022; Vieira et al., 2023). RMSD was employed to study the conformational stability and structural integrity of the molecules present in the protein–ligand complexes (Ghosh et al., 2021; Rudrapal et al., 2022). The RMSD of the target protein was gradually increased from 0 to 33.5 ns and reached a peak of 0.28 nm, followed by a sudden decrease. Subsequently, the RMSD values showed continuous fluctuations and reached up to 0.27 nm at 58.4 nm, before a gradual decrease until 92.8 ns. Following that, the RMSD stabilized with minor fluctuations; the protein backbone was normalized after 149.7 ns and remained stable throughout the simulation with minor fluctuations. The RMSD profile of the protein–TJ-M2010-5 complex was more variable in the beginning, reaching up to 0.25 nm at 13.4 ns, followed by minor fluctuations throughout the MD simulation for 200 ns, indicating the structural instability of the complex. The protein–punicalagin complex first displayed significant fluctuations in RMSD and reached a maximum value of 0.25 nm at 48.2 ns and 145.8 ns. After that, the complex stabilized with some minor deviations. The protein–chebulagic acid complex exhibited oscillations in RMSD from 0 to 107.7 ns and reached a maximum value of 0.25 nm at 69.9 nm. Subsequently, the values stabilized with slight deviations until 145.8 ns and remained stable throughout the 200 ns. The RMSD exhibited by the protein–sennoside C complex was the highest of all the other complexes and the target protein, which reached a maximum value of 0.29 nm at 74.8 ns. The trajectory showed continuous oscillations until 146.8 ns, followed by the complex being partially stabilized with minor fluctuations (Figure 4A). Overall, punicalagin and chebulagic acid stabilized the backbone of the target protein compared to other ligands.

Graphical plots for the top three complexes, the reference complex, and the target protein alone. (A) RMSD shows the interaction, stability, and conformational changes of the molecular structures of the different complexes (p-value = 0.0227). (B) RMSF shows the positional fluctuation of individual residues of the protein in different complexes (p-value = 0.0118). (C) SASA shows the degree of the protein surface exposed to the solvent present in the medium (p-value = 0.0130). (D) Radius of gyration (Rg) measures the compactness, folding, and unfolding of the protein structure in different complexes (p-value = 0.0465). (E) H-bond shows the number of H-bonds formed in the complexes during dynamic conditions (p-value = 0.0118). All of the above-mentioned p-values are <0.05, indicating statistical significance. Note: one-way ANOVA with a Brown–Forsythe test was used for comparing each parameter individually in duplicate for significance.

RMSF analysis helps identify regions of increased flexibility and stability within proteins. The RMSF profile of the target protein varied significantly between 2–5, 40–50, 88–96, 105–117, and 123–129 residues of amino acids. For the protein–TJ-M2010-5 complex, major fluctuations in the RMSF were found between the 2–5, 43–50, 86–96, 105–117, 123–129, and 141–142 residues of amino acids. The complex showed more fluctuations in its residues than the target protein. For the protein–punicalagin complex, major fluctuations were observed between the 2–5, 50, 89–95, 106–107, 112–115, 123–125, and 129 amino acid residues. The complex exhibited fewer fluctuations compared to the target protein, showing the stabilizing effect of punicalagin when bound to the target protein. Similarly, the protein–chebulagic complex displayed fluctuations between 2–5, 41–50, 90–96, and 123–127 amino acid residues, which have significantly lower RMSF than the target protein. The protein–sennoside C complex exhibited fluctuations between the 2–5, 38–50, 89–95, 114–115, 123–126, and 142 amino acid residues. This complex showed lower fluctuations than the target protein, although the fluctuations were slightly higher than those of punicalagin and chebulagic acid complexes (Figure 4B).

SASA was analyzed to gain insights into the structural and functional properties of biomolecules by quantifying the protein surface exposed to the solvent present in the environment. Effective ligand binding was expected to minimize the solvent accessible surface area of the protein. The target protein exhibited a maximum SASA value of 95.450 nm^2^. In contrast, the reference complex showed a maximum SASA value of 103.325 nm^2^, which is higher than that of the target protein, indicating that ligand binding was disturbed by MD simulation. The maximum SASA values recorded for the protein–punicalagin, protein–chebulagic acid, and protein–sennoside C complexes were 100.666, 96.848, and 101.003 nm^2^, respectively. Subsequently, the values were reduced during a further simulation for 200 ns, suggesting potential protein stabilization by ligands over the course of the MD simulation (Figure 4C). Overall, punicalagin and chebulagic acid complexes exhibited lower SASA values than the reference drug. On the other hand, sennoside C showed elevated SASA values, similar to the reference drug. This data trend shows a strong correlation with the docking results, in which punicalagin and chebulagic acid exhibited the lowest binding energies, which ultimately implies stable binding with the target protein. Rg shows the compactness of the molecular system, thereby providing insights into the stability, folding, and unfolding conformational changes in the protein. The Rg of the target protein showed minor fluctuations initially and gradually increased, reaching a maximum of 1.52 nm on 98.4 ns. Subsequently, the system was stabilized after 148.4 ns and remained the same throughout the simulation. Conversely, the Rg value of the reference complex displayed continuous fluctuations throughout the simulation and reached a maximum of 1.52 nm at 98.4 ns. This suggests that ligand binding negatively influences the rigidity of the protein structure. Initially, the Rg values of the punicalagin, chebulagic acid, and sennoside C complexes oscillated and reached maximums of 1.53, 1.50, and 1.51 nm, respectively. Following that, the Rg values of punicalagin and chebulagic acid complexes were stabilized after 147.7 ns and 148.8 ns, respectively. Notably, both punicalagin and chebulagic acid complexes exhibited better structural rigidity than the other complexes (Figure 4D). Furthermore, the stability of the complexes was assessed by studying the number of intermolecular hydrogen bonds, which helps maintain the complex intact. In this study, the reference complex did not form any hydrogen bonds throughout the simulation, indicating an unstable binding interaction. The H-bonds formed by the punicalagin, chebulagic acid, and sennoside C complexes varied between 1–7, 1–9, and 1–8, respectively, throughout the simulation (Figure 4E). Notably, the chebulagic acid and sennoside C complexes displayed an increased and decreased number of H-bonds in post-simulation compared with the binding energy from the docking results. An increase in bond number represents the additional conformational states that were explored by the ligand, thereby enhancing its interaction and stability under dynamic conditions, whereas a decline in bond number suggests a weaker binding of the ligand. The statistical significance tests were carried out individually for all the parameters using one-way ANOVA with the Brown–Forsythe test, and all the parameters were statistically significant (p-values are <0.05).

Molecular mechanics Poisson–Boltzmann surface area

Binding free energy (ΔG_bind_) was computed using the sum of SASA, van der Waals, electrostatic energy, and polar and non-polar solvation energies. MM-PBSA helps in comparing the binding affinities of different proteins and ligands. A more negative ΔG_bind_ indicates a stronger binding affinity between the ligand and the target protein. The results of MM-PBSA provided insights into the energy contributed to the binding of the target protein and ligands. The average energy exhibited by the protein–punicalagin complex was −1,974.13 kcal/mol, whereas the protein and ligand alone displayed −2,320.87 kcal/mol and 369.52 kcal/mol, respectively, leading to a total ΔG_bind_ of −22.78 kcal/mol. Similarly, the protein–chebulagic acid complex exhibited an average energy of −2,234.98 kcal/mol, with the unbound protein and ligand displaying −2,353.13 kcal/mol and 145.02 kcal/mol, respectively, yielding the total ΔG_bind_ of −26.88 kcal/mol. The protein–sennoside C complex showed an average energy of −2,092.39 kcal/mol, the protein alone exhibited −2,316.08 kcal/mol, the ligand alone exhibited 23.51 kcal/mol, and the total ΔG_bind_ was −15.83 kcal/mol. For comparison, the reference complex was also analyzed for its ΔG_bind_, which was approximately −4.24 kcal/mol. The energy displayed by the protein–TJ-M2010-5 complex was −2,374.9 kcal/mol, the protein alone exhibited −2,334.55 kcal/mol, and the ligand alone exhibited −36.12 kcal/mol. Overall, punicalagin and chebulagic acid, among the other ligands, demonstrate the most favorable thermodynamic integrity. The MM-PBSA values of the polyphenols were consistent with the SASA and H-bond results, which suggest a reduction in the exposure of the protein surface toward the solvent upon ligand binding and stable intermolecular hydrogen bonding throughout the MD simulation of 200 ns (Table 3) (Figure 5).

MM-PBSA analysis of the 200°ns trajectory to estimate the binding free energies of the top three protein–ligand complexes, along with the reference complex.

Drug likeness and pharmacokinetic profile

SwissADME and pkCSM were utilized to study Lipinski’s rule of five and the pharmacokinetic properties of the ligands (Ghosh et al., 2021). The synthetic reference drug exhibited no violations of Lipinski’s rule of five, whereas all three Triphala-derived ligands showed three violations: 1. molecular weight exceeding 500 Da, 2. nitrogen-oxygen bonds exceeding 10, and 3. nitrogen–hydroxyl bonds surpassing 5. Specifically, the molecular weights of punicalagin, chebulagic acid, and sennoside C were 1,084.72, 954.66, and 848.76 Da, respectively. All three polyphenols exceeded the threshold for nitrogen–oxygen bonds and nitrogen–hydroxyl bonds up to 10 and 5, respectively. The normal bioavailability score range was 0–10, and all the polyphenols showed poor bioavailability (0.11–0.17), while the reference compound exhibited a significantly higher bioavailability score of 0.55 (Martin, 2005; Daina et al., 2017). All the polyphenols displayed a water solubility of −2.892, while the reference compound had a solubility of −3.812, suggesting that the ligands ranged from moderately soluble to soluble (Daina et al., 2017). Further pharmacokinetic evaluation indicated that punicalagin, chebulagic acid, and TJ-M2010-5 demonstrated good intestinal absorption (100%, 93.835%, and 91.865%, respectively), while sennoside C showed a very poor intestinal absorption rate of 0% (Zhao et al., 2002). Even though sennoside C exhibited significant binding affinity toward the TIR domain, this major limitation, the lack of intestinal absorption, reduces its oral efficacy, and it can be overcome by advanced drug delivery platforms such as nanoparticle encapsulation (Sun et al., 2024).

The blood–brain barrier (BBB) membrane permeability revealed that TJ-M2010-5 exhibited a permeability of 0.287 log BB, whereas Triphala polyphenols showed very poor permeability (<-1) compared to the reference compound. Furthermore, none of the ligands met the range of −0.89 to 3.32 log mL/min/kg in total clearance (Zhivkova and Doytchinova, 2013). Cytochrome P450 (CYP450) enzymes play a crucial role in drug metabolism, particularly with CYP2D6/CYP3A4. An ideal drug should function as a CYP substrate for proper metabolic processing. The reference drug functions as a CYP substrate, whereas all the other tested ligands failed to act as a CYP substrate. Additionally, a drug candidate should not inhibit CYP450 enzymes. All the polyphenols did not show any inhibitory properties toward CYP450 enzymes; in contrast, the reference drug demonstrated inhibitory effects over these enzymes. Organic cation transporter 2 (OCT2) is a renal transporter that plays a crucial role in the excretion of drugs. An ideal drug should serve as an OCT2 substrate for proper excretion. However, all three polyphenols failed to exhibit the property of being OCT2 substrates; in contrast, the reference drug acted as an OCT2 substrate. Toxicity assessments showed that the Triphala-derived ligands did not exhibit hepatotoxicity or AMES (mutagenic) toxicity, whereas the reference compound displayed hepatotoxicity. These computational data align with the previously reported in vivo toxicity studies conducted in Sprague–Dawley rats, which suggest that the Triphala extract administration is not associated with hepatotoxicity. This was evidenced by the normal serum levels of hepatic biomarkers, including aspartate aminotransferase (AST), alanine aminotransferase (ALT), and alkaline phosphatase (ALP) levels (Arpornchayanon et al., 2022). The predicted maximum tolerated doses for punicalagin, chebulagic acid, and sennoside C were 0.438, 0.438, and 0.458 log mg/kg/day, respectively, while the reference drug had a higher tolerated dose of 0.505 log mg/kg/day (Pires et al., 2015). Rat oral acute and chronic toxicities were also evaluated with the pkCSM predictive model; the lethal dose LD50 for all the polyphenols was commonly in the acute toxicity range of 2.482 mol/kg, whereas the reference compound had a slightly higher LD50 value at 2.66 mol/kg. The lowest observed adverse effect level (LOAEL) in rat models was predicted for all the Triphala-derived polyphenols and the reference drug to study the oral chronic toxicity. The LOAEL values for punicalagin, chebulagic acid, and sennoside C were 12.221, 10.72, and 7.099 (log mg/kg_bw/day), respectively (Table 4) (Pires et al., 2015). TJ-M2010-5 exhibited a significantly lower LOAEL of 1.112 compared to the polyphenols, suggesting a higher potential for chronic toxicity. Although these polyphenols are broadly recognized for their anti-inflammatory properties, and numerous in vitro and in vivo studies have demonstrated their efficacy in inflammatory models, significant limitations still persist. To overcome these challenges, prodrug strategies can be employed to enhance their ADME profiles. Additionally, advanced drug delivery platforms such as nanocarrier-based formulations offer promising approaches to improve the target delivery and bioavailability of the polyphenols, thereby increasing their therapeutic potential (Rautio et al., 2008; Sun et al., 2024).

Discussion

The first line of defense in the human body is provided by the innate immune response, which is a nonspecific protective mechanism that involves numerous signaling pathways. Among these, NF-κB is one of the most important pathways, which involves the adaptor protein MyD88. MyD88 mediates the downstream signaling of the NF-κB pathway, leading to the subsequent release of pro-inflammatory cytokines. MyD88 comprises three domains, namely, the death domain (DD) at the N-terminal, the Toll/interleukin-1 receptor (TIR) domain at the C-terminal, and an intermediary domain (INT) (Ohnishi et al., 2009; Saikh, 2021). The TIR domain (155–296 aa) facilitates MyD88 to interact with the TIR domain of Toll-like receptors (TLRs) or interleukin-1 receptors (IL-1Rs) to form higher-order complexes by recruiting more TIR domains of MyD88 to the vicinity. DD (54–109 aa) activates interleukin-1R (IL-1R)-associated kinase (IRAK) 1–4, subsequently undergoing various phosphorylation procedures and releasing the transcription factor NF-κB. The INT (110–155 aa) domain forms a link between the TIR domain and DD. Studies show that INT is essential for IRAK phosphorylation, and its absence is associated with the inability of MyD88 in downstream signaling (Clabbers et al., 2021). Since the TIR domain of MyD88 plays an important role in the triggered immune response, many synthetic drugs have been used for targeting this domain, which were found to have harmful side effects. Polyphenols are naturally occurring phytochemicals in all plant-based foods, known for their anti-inflammatory properties and toxin-free nature. A search for a toxin-free Ayurvedic drug with significant anti-inflammatory properties has led to the drug named Triphala, which is made from three fruits. These fruits are rich in polyphenols that can work synergistically to provide health benefits, including anti-inflammatory properties. Hence, the 33 major polyphenols were selected, and their crystal structures were optimized with the MMFF94 force field for docking studies with the TIR domain of MyD88. Our docking studies revealed that the polyphenols had good binding affinity toward the TIR domain of the adaptor protein, which was expressed in binding energy. In particular, punicalagin, chebulagic acid, and sennoside C demonstrated the highest binding affinities for the target protein, with more negative scores. In molecular docking, high negative scores indicate efficient binding; however, the interaction, stability, and conformational changes of these complexes were further confirmed by studying their non-covalent bonds and by subjecting them to an MD simulation for 200 ns (Ghosh et al., 2021; Rudrapal et al., 2022). MD simulation is a computational approach that aids in the real-time simulation of the movements of atoms, providing insights into the kinetic and thermodynamic characteristics of biomolecules (Gelpi et al., 2015; Hollingsworth and Dror, 2018). Key integrity parameters, such as RMSD, RMSF, SASA, Rg, and H-bonds, were evaluated for the stability and dynamicity of atoms present in protein–ligand complexes. Notably, in our study, the punicalagin and chebulagic acid complexes exhibited a stable graph of RMSD, SASA, and Rg, while the protein–TJ-M2010-5 and protein–sennoside C complexes demonstrated more fluctuations in the RMSD, SASA, and Rg, suggesting structural instability compared to other complexes and the target protein. In general, the RMSD is a key metric that helps in assessing the stability of molecular structure, and a higher RMSD value indicates less stability. SASA measures the extent of protein exposure to the solvent. A high SASA value indicates that a large portion of the protein is exposed to the solvent, suggesting that the ligands are not properly bound to the protein (Arooj et al., 2020; Ghosh et al., 2021; Rudrapal et al., 2022). Rg provides insights into the compactness and rigidness of atoms in biomolecules. Higher Rg values indicate the increased flexibility of a protein with weaker interactions between the protein and ligand. RMSF provides insights into the atomic-level modifications in the protein chain, which aids in studying the flexibility and stability of residues present in the protein. RMSF has a direct influence on the affinity and specificity of ligand binding. A protein–ligand complex with lower RMSF denotes the strong interaction between them (Arooj et al., 2020; Ghosh et al., 2021; Rudrapal et al., 2022). Our study showed that the RMSF of the reference drug complex had higher fluctuations than the polyphenol complexes. Overall, our MD trajectories demonstrated that the binding of polyphenols increased the stability of the protein structure more than the reference drug. Notably, the strength of ligand binding was governed by the non-covalent hydrogen bonds. Initially, the proteins may display numerous intramolecular H-bonds in their naïve state to stabilize their secondary and tertiary structures. Upon binding of ligand molecules, some of the intramolecular bonds are degraded and replaced by intermolecular H-bonds (Arooj et al., 2020; Ghosh et al., 2021; Rudrapal et al., 2022). Intermolecular H-bonds were quantified after docking and MD simulation. Notably, the H-bond in the punicalagin complex remained the same after docking and MD simulation, whereas the chebulagic acid complex showed increased H-bonds after MD simulation. Additional H-bond formation during MD simulation suggests that the protein–ligand system has explored multiple conformations under dynamic conditions.

The binding free energy was calculated with MM-PBSA, which provided insights into the energy contributed to the binding of the target protein and ligands (Rifai et al., 2019; Murugesan et al., 2021). In general, the higher the binding free energy, the lower the binding affinity. Our findings revealed that the reference drug has a higher binding free energy than the tested polyphenols. Meanwhile, the chebulagic acid and punicalagin complexes displayed very low binding free energies in MM-PBSA, which correlates with the other MD trajectories. Overall, MD simulation data corroborated the findings of molecular docking, suggesting that the binding of punicalagin and chebulagic acid has enhanced the stability and integrity of the target protein, with significant stabilization observed after 145–150 ns, as evidenced by a peak reduction in RMSD, Rg, and SASA values. In addition, evaluation of the drug-likeness and pharmacokinetic properties was essential for knowing the druggability of a therapeutic compound to avoid drug rejection during the clinical trials (Aljarba et al., 2022). Here, the selected polyphenols exhibited limitations in ADME properties, particularly in absorption and metabolism. Unlike the reference compound, the polyphenols have not exhibited hepatotoxicity or mutagenicity. Incorporation of sustained drug delivery systems can improve the pharmacokinetic performance of these polyphenols. However, direct experimental evidence of these polyphenols interacting with the TIR domain of the MyD88 adaptor protein was not conclusively demonstrated. Studies on in vitro and in vivo inflammatory models concluded that using Triphala extracts can suppress the inflammatory biomarkers, including TNF-α, IL-1β, and IL-6, suggesting the indirect modulation of MyD88-mediated NF-κB expression (Kalaiselvan and Rasool, 2016). Further in vitro studies are needed to elucidate the precise mechanism of polyphenol interactions with MyD88. Overall, our findings provide strong evidence for supporting the hypothesis that Triphala-derived polyphenols modulate the triggered immune response by interacting with the TIR domain of MyD88, thereby inhibiting homodimerization, which may further suppress the downstream immune signaling pathways. This result highlights the therapeutic potential of Triphala-derived polyphenols in the management of inflammatory diseases.

Conclusion

In conclusion, this study employed a computational approach to provide insights into the molecular mechanisms underlying the anti-inflammatory potential of Triphala-derived polyphenols. The findings highlighted the non-toxic nature of polyphenols and their favorable interaction profiles with the TIR domain of the MyD88 adaptor protein. Compounds such as chebulagic acid and punicalagin have shown the maximum interaction with the target protein, which could act as promising scaffolds for drug design. However, to comprehensively evaluate the proposed hypothesis, further in vitro investigations using cell lines such as RAW 264.7 or THP-1 are required to validate the mechanism of NF-κB inhibition and cytokine reduction. Additionally, in vivo inflammatory models could provide mechanistic insights into the pharmacodynamics of these compounds, thereby supporting the formulation of potential therapeutic agents to combat chronic inflammatory conditions.