Original Research•

Unveiling molecular insights: in silico exploration of TLR4 antagonist for management of dry eye syndrome

Abstract

Background Dry eye disease is the most commonplace multifractional ocular complication, which has already affected millions of people in the world. It is identified by the excessive buildup of reactive oxygen species, leading to substantial corneal epithelial cell demise and ocular surface inflammation attributed to TLR4. In this study, we aimed to identify potential compounds to treat of dry eye syndrome by exploring in silico methods.

Methods In this research, molecular docking and dynamics simulation tests were used to examine the effects of selected compounds on TLR4 receptor. Compounds were extracted from different databases and were prepared and docked against TLR4 receptor via Autodock Vina. Celastrol, lumacaftor and nilotinib were selected for further molecular dynamics studies for a deeper understanding of molecular systems consisting of protein and ligands by using the Desmond module of the Schrodinger Suite.

Results The docking results revealed that the compounds are having binding affinity in the range of −5.1 to −8.78 based on the binding affinity and three-dimensional interactions celastrol, lumacaftor and nilotinib were further studied for their activity by molecular dynamics. Among the three compounds, celastrol was the most stable based on molecular dynamics trajectory analysis from 100 ns in the catalytic pockets of 2Z63.pdb.pdb. Root mean square deviation of celastrol/2Z63 was in the range of 1.8–4.8 Å.

Conclusion In particular, Glu376 of TLR4 receptor is crucial for the identification and binding of lipopolysaccharides (LPS), which are part of Gram-negative bacteria’s outer membrane. In our investigation, celastrol binds to Glu376, suggesting that celastrol may prevent the dry eye syndrome by inhibiting LPS’s binding to TLR4.

What is already known on this topic

  • The innate immune system mediates inflammation, which damages the ocular surface by this factor toll-like receptors are upregulated and cause dry eye disease.

What this study adds

  • Celastrol binds with the Glu376 of 2Z63.pdb, which is important for the recognition and binding of lipopolysaccharides to TLR4 in Gram-negative bacteria.

How this study might affect research, practice or policy

  • Through this research work, we found celastrol has good binding affinity and stability towards TLR4 receptor thereby treating dry eye disease

Introduction

Dry eye syndrome is a leading multifractional ocular disease, which affects millions of individuals worldwide.1 Untreated chronic dry eye disease (DED) has a high risk for ocular surface infection, immunoinflammatory and corneal epithelial cell damage, which may result in poor vision.2 Microbial pathogens on the ocular surface are recognised by lipopolysaccharide (LPS) and activate immunological defence, which enables the induction of rapid and extensive macrophages by binding TLR4 on the cell surface, and inflammatory processes are triggers through the activation of synthesis of proinflammatory cytokines and the NF-κB.3 Toll-like receptors are a well-known family of glycoprotein receptors, and all TLRs have been reported to be expressed in the cornea and conjunctiva, although some differences in subcellular localisation remain.4 Corneal inflammation and damage induced by DED can trigger the generation of endogenous TLR4 ligands, activating the immune system. The synergy of biological systems and computational approaches enables the exploration of drug development, facilitating the swift acquisition of structural, chemical and biological data. This integrated approach enhances our comprehension of potential drug targets. This combination allows us to understand potential drug targets. Researchers have used a variety of methods, such as virtual screening, molecular docking and interaction of various immune systems, to develop therapeutic agents.5 In this study, we performed in silico method to analysis the interactions between three compounds, namely nilotinib, lumacaftor and celastrol against the receptor structure, namely the ‘crystal structure of the TV8 hybrid of human TLR4 and hagfish’ (PDB ID: 2Z63.pdb). The primary objective of this study was to gain insights into the binding modes and dynamics of these protein-ligand interactions through a combination of molecular docking and dynamics model. Getting knowledge of the interactions between proteins and small molecules is crucial for the development of novel therapeutic agents and drug design strategies. By employing advanced computational techniques, we aimed to shed light on the structural and dynamic aspects of the protein-ligand complexes and elucidate their functional implications. Studies have shown that topical application of TLR agonists on the corneal epithelium can cause intense facial inflammation.6 Corneal inflammation signal pathway can be stopped by tlr4 antagonist by blocking the binding of ligands to receptor and interfering with intracellular signalling pathways to prevent signal transduction to stimulation with LPS,7 suggesting the TLR antagonists have potential therapeutic effects in modulating corneal inflammation. Therefore, in the current study, we report the successful in silico screening of ligand and structure-based molecular docking and dynamic in the discovery of novel TLR4 antagonists.8

Methodology

Compound selection

1000+ herbal and synthetic compounds are selected based on anti-inflammatory activity from that we shortlisted 30 compounds based on our preference such as toxicity, affinity, docking score, molecular weight, solubility, for ocular formulation. Our screening process commenced with the meticulous selection of a meticulously curated list comprising 30 compounds, which encompassed a gamut of chemical entities, including nandrolone phenpropionate, adapalene, antrafenine, nilotinib, difenoxin, irinotecan, lumacaftor, cabozantinib, pranlukast, ergocalciferol, sparstolonin B, atractylenolide I, oscillatoria, celastrol, loteprednol, varenicline, epigallocatechin gallate, kaempferol, ferulic acid, meclizine, lactoferrin, lifitegrast, xanthohumol, resveratrol, quercetin, vorapaxar, atractylenolide I, resatorvid, ciclosporin and hispidulin.

Preparation of ligands

To pave the way for in-depth investigations, we commenced by procuring the two-dimensional (2D) structure in SDF format for all 30 compounds. Subsequently, employing the versatile OpenBabel framework within a Linux environment, we adroitly transformed these 2D structures into their three-dimensional (3D) counterparts in SDF format. In pursuit of leveraging AutoDock Vina,9 we seamlessly transmuted these 3D structures into pdbqt format. A meticulous compilation of the ligand names was documented within a text file, serving as a crucial prerequisite for their deployment in the AutoDock Vina virtual screening tool.

Preparation of receptor

The focal point of our analysis was the formidable receptor structure, designated as the human TLR4 and hagfish (PDB ID: 2Z63.pdb). Procuring the requisite protein data bank (PDB) structure, we painstakingly orchestrated a series of meticulous adjustments through the utilisation of MGL Tools V.1.5.7. These meticulous modifications encompassed the judicious rendering of hydrogen bonds, polar charges and a meticulous definition of the grid box, imperative for facilitating the optimal positioning of the ligands. Pertaining to the latter, the discernment of the optimal coordinates within the x, y and z axes was a crucial determinant. This judiciously ascertained information was instrumental in generating the indispensable configuration file, mandatorily required for the ensuing Vina docking operations.

Molecular docking

Equipped with the meticulously prepared receptor, ligands meticulously converted into pdbqt format and the all-important configuration file at our disposal, we initiated molecular docking process by using AutoDock Vina.9 The grid file was introduced with XYZ 8.631347, –33.659551 and −12.366184 coordinates, respectively. The prepared ligands were docked in the binding site of 2Z63.pdb (figure 1) to evaluate the binding affinities of these compounds against the receptor structure. The validation of docking was performed by redocking.

Figure 1
Request permissionrequest permission icon
Figure 1

Overview of binding sites on human TLR4 and hagfish (PDB ID: 2Z63.pdb).

Molecular dynamics

To study the dynamic behaviour and stability of protein–ligand complexes, molecular dynamics (MD) simulations were performed using the Desmond module of the Schrödinger software. First, the protein–ligand complex obtained from the docking step was prepared for MD simulations. The structures were solvated in an explicit water box, ensuring that the entire complex was surrounded by water molecules10 to mimic the physiological environment. Counterions are added to neutralise all charges of the system. The system energy reduced to eliminate any unfavourable contacts and relax the initial structure. Following energy minimisation, the system underwent a stepwise equilibration process to allow the system to reach a stable state. The equilibration was performed in two stages: NVT (constant number of particles, volume and temperature) and NPT (constant number of particles, pressure and temperature) ensembles. During the NVT equilibration, the system temperature gradually increased from an initial temperature to the desired simulation temperature. The protein and ligand atoms were harmonically restrained to maintain the overall structure while allowing the solvent and ions to adjust. In NPT equilibration, the system pressure was gradually equilibrated. Here, the position restraints on the protein and ligand were further reduced to allow for greater flexibility in the system. The pressure was maintained at the desired value using the Martyna-Tobias-Klein barostat,11 which adjusted the system volume to maintain the set pressure. Once equilibrated, the production MD simulations were initiated for 100 ns to capture sufficient dynamics of the protein–ligand interactions.12 The simulations were performed under periodic boundary conditions to simulate an infinite system, allowing molecules to move freely. The all-atom optimized potentials for liquid simulations (OPLS-AA)force field13 was employed to describe the protein and ligand parameters during the simulations. Throughout the MD simulations, the coordinates and energies of the system were recorded at regular intervals. The resulting trajectories were then analysed to gain insights into various aspects of the protein–ligand complexes.

Results and discussion

Docking studies

A meticulous analysis of the obtained results ensued, wherein we judiciously identified the ligand boasting the most diminutive binding affinity energy value among the compounds under scrutiny. Employing meticulous precision, we collated and meticulously recorded these values within an exclusive file, thereby facilitating their facile interpretation. Intriguingly, on meticulous investigation, we uncovered an inherent structural asymmetry within the drug’s size, starkly evident within its own 3D pdbqt file. Post the exclusion of ‘Lactoferrin’, our exacting analysis revealed an elite cohort of 15 compounds, each distinguished by binding affinity energy values ranging from a formidable −5.1 to a staggering −8.7 shown in online supplemental table 1. Superimposed image of docking and redocking of celastrol with target protein receptor TLR4 (PDB ID=2Z63) is represented in online supplemental figure 1

3D interactions

The 3D interactions as mentioned in figure 2 of celastrol in complex with human TLR4 and hagfish reveal that the compound is forming two hydrophobic bonds with binding site residues Glu254 and Ans227 and one H-bond with Glu229. Lumacaftor is forming one hydrophobic (Glu225) and two H-bonds (Leu283 and Thr284) with the active site residues of human TLR4 and hagfish. The most interactions were found with nilotinib. Nilotinib formed Six H-bonds and five hydrophobic bonds with Val310, Asn309, Glu286, His229, His256, Ile285, Thr284, Glu254 and Arg227. These bonds determine the binding affinity and selectivity of these ligand for human TLR4.

Figure 2
Request permissionrequest permission icon
Figure 2

3D interaction diagram of top molecules in complex with 2Z63.pdb. 3D, three dimensions.

MD simulation

In this study, we performed computational analysis to investigate the interactions between three herbal compounds, namely nilotinib, lumacaftor and celastrol against human TLR4 and hagfish (PDB ID: 2Z63.pdb). To start the simulations, we used the docking method to obtain an initial binding poses of the compounds in the protein body. Various analyses, such as root mean square deviation (RMSD), root mean square fluctuations (RMSF), protein–ligand interactions, ligand profiling, were done post MD simulation for comprehensive analysis.

2Z63/celastrol complex

2Z63/celastrol complex presents the results of an MD simulation depicting the complex behaviour of the protein in compound within its allosteric binding pocket. The simulation was conducted over a duration of 100 ns. Throughout the simulation, the protein remained bondedly stable in ligated shape, and really fewer primary deviations in RMSD have been determined in the course of simulation as shown in figure 3. The RMSD values of the protein and the ligand ranged from 1.8 to 4.8 Å and 3.2 to 5.6 Å, respectively, indicating that they are all stable. The RMSF analysis (figure 3) shows that protein residues fluctuate in the range from 3.0 to 5.4 Å, with the best peaks occurring in the loop area and N-terminal and C-terminal zones, which expected greater flexibility due to their non-interactive nature. The binding strength of the ligand to the protein is indicated by the low RMSF values observed for the binding sites. Secondary structure analysis shows that the α-helices and β-sheets in the protein are inflexible compared with the unstructured region, where fluctuations primarily exhibited at the loop surfaces. The residues with higher peaks in the RMSF graph corresponded to these surface’s loop and terminal areas, as determined by MD simulation shown in figure 3. The simulation analysis also investigates the changes of secondary structure elements (SSEs) during the interaction of the 2Z63/Celastrol complex.14 Online supplemental figure 2 illustrates the changes in SSEs donated by the residual index over 100 ns simulation time. The percentage of SSEs remained constant at 24.34%.

Figure 3
Request permissionrequest permission icon
Figure 3

(A) RMSD spectra of the 2Z63.pdb represents the RMSDs per frame (100-ns) of the complex model trajectories of 2Z63/Celastrol (B) RMSF plot of residual C atom of 2Z63.pdb shows the fluctuations that come about when Celastrol occupy the 2Z63.pdb binding pocket. Ligand interaction properties show its interaction ratio (C) and interaction strength (D) of Celastrol, with residues 2Z63.pdb. RMSD, root mean square deviation; RMSF, root mean square fluctuation.

During the 100 ns simulation, the amino acid residues in the allosteric region of the protein participated in series of consistent interactions with the compound, as shown in figure 3. H-bonds, in particular, play a crucial position in ligand binding due to their influence on drug selectivity, metabolism and adsorption. The residue Glu376 was determined to have tremendous H-bond interactions at some stage in the 100 ns simulation. Lys402 and Ile450 formed hydrophobic bond with the compound for over 75% of the simulation and it also formed water bridges with His426 for 75%. Tyr403 and Tyr451 formed bond for less than 25% of the simulation with the compound. From the protein ligand contact (figure 3), it is evident that Glu376 formed hydrogen bond for 99% of the simulation time with hydroxyl group attached to the aromatic group, Lys402 formed pi-cation for 59% with the aromatic ring and His426 formed water bridge for 71% with the -OH of carboxylic acid group of Celastrol. The ligand torsions plot (online supplemental figure 3) describes the conformational changes of each rotational bond in the ligand during the 100 ns simulation. The plot diagram shows the 2D structure of the ligand with colour-coded rotational bonds, accompanied by dial plots and bar plots representing the torsion angles for each bond. The ligand properties such as RMSD, radius of gyration (rG), molecular surface area (MSA), solvent accessible surface area (SASA) and polar surface area (PSA) are given in online supplemental figure 4, which are found within the range 0.2–0.6 Å, 4.32–4.48 Å, 378–387 Å2, 300–450 Å2 and 150–165 Å2, respectively.

2Z63/nilotinib complex

In figure 4, the RMSD plot of the 2Z63/nilotinib complex shows the stability of the protein over 100 ns of simulation, with first-order changes over approximately 4. Throughout the simulation, the protein remained unstable in ligated shape, and important deviations in RMSD observed throughout simulation, this instability was observed after 10–15 ns and of the simulation length. Analysing the RMSF of C-alpha atoms of protein residues found out ordinary changes, similar to those observed in the apo-form. However, specific residues around positions 250 and 400 exhibited higher RMSF values, reaching up to 4.0 Å. Throughout the simulation, the secondary structural elements (SSEs) of 2Z63.pdb, inclusive beta strands and alpha helices, were monitored. Online supplemental figure 5 illustrates the changes in SSEs over time, denoted with the aid if residue index. The evolution of the SSE percentage for every residue in 2Z63.pdb can be observed.

Figure 4
Request permissionrequest permission icon
Figure 4

(A) RMSD spectrum of 2Z63.pdb represents the average RMSD for each stem (100-ns) of the 2Z63/nilotinib complex model trajectory. (B) RMSF plot of C atoms of residues 2Z63.pdb Changes vary. Shown is what occurs when nilotinib occupies the 2Z63.pdb binding pocket. Ligand interaction properties of nilotinib; It shows the conformational ratio (C) and interaction strength (D) with residues of 2Z63.pdb. RMSD, root mean square deviation; RMSF, root mean square fluctuation.

In the case of the 2Z63/nilotinib complex, the residues Asn409, His431, Asn433 and His458 have been determined to be most essential in terms of H-bond interactions and these formed H-bonds for more than 50% of the simulation time. The stacked bar graphs in figure 4 constitute the duration of each specific interaction. Phe408 formed Pi-Pi staking and Pi-cation with Nilotinib for 32% and 39%, respectively. There is a greater number of bonds formed in the 2Z63/nilotinib complex but the interaction time is less than 25% which indicates the disability of nilotinib in complex with 2Z63.pdb.

The simulation starts from the centre of the radial view, and the evolution of the torsion angle is depicted radially outward. The result of the torsion angle (online supplemental figure 6) is shown using bar graphs, showing the complete data on the dial plot. Additionally, if torsional potential data are available, the graph shows the potential energy associated with the rotational linkage (and the associated torques). Potential values are shown on the left y-axis and are measured in kcal/mol. Ligand properties such as RMSD, rG, MSA, SASA and PSA presented in online supplemental figure 7 are in the range of 1.5–4.5 Å, 1. 5.5–6.5 Å and 480–1, respectively, 496 Å2, 450–900 Å2 and 105–135 Å2.

2Z63/lumacaftor complex

In figure 5, the RMSD plot of the 2Z63/lumacaftor complex illustrates that the protein-maintained stability in its ligated form throughout the simulation, with minimal deviations in RMSD. RMSD values for the protein and the ligand ranged from 1.0 to 3.5 Å and 5 to 8.5 Å, respectively, demonstrating their overall stability. However, protein and ligand become unstable at 80 ns. This pattern suggests that the ligand remained tightly bound throughout simulation. Peaks in the RMSF (figure 5) plots represent protein regions exhibiting the highest conformational changes during the simulation. In contrast, the tails of the protein (N and C termini) tend to change more frequently. Residues around 250–300 displayed the most significant RMS changes even when not in contact with ligands. This area corresponds to circular area. It has also been found that components such as β-strands and α-helices are generally stiffer than the disordered regions and show less flexibility than the loop region. Protein during simulation the SSEs, including alpha helices and beta strands, was monitored. Online supplemental figure 8 displays the distribution of SSEs along the protein structure according to residue index. Some α-helices and β sheets existed for certain periods and then transitioned into loop regions of the secondary structure. This continuous change was observed during 100 ns simulations, demonstrating the change in this protein and its ability to interact with other proteins in the cell.14

Figure 5
Request permissionrequest permission icon
Figure 5

(A) RMSD spectra of the 2Z63.pdb represents the RMSDs of each frame (100-ns) of the sample trajectories for 2Z63/Lumacaftor complex model orbital (B) RMSF plot of residual C atom of 2Z63.pdb shows the variable changes when 2Z63.pdb binding pocket is occupied by lumacaftor. Ligand interaction properties of Lumacaftor show the interaction coefficient (C) and interaction e-energy (D) with 2Z63.pdb residues. RMSD, root mean square deviation; RMSF, root mean square fluctuation.

Figure 5 depicts the contacts formed between protein and lumacaftor during the simulation. Key H-bond interactions were identified, with Gly120 and Asn143 playing pivotal roles. Hydrophobic bonds formed by Tyr170 and Pro145 occurred for less than 30% of the simulation time while water bridges involving Asn143 and Ala118 persisted for more than 40% of the simulation time with lumacaftor. The ligand torsions plot (online supplemental figure 9) gives bits of knowledge of each rotational bond change in the compound during the simulation force provided. In online supplemental figure 9, the top panel presents a 2D representation of the ligand with colour-coded rotatable bonds. Each rotatable bond torsion is accompanied by a dial plot and bar plots of the same colour.14 The dial plot shows the torsion angle changes extending radially from the centre during the simulation. The bar chart shows the desired value for the torsion angle obtained from the dial plot chart. Histograms and torsional potential correlations provide information about the conformations a ligand undergoes to maintain its shape when bound to a protein. The ligand properties presented in online supplemental figure 10 such as RMSD, rG, MSA, SASA and PSA are found to be in the range of 1–3 Å, 4.8–5.4 Å, 396–404 Å2, 280–400 Å2 and 165–178 Å2, respectively.

Conclusion

This exhaustive endeavour, characterised by the astute employment of Auto Dock Vina, has witnessed the meticulous virtual screening of 30 compounds against the backdrop of the formidable ‘crystal structure of the TV8 hybrid of human TLR4 and hagfish’ (PDB ID: 2Z63.pdb). Through the attainment of the top nine compounds, characterised by remarkable binding affinity energies, our endeavours have shed unprecedented light on potential drug candidates boasting a commendable ability to seamlessly interact with the aforementioned receptor structure. Future explorations and incisive analyses stand poised to harness these insightful findings, potentially fostering ground breaking advancements in the dynamic realms of drug discovery and development. From the MD data, we can conclude that among the three complexes 2Z63/Celastrol complex is showing stability for the 100 ns of simulation time with very few fluctuations whereas 2Z63/Nilotinib has the least stability as the ligand RMSD reached to 32 Å in the 100 ns simulation time. This may be because it was close to amino acids constituting the active site, and because of that it may result in unstable complex as MD simulations try to find the best position by rearranging the amino acids. Glu376 of TLR4 is specifically important for the recognition and binding of LPS, a component of the outer surface of Gram-negative bacteria and in our study celastrol is binding with Glu376 which indicates that celastrol may inhibit the binding of LPS to TLR4 thereby preventing the dry eye syndrome. The results from the in silico data indicate that celastrol can be formulated for the further studies in the treatment of dry eye syndrome.

Overview

Article metrics
Altmetricopen-url
Dimensionsopen-url