Exploring the interaction between 3-D structure of TLR 9 and prostaglandin analogues

1.1 TLRs and autoimmune diseases

As TLRs detect different molecular patterns of pathogens, they may also recognize endogenous proteins and nucleic acids as foreign particles. It is speculated that TLRs’ function in autoimmunity is the consequence of their continuous activation or dysregulation, as TLR-mediated pathways Eventually, pro-inflammatory cytokines are produced. and up-regulation of antigen-presenting cells (APCs) (Hosseini et al. 2015). Consequently, TLRs play important part in the pathogenicity of coronary arterial disease (CAD), systemic Lupus erythematosus (SLE), Rheumatoid Arthritis (RA), and Psoriasis (Frantz, Ertl, and Bauersachs 2007; Hosseini et al. 2015).

The pathogenesis of psoriasis (formation of psoriatic plaque) begins with binding a small peptide, such as Cathelicidin LL37 (Scheenstra et al. 2020), to damaged nucleic acids. Subsequently, it stimulates TLR9 (Act via MYD88 and TRAF6) on plasmacytoid dendritic cells (pDCs) that results in the production of interferons (IFN-α and IFN-β), interleukins (IL-12p40), and tumor necrosis factors (TNF) which induce the auto-inflammatory responses (Takeuchi and Akira 2010). Therefore, TLRs could be considered as promising targets in the treatment of autoimmune disorders.

TLRs influence innate and adaptive immune responses. Self-components activate TLRs inappropriately, causing sterile inflammatory and autoimmunity. Numerous processes cause autoimmunity, including autoreactive immune cell subsets and immunological tolerance loss (Tabeta et al., 2004). Genetic and environmental variables that lack adaptive immune response control to self-antigens cause organ-specific autoimmune disorders (Patra et al., 2020). Organ-specific autoimmunity, type 1 diabetes, and Crohn's disease patients had PRR overexpression (Basith et al., 2011). PAMPs in tissues after infection have been associated to autoimmunity in various studies. These are among many studies showings that PRR stimulation on innate immune cells by PAMPs or pathogens can dysregulate self-tolerance and activate autoreactive T- and B-cells. TLRs also recognize host-derived endogenous ligands that have changed or accumulated in non-physiologic compartment (Vidya et al., 2017). DAMPs, such as HMGB1, saturated fatty acids, and amyloid β, are produced from injured tissue or apoptotic cells and can cause chronic or acute inflammation. TLRs can bind to self-molecules and cause autoimmune diseases (Reuven et al., 2014.). Table 1 lists endogenous ligands and TLRs.

1.2 Prostaglandins and TLRs modulation

Prostaglandins (PGs) are lipid autacoids involved in biological processes, such as inflammation, blood flow, and fertility (Ricciotti and FitzGerald 2011). The relationship between PGs and TLRs has been established, discussed, and demonstrated in some literature. Some studies revealed that PGs could modulate the expression of TLR2 on immune cells and control inflammatory events in the brain. 15-deoxy-Delta12,14-PG J(2) (15d-PGJ(2)), 15d-PGD(2), and PGD(2) downregulate TLR2 expression leading to the interference of brain inflammation (Yoon et al. 2008). Prostaglandin E2, which limits mRNA translation, also decreases TLR4 expression on alveolar macrophages (Degraaf et al. 2014). In a recent study on developing cancer vaccines, PGE2 in combination with intracellular TLRs and TLR3/7/8 agonists were used in the DC maturation cocktail. This mixture improved DCs maturation and cytokine production but inhibited DCs-mediated antigen presentation to CD8 + T cells (Gierlich et al. 2020), showing that while prostaglandins facilitate maturation of DCs, they abolish DCs’ antigen-presenting ability, probably through direct interaction with intracellular TLRs.

Therefore, it has been proved that PGs have a significant role in the modulation of TLRs; however, the molecular mechanism of PGs activity on intracellular TLRs regulation, such as TLR9, seems to be via direct PG-TLR interaction. Accordingly, PGs could render therapeutic effects on autoimmune diseases through direct blockage of intracellular TLRs.

1.3 The crystal structure of TLR9

The crystal structure of human TLR9 is not present in any of the relevant databases. Nonetheless, the conserved functional sites shared by various species (bovine, mouse, and horse) have made modelling human-TLR9 simple, allowing the molecular interaction to be approximated by homology. The sequence similarity between horse-TLR9 and human-TLR9 is 84 percent, making it an ideal template for modelling an accurate three-dimensional protein structure (Umeharu Ohto et al. 2015).

CpG-DNA primarily interacted with horse-TLR9 at the leucine-rich repeat (LRR) region that extended from the LRR N-terminal domain (LRRNT) to the LRR10. Therefore, targeting this region can stop the activation of immune cells and block the immunomodulatory effect of TLR9. LRR2 (residues 87–110) LRR5 (residues 167–190), LRR8 (residues 243–268), LRR11 (residues 333–356), LRR18 (residues 545–572), and LRR20 (residues 600–622) generate an insertion loop that is essential for both recognising nucleic acid fragments (LRR2, 5, and 8) and involvement in protein–protein interactions (LRR2, 5, 8, 11, 18, and 20) (Fig. 1). Umeharu Ohto et al. have validated the role of the inserted loop by introducing a mutation that reduced or eliminated the function of TLR9, demonstrating that this region is essential for protein–ligand interaction (Ishida et al. 2018). The ability to exist as a monomer and a heterodimer is contingent on binding to CpG ligand. TLR 7–9 have a long insertion loop (Z-loop) between LRR14 and LRR15, unlike other TLRs. According to numerous studies, the Z-loop is in charge of the proteolytic processing required for the formation of functional TLR9. Proteolytic processing of an insertion loop (Z-loop) is required for homodimerization after binding to an unmethylated CpG ligand, which results in its activation. (Tanji et al., 2015).

Close

Fig. 1

In the modelled human TLR9, the important LRRs regions, such as insertion loops and CpG DNA contact regions, as well as the Tyr132 and Tyr208 residues, which have been mentioned in uniport for their important role in binding activity, appear bluish. In regions 1–30 (splices variant), isoforms 1 and 2 appear to have different sequences, which is beyond the scope of this study. The transmembrane domain (residues 819–839) and the cytoplasmic TIR domain (residues 840–1032) that interact with MYD88 to cause cytokine secretion will also be excluded from the study. To the right, the conserved region after alignment of TLR9-human in orange and TLR9-horse in grey (liganded form with CpG DNA in red and blue).

Export to PPT

1.4 Generation of a functional receptor

In macrophages and dendritic cells, TLR9 normally exits the endoplasmic reticulum and travels to endolysosomes. TLR9′s ectodomain is cleaved in the endolysosome (Z-loop proteolytic cleavage), leaving no full-length protein in the chamber where ligand is recognised. Although both the complete and sliced forms of TLR9 can bind ligand, only the processed form activates MyD88, implying that the fragment receptor, rather than the complete form, is functional. In addition, situations that inhibit receptor proteolysis, such as forced TLR9 surface localization, make the receptor inoperable. We believe that the ectodomain cleavage of TLR9 is a strategy for limiting receptor activation and preventing TLR9 from responding to self-nucleic acids. (Ewald et al. 2008).

2.1 Homology modelling and binding site studies

The crystal structure of the extracellular domain of TLR9 was modelled by the SWISS-MODEL web tool (Schwede et al. 2003), horse-TLR9 was the template for the modelling process (https://www.uniprot.org/uniprot/Q2EEY0). The modelled structure of TLR9 was validated upon matching with the horse-TLR9 structure carried by the Needleman-Wunsch algorithm and blossom-62 matrix on the UCSF Chimera tool (Pettersen et al. 2004). Further verification of the modelled protein was conducted by the PROCHECK server, the ERRAT2 server, and the molecular dynamic (MD) simulation (Laskowski et al., 2006; Schreiner et al. 2012; Hospital et al. 2015). Chi1-Chi2 and Ramachandran plots (Cowell et al. 2004) were also used to determine the adjacent residues’ torsion angles and predict the possible geometry of the secondary structure for all residues.

Three different strategies were used to choose the binding site. First, by looking for conserved binding sites in other proteins (especially horse TLR9) (Fig. 1); second, by using the FT Map webserver (Kozakov et al. 2015) to estimate the ability of residues to form hydrophobic and H-bond interactions (Fig. 2); and third, by performing multiple dockings against all of the protein's possible hot spots.

Close

Fig. 2

The chart is a demonstration of the contact rate of protein residues that translated into prediction of H-bond interaction and non-H-bond interaction that is done by FT Map tool for small molecule mapping service.

Export to PPT

2.2 Molecular docking

The protein and ligands were prepared for docking by adding hydrogen atoms and charges. Charges were modelled by AMBER ff14SB for the residues and AM1-BCC (Austin model with bond and charge correction) for the ligands (Maier et al. 2015; Jakalian, Jack, and Bayly 2002). Checking on missing residues (missing structures and loops) was examined by the Salilab modeller (https://salilab.org/modeller/). The docking process and estimated binding affinities were handled by the FlexX Docking tool that is built-in to SeeSAR software (https://www.biosolveit.de/products/seesar/). Up to 100 poses were produced for each ligand (Irwin and Shoichet 2005; Wishart et al. 2006; Kim et al. 2016; https://enamine.net/) and the docking results were confirmed via re-docking the ligands. Visualization was carried out by Chimera and Discovery Studio (Biovia, 2021).

2.3 Molecular dynamic simulation

. The MD simulation process was carried out for the ligand–protein complexes that had given the best-estimated affinity scores. The GROMACS 2018.2 package (Van Der Spoel et al. 2005) was used along with Charmm27 to parametrize the residues. The system was solvated by applying the TIP3P water model (Mark and Nilsson 2001). Swissparm (Zoete et al. 2011) were employed to generate the topology and parameters of three ligands: Bimatoprost (BI), Butaprost (BU), and Misoprostol (MI). A triclinic periodic box (periodic boundary condition) was generated, and the distance between the edge of the box and the complex of protein–ligand was set to 1.0 nm. Cl-and Na ions (0.15 M) were added to achieve a neutral concentration. Energy minimization was run using the steepest descent to make the system's potential energy reach −1.4*106 kJ/mol. The energy step size was set to 0.01 nm, and a maximum of 50,000 steps was introduced. The temperature represented by the NVT (isochoric isothermal) ensemble was equilibrated to 310 K. The time of equilibration was fixed to 1 ns with a time step of 2 fs. Afterward, NPT (isothermal-isobaric) equilibration was processed for 1 ns with a time step of 2 fs. Further settings of the NPT ensemble were 0.1 nm for both Van der Waals (VDW) and electrostatic interactions. The Modified Berendsen thermostat for temperature coupling (310 K), the Particle Mesh Ewald (PME) for calculating the long-range electrostatics, and the Parrinello-Rahman as a Barostat were also used. The MD simulations were performed for 100 ns with the same settings mentioned for NPT. The free binding energy, VDW interactions, electrostatic interactions, solvent accessible surface area (SASA), the contribution energy of the residues, and the polar solvation energy were carried by the MM/PBSA method (Genheden and Ryde 2015). The root means square deviation (RMSD), root mean square fluctuation (RMSF), the radius of gyration (Rg), and the principal component analysis tools that are built-in to GROMACS software were used to monitor the conformations of the protein. (Fuglebakk, Echave, and Reuter (2012); Lobanov, Bogatyreva, and Galzitskaya (2008)).

3.1 Homology modelling and binding site results

The Procheck server has suggested a high-quality modelled 3D structure. The overall quality factor was 79.15, which is very close to the quality of the template (horse-TLR9 scored 77.34). The modelled protein was also considered qualitatively high based on Ramachandran, and Chi1-Chi2 plots. Tyrosine’s (hydrogen bond donor–acceptor and π-π interaction), prolines (hydrophobic interaction and H-bond acceptor), arginine’s (hydrogen bond acceptor/donor), and methionine’s (hydrophobic interaction) were involved in the formation of the binding site (Erusappan et al., 2021). The FT Map results regarding probing the protein hot spots (Fig. 2) propose that the region between LRR5 and LRR8 (residues 167–268) is favourable for forming both H-bonds interactions and non-H-bonds interactions and can be a target to form stable and bounded complexes (Kozakov et al. 2015).

Furthermore, the binding pocket is composed mainly of two cavities with a specific area linked to each other to form a space for occupation (Fig. 4). In general, ligands characterized by extended linear structures could have the advantage of occupying the binding site, with some consideration for the number of hydrogen bond acceptors and donors.

3.2 Molecular docking

Thousands of ligands were docked into the binding pocket of TLR9 to conduct virtual screening. However, the prostaglandin analogues were uniquely producing remarkable interactions with well-estimated affinities (Fig. 3) for the proposed biding site; the results are not a coincidence since many studies have reported prostaglandins’ vital role in influencing the TLR family. Furthermore, the high structural similarity of prostaglandin analogues, including a central ring with two extended hydrocarbon arms, led to close docking results, which can verify the docking process and confirm that prostaglandin analogues are likely to be promising blockers of human-TLR9 (Tu et al., 2022).

Close

Fig. 3

An illustration of the prostaglandin analogues structures with their estimated affinity in a range of micromolar.

Export to PPT

Close

Fig. 4

Visualization and prediction of H-bonds and non-H-bonds carried by Chimera and Discovery Studio. The binding site resembles a groove that made up from the connection of LRR5, 6, 7, and 8.

Export to PPT

As shown in Fig. 4, the phenyl group of BI was the most interactive molecule among the other two compounds because it penetrated deeply into the binding pocket to form a difficult-mediated interaction with Pro269, Val233, and Arg231. When combined with nearby amino acids like Arg185 and Cys184, the amide terminal's interactive oxygen could form H-bonds (Table 2). Overall, BI formed 5H-bonds with Tyr229, Met266, Arg185, and Cys184; however, Tyr229′s H-bond and one of Met266′s failed to withstand the 100 ns MDs (Fig. 5). Tyr208 has successfully undergone van der Waals interaction with BI, a crucial residue in the binding site. BU had the second-best outcome as a result of the terminal hydrophobic moieties that cyclobutane represented and a propyl group that made favorable hydrophobic interactions with Met266, Tyr208, and Arg231. The terminal ester group simultaneously formed H-bonds with Arg185, Cys184, Arg185, and Cys178 (2 bonds). None of these H-bonds remained strong during the MDs, with Arg185having the weakest bond (dos Santos Nascimento et al., 2022). It is thought that BU's loss of binding energy is caused by a collision with Gln186 (the oxygen atom of the ligand has clashed with the nitrogen atom of the residue). However, BU's affinity did not wholly rely on H-bonds, but VDW, electrostatic, and hydrophobic forces also showed remarkable engagements in maintaining the respective ligand–protein interaction. With a score of 0.153–15.2 µM and the development of around four H-bonds with Cys184 and Arg185 (3 bonds), MI took third place among estimated affinity values. The H-bonds with Cys184 and one with Arg185 remained reasonably strong during MDs. Tyr180 is predicted to interact hydrophobically with the terminal carbon chain of MI, while the introduction of H-bonds by 3-hydroxycyclopentanone (the core) and its terminal ester enhanced protein–ligand interactions. It is believed that the prostaglandins' extended linear structures, which range in length from 14 Å to 16 Å, will have a positive effect on how well they bind to the protein.

Close

Fig. 5

A) The RMSD calculated for the backbone of Apo and the complexes receptor-BI, MI, and BU. B) The per-residue RMSF calculated for the backbone of Apo and the complexes receptor-BI, MI, and BU. C) SASA calculated for the protein of Apo and the complexes receptor-BI, MI, and BU. D) Radius of gyration calculated for protein of Apo and the complexes receptor-BI, MI, and BU. E) H-bonds distance of receptor-BI complex. F) H-bonds distance of receptor-BU complex. G) H-bonds distance of receptor-MI complex.

Export to PPT

Prostaglandin G2′s estimated low affinity could be primarily attributed to the 2,3-dioxabicyclo [2.2.1] heptane moiety's (Fig. 3) improper interaction with the binding site. As a result, it changed the molecule's overall conformation, lengthening the structure to 9.4 and reducing the interaction spectrum. Due to the lack of the -bond in the hydrocarbon arms and the substitution of the hydrogen donor–acceptor hydroxyl group for the ketone group, Unoprostone's structure flexibility had a negative impact on the interaction with the receptor.

3.3 Molecular dynamic simulation results

MD simulations were run separately for compounds BI, BU, and MI combination with TLR9 to discover their possible interactions and conformations in the designed physiological environment. Upon taking the last 100 frames from the MD simulation to calculate the free binding energy by MM/PBSA method, the VDW interactions were the dominant intermolecular forces supporting the interaction between ligands and the receptor (Table 3). Despite the considerable resistance against the binding of prostaglandins with the protein caused by the polar solvation energy, VDW forces extensively participated in supporting the total binding energy for BI, MI, and BU by creating the energy amount of −190.5 kJ/mol, −155.6 kJ/mol, and −135.6 kJ/mol, respectively. The results regarding solvent-accessible surface area (SASA) energy, the energy of binding between solvent and the accessible surface of solute, were encouraging, as shown in Table 3.

The RMSD calculation has revealed the Apo-protein has more conformational change than Holo-protein upon occupation of the ligands, and the least conformation change was detected by the complex with BI which cause the receptor to be rigid and more stable, the RMSD graph for BI was fluctuated around 2.1 Å while the Apo form around 3.5 Å and does not seem the system in case of apo will reach convergence soon (see Fig. 5). The RMSF calculation demonstrated a high fluctuated region represented by LRR14 and LRR15 which collectively represent the residue sequence section 415–494 that where Z-loop occur (see Fig. 5). The folding of the protein is increasing slightly in the existence of BI and BU as observed by SASA calculation and support the RMSD in that Apo form has more conformational change than others which then empowered by the compactness calculation which is conducted by radius of gyration which is showing the Apo form to have more fluctuation in compactness (unfolded) (see Fig. 5).

For the backbone of apo and hollo protein, the principal component analysis or the essential dynamics of the protein were identified and calculated (Fig. 6). It is a procedure to track the protein's motion in order to comprehend its biological purpose. PCA assists in shortening the MDs trajectory to locate protein configurational spaces with a few degrees of freedom where anharmonic rhythms can occur. In the simulation, the conformational dynamics of the samples (Apo, BI, BU, and MI) each of which has slightly different conformational behavior, as illustrated by the projections of the backbone coordinated onto the plane described by the first principal plane (the first two eigenvectors). Wide scattering on the PCA plane (Fig. 6B), particularly on the first eigenvector (x-axis), suggests that the apo-protein has the most flexibility, whereas the second eigenvector (y-axis) is largely overlapping with other complexes; at the molecular level, the Z-loop is primarily responsible for these differences; except for Apo-protein, the entire protein structure fluctuates in a regular pattern. Without a shadow of a doubt, this kind of motion (the flexibility of the Apo-form) can be considered as important for it to carry out its normal function and interact with other proteins that are relevant within the compartments of the cells. High fluctuation up to 6 Å of the Z-loop appeared on the RMSF plot (Fig. 6A) of the eigenvector one downwardly upon binding of MI followed by BI and BU and lastly by apo-protein (so, the high scattering on Fig. 6B is not due to Z-loop only). In vector two the fluctuation up to 4 Å of the Z-loop continued for BU and BI, whereas MI has new regions of fluctuation account for LRR1-LRR5. Fig. 6c shows 50 aligned frames for each sample of eigenvector one and two. The difference in the behavior of the Z-loop when comparing the complexes and the Apo-form is evidence of a change in the level of activation of the receptor. This is because the Z-loop is responsible for the proteolytic activity. It should be noted that eigenvector one concentrated on the Z-loop fluctuation while eigenvector two concentrated primarily on the insertion loop fluctuation. Furthermore, a dynamic process of change has occurred at the CpG DNA binding site. In contrast to the apo-form and other complexes, the binding site of the BI-receptor has become more rigid, while the MI-regions receptor's LRR1–LRR6 (which represent the binding site), LRR18, and LRR20 (which represent insertion loops for protein–protein interaction) have fluctuated the most (Ullah et al., 2020).

Close

Fig. 6

A) The RMSF for the eigenvector 1 and eigenvector 2 for all the occupied and unoccupied protein. B) 2D projection of the trajectories. C) 50 frames of the complexes and the apo-protein for each vector where extracted.

Export to PPT