Results

Dynamics of apo- and CpG ssDNA and IM-bound TLR9 dimer complexes

Initial sequence homology characterization with other TLR9 ECDs, mapped onto the modeled human TLR9 structure, is shown in Figure 1B. The modeled TLR9_apo had an root-mean-square deviation (RMSD) value of 0.55 Å with mouse TLR9 and 0.45 Å with equine TLR9. The figure demonstrates the conservation of ligand-exposed and ligand-unexposed regions of the LRRs. We observed 26 LRR repeats in the human TLR9 model, similar to mouse and horse TLR9 structures previously observed.21,25 We then, prepared three different system, i.e., TLR9_apo, TLR9_CpG, and TLR9_CpG_IM and ran three independent all-atom simulations of 300 ns each (totaling to 900 ns for each system). Post simulation, the protein backbone RMSD was calculated to assess stability of the protein in the simulations. The backbone RMSD distribution of the three complexes show small deviations and unimodal distribution with average values of 0.23 ± 0.026, 0.25 ± 0.035, and 0.24 ± 0.029 nm for TLR9_apo, TLR9_CpG, and TLR9_CpG_IM, respectively, suggesting stability of the simulated complexes (Figure 2A). The stable deviation in TLR9_apo during the simulations suggests that TLR9 exists as a preformed dimer, at least in the ectodomain (ECD). This is consistent with a previous study showing that human TLR9 exists as a preformed dimer.31 The RMSD distributions of the ligand CpG in TLR9_CpG and of both ligands, CpG and IM, in TLR9_CpG_IM complexes were clustered at 0.35 ± 0.07 and 0.37 ± 0.06 nm, respectively. The minimal deviation in the ligand RMSD indicated stability and possible interactions with the TLR9 homodimer (Figure 2B).

Figure 2: Structural dynamics of TLR9 complexesDensity distributions of RMSD for TLR9 backbone (A) and both (ssDNA and IM) ligands (B). Density distribution of radius of gyration (Rg) for TLR9 backbone (C). Density distribution of solvent-accessible surface area (SASA) for TLR9 (D). RMSF of TLR9 backbone (E). The average values along with standard deviation are also given.

The radius of gyration (Rg) of the simulated complexes was assessed to evaluate the rigidity and compactness of the protein complexes during the simulation period. The average Rg values of the TLR9_apo, TLR9_CpG, and TLR9_CPG_IM were 3.8 ± 0.01, 3.8 ± 0.008, and 3.8 ± 0.01 nm, respectively. The Rg of TLR9_apo, TLR9_CpG, and TLR9_CpG_IM complexes showed a unimodal distribution and had a broader distribution of Rg values (Figure 2C). We also calculated the solvent accessible surface area (SASA) to determine any major folding or conformational changes. The average SASA value was 636 ± 6.7, 638 ± 6.57, and 634 ± 6.35 nm for TLR9_apo, TLR9_CpG, and TLR9_CpG_IM, respectively, indicating no major conformational events (Figure 2D). Further, the per-residue fluctuations of the protein backbone were calculated with the root-mean-square fluctuation (RMSF) tool. The RMSFs of all three complexes were similar, except for parts of the Z loop (441–469 aa) region as well as neighboring residues of Z loop of each TLR9 monomer (Figure 2E). Notably, the Z-loop fluctuation was approximately 2-fold higher in the TLR9_CpG and TLR9_CpG_IM complexes than in TLR9_apo. Z-loop dynamics are crucial for binding and adapting CpG and IM to TLR9,32 and our results support this hypothesis. In addition, we observed similar per-residue fluctuations for both the monomers of TLR9, demonstrating a positive correlation between the dynamics of the monomers (Figure 2E). The simulations show that all TLR9 complexes remain structurally stable, with ligand binding primarily enhancing Z-loop flexibility, supporting its central role in accommodating CpG and IM while preserving the integrity of the preformed TLR9 dimer.

Hydrogen bonds and interaction energy of apo- and CpG ssDNA and IM-bound TLR9 dimer complexes

Furthermore, to check the intra- and inter-molecular hydrogen bonds within TLR9 complexes, we analyzed hydrogen bonds using a cutoff distance 0.35 nm and an angle of 30° between TLR9 monomers, as well as between TLR9 and CpG. The average hydrogen bonds and interaction energy observed between TLR9 monomers (TLR9 and TLR9∗) in TLR9_apo, TLR9_CpG, and TLR9_CPG_IM were 16 ± 5 and −1,537 ± 223 kJ/mol, 16 ± 4 and −1,591 ± 146 kJ/mol, and 13 ± 3 and −1,388 ± 139 kJ/mol, respectively (Figures 3A and 3B). From the results, the hydrogen bonds demonstrate the intactness of the TLR9 dimer in all simulated protein complexes. The distribution of H-bonds was broader for TLR_apo which reduced in TLR_CpG and TLR9_CpG_IM complexes. This heterogeneity in interaction may be a result of CpG and IM binding that reorganizes the dimer interface interactions. The interaction energy further suggests the stabilization of TLR9 dimer upon CpG binding. Although the hydrogen bonds remained similar in TLR9_apo and TLR9_CpG, there was a modest increase in their interaction energies which could be as a result of structural stabilization of TLR9 monomers upon CpG motif binding. By contrast, TLR9_CpG_IM had significantly fewer inter-molecular hydrogen bonds and reduced interaction energy between its monomers, indicating the reallocation of H-bonds between TLR9-TLR9∗ dimer to IM that binds in the dimer interface (Figure 3A).

Figure 3: H-bond analysisThe distribution of H bonds (A) and interaction energy (B) between TLR9 and TLR9∗. The distribution of H bonds between TLR9 and CpG ssDNA (2:2 stoichiometry) (C), between TLR9 and IM (2:2 stoichiometry) (D). Distribution of intra-H bonds within the TLR9 dimer (E). The average values along with standard deviation are also given.

The average number of hydrogen bonds between CpG and TLR9 (2:2) in TLR9_CpG and TLR9_CpG_IM complex was 27 ± 5 and 26 ± 6, respectively (Figure 3C). There was no increase in hydrogen bonds between CpG and TLR9 in TLR9_CpG_IM complex and TLR9_CpG complex which suggests that the binding of IM does not affect CpG binding to TLR9. In addition, the hydrogen bonds between TLR9 and IM from the TLR9_CpG_IM complex was around 20, indicating the strong interaction of IM with individual TLR9 monomers (Figure 3D). The intra-molecular hydrogen bonds within each TLR9 monomer did not vary much among the three complexes TLR9_apo, TLR9_CpG, and TLR9_CpG_IM (Figure 3E).

We analyzed the H-bond occupancy (in %) between TLR9 and TLR9∗ (∗ are residues in TLR9∗, hereafter) over 300 ns and observed H-bond between residues Q797∗-A736 (59%), R348∗-G565 (58%), Q797-A736∗ (53%), and H260-Q562∗ (52%) in TLR9_apo complex. For TLR9_CpG, we calculated Y180-E616∗ (70%), Q797-A736∗ (55%), Q797∗-A736 (54%), R348∗-G565 (49%), and K738∗-G791 (49%). The IM bound complex TLR9_CpG_IM showed Y180-E616∗ (66%), Q797-A796∗ (60%), Y180∗-E616 (60%), and H260-G563∗ (48%). We note that Y180-E616∗ had a maximum of 18% occupancy in TLR9_apo, while this pair in TLR9_CpG and TLR9_CpG_IM show more than 60% H-bond occupancy. Collectively, these results suggest that CpG and IM binding induce rearrangements within the ectodomain that may represent early structural events preceding downstream signaling, although our simulations cannot directly capture TIR-domain dynamics.

Principal components, free energy landscape, and dynamic cross-correlation matrix analysis

Following stability validation, we assessed variations in the global motions of TLR9 complexes, particularly those induced by CpG and IM binding. To this end, the principal components (PCs) of the concatenated trajectory using protein backbone atoms were analyzed. The top three PCs contributed substantially to the collective motions of all systems, accounting for ∼50% of the total variance (Figure S1). Phase-space projections (PC1 vs. PC2, PC2 vs. PC3, and PC1 vs. PC3) revealed periodic transitions separated by minimal energy barriers. Because the top PCs exhibited the largest amplitudes, porcupine plots were generated to visualize the directionality of motion (Figure 4). Across all complexes, the Z loop displayed the most pronounced conformational fluctuations, whereas the remainder of the ectodomain exhibited comparatively modest motions, consistent with the RMSF profiles described earlier (Figure 2E).

Figure 4: Conformational changes during simulationsPorcupine plots depicting the regions of the dominant motions from the top principal components (PC1, PC2, and PC3). Orange arrow vectors indicate the direction and relative magnitude of residue displacements, highlighting the dominant collective motion predicted by PCA.

We computed the free energy landscape (FEL) of the three simulated complexes to identify the meta-stable states of these ligand complexed proteins and characterize conformational heterogeneity. The FELs in Figure 5 show five, seven, and five low-energy populated states for TLR9_apo (Figure 5A), TLR9_CpG (Figure 5B), and TLR9_CpG_IM (Figure 5C), respectively. Visual inspection indicated that most free-energy differences originated from loop regions, which constitute a substantial portion of the receptor. The dimer interface was predominantly lined with hydrophilic and charged residues, along with a smaller number of hydrophobic (Leu) residues, which formed hydrogen bonds with backbone carbonyls of asparagine, histidine, and alanine (Figure 6; Figure S2). The increased number of metastable states in the TLR9_CpG complex likely reflects subtle conformational adjustments upon CpG binding, whereas the reduced number of states in the TLR9_CpG_IM complex suggests that IM binding stabilizes the CpG-bound dimer.

Figure 5: Free-energy landscapesFEL plots for TLR9_apo (A), TLR9_CpG (B), and TLR9_CpG_IM (C) using PC1 and PC2 coordinates. Corresponding structures from the metastable states are superimposed (middle). Representative structures from TLR9_apo, TLR9_CpG, and TLR9_CpG_IM are overlaid to show regions with variance.

Figure 6: TLR9 Dimer interfaceThe residues considered within 4 Å of each TLR9 monomer at the dimer interface TLR9_apo (A), TLR9_CpG (B), TLR9_CpG_IM (C) are shown.

We further analyzed the dynamic cross-correlation matrix (DCCM) plots of the three structures to determine the correlated motion of TLR9. The DCCM plots for TLR9_apo demonstrated a weak anticorrelation between the 100–300 aa region and 700–800 aa region of TLR9 monomers and a stronger anti-correlation between the 100–300 aa region of TLR9 and 700–800 aa region of TLR9∗ (Figure 7A). The anti-correlation suggests that N terminus of TLR9 moves in an anti-parallel fashion, weakly to C terminus of TLR9 (see red boxes in Figure 7A) and strongly to the C terminus of TLR9∗ (see orange boxes in Figure 7A). The anti-correlation is stronger at the C terminus TLR9∗ because of the cumulative motion of N terminus TLR9 and TLR9∗. There is positive correlation between the N terminuses of TLR9 and TLR9∗ which suggest that they move in the same vector direction (see yellow box in Figure 7A). These correlations are consistent with the preformed dimer arrangement reported previously and reflect ectodomain-level motions observed in our simulations.31 The results are also in line with the PCA analysis where the above-mentioned regions had opposite global motions in TLR9_apo and TLR9_CpG (Figure 4).

Figure 7: Dynamic cross-correlation matrix of all three systemsTLR9_apo (A), TLR9_CpG (B), and TLR9_CpG_IM (C). Map shows that the N terminus of TLR9 moves anti-parallel to the C terminus of TLR9 (red boxes) and even more strongly to the C terminus of TLR9∗ (orange boxes). The N terminuses of TLR9 and TLR9∗ display positively correlated motion (yellow box), indicating movement in the same direction. Schematic representation showing the possible monomer motion and how the DCCM relationships hints on a hinge-like interaction interface between the central regions of TLR9 and TLR9∗ (D).

In contrast, for TLR9_CpG, the anti-correlation observed earlier, between the 100–300 aa and 700–800 aa regions of individual monomers increased upon CpG binding while the anti-correlation with 700–800 aa of TLR9∗ reduced (Figure 7B). The N terminus positive correlation at TLR9 and TLR9∗ was observed in TLR9_CpG complex too. We also observed that the region downstream of Z loop was strongly correlated within the monomer after CpG binding. As expected, these correlations in complexes (TLR9_apo, TLR9_CpG_IM) were not observed for TLR9_CpG_IM complex, which might be as a result of relaxed dimer conformations upon IM binding (Figure 7C). The global motions observed in the PCA analysis corroborated the DCCM results and schematic view of the correlation/anticorrelation motions observed from DCCM were shown on Figure 7D. We also calculated the distance between the A765-A765∗ (C-ter LRR26) residues of the TLR9-TLR9∗ complex to analyze if the monomers moved away from each other at the C-terminal transmembrane domain. We did not observe any large changes, implying a stable dimer and TLR9_apo had heterogeneous distribution, followed by TLR9_CpG and then TLR9_CpG_IM (Figure S3). However, the decreasing trend in distance distribution and the decrease in heterogeneity might be a result of the monomers moving toward each other when CpG and IM bind to the dimer. Overall, CpG and IM binding reshape TLR9 dynamics by concentrating global motions in the Z loop, reorganizing long-range correlated movements, and stabilizing the dimer into fewer, more compact conformational states without disrupting its overall integrity.

Molecular interaction between TLR9 and TLR9∗, CpG DNA, and IM in metastable complex states

We further picked representative metastable states from each FEL cluster identified in the previous section and analyzed the inter-monomeric (TLR9-TLR9∗) and ligand protein interactions. These interactions provide critical insight into how ligand binding reshapes TLR9 dynamics. Our simulations revealed several residues not previously implicated in human TLR9 activation, including H612, H641, T642, and K690. In the apo state, inter-monomer contacts such as P105–H641∗/T642∗, H641–P105∗, T642–P105∗, and D618/H107-H107∗/D618∗ were observed in addition to hydrogen-bond interactions.

We then examined states 1, 2, 4, and 5 as representatives of the four FEL clusters in TLR9_CpG complex (Figure 8A). The 10 bp ssDNA was stabilized by conventional hydrogen bonds, electrostatic interactions, and hydrophobic contacts. Residues S104, P105, and M106 in LRR2 formed hydrogen bonds, while F108 formed a pi-pi T-shaped interaction with 4C/4C∗ of the ssDNA across all states. In addition, R74 and H76 in LRR1 formed electrostatic interactions with 3G/3G∗ and 4C/4C∗ of the ssDNA, and W47 (N-terminal LRR) together with W96 (LRR2) established pi-pi stacking interactions with 6 T/6 T∗ and 7 T/7 T∗ of the ssDNA. K181 (LRR5), K207 (LRR6), and K292 (LRR9) interacted electrostatically with 10 T/10 T∗, whereas Y208 and E616 (LRR20) formed hydrogen bonds with 10 T of the ssDNA. In clusters 1, 4, and 5, K181 and K207 also stabilized 9 T of the ssDNA. In the apo form, K181 (24% H-bond occupancy)/Y180 (18% H-bond occupancy) interacted with E616∗ but upon ligand binding, K181 shifted to form hydrogen bond with the CpG ssDNA while Y180-E616 H-bond occupancy increased to 70%.

Figure 8: CpG DNA and IM binding site residues of TLR9The residues considered within 4 Å of CpG DNA complex (A), and CpG and IM (B) binding site of TLR9. The single representative low energy complex was considered for clarity from FEL, hence, some of the residues discussed in the main text may not be present in the figure.

Surprisingly, K690 (LRR23) of TLR9 engaged the CpG ssDNA ligand bound to TLR9∗ electrostatically. Whereas H612 (LRR20), H641, and T642 (LRR21) formed hydrogen bonds with 5G, 6T, 8T, and 9T of the CpG ssDNA bound to TLR9∗. These findings suggest a redistribution of inter-monomer interactions upon ligand binding, reflecting structural adaptation within the ectodomain. We also observed from the results that across all states, ssDNA was primarily bound to one monomer but stabilized by interactions from both.

For the TLR9_CpG_IM complex, we analyzed metastable states 1, 2, and 3 (Figure 8B). All CpG-TLR9 interactions observed in the TLR9_CpG complex were retained, and residues H612, H641, and K690 continued to interact with CpG bound to TLR9∗. Thus, CpG ssDNA remained stabilized by both monomers even after IM binding.

In the modeled complex, the IM adopted an orientation in which its 5′ end inserted into the TLR9 dimer while the 3′ end extended outward. RMSD analysis supported the stability of this configuration (Figure 2B). IM actively interacted with R348, S350, and R426 of TLR9 and interestingly, residues Y536, G565, and N567 of TLR9 contacted with the C and T nucleotides of IM bound to TLR9∗. These shared interactions between CpG, IM, and both monomers support cooperative binding and are consistent with the 2:2 stoichiometry of the TLR9_CpG complex and the 2:2:2 stoichiometry of the TLR9_CpG_IM complex observed previously.24 In TLR9_apo, residues R348-D534∗/Y536∗/G565∗/H566∗, F402-M561∗, and R426-S509∗/D534∗ of TLR9 formed pi-cation interactions, electrostatic interactions, or hydrogen bonds, but onlyR348-D534∗/Y536∗/G565∗/H566∗ and R426-D534∗ were observed in CpG bound state. Upon IM binding, however, R348, F402, and R426 re-engaged through pi-pi, electrostatic, and hydrogen-bond interactions with the 2C nucleotide of IM, consistent with interactions previously reported in mouse TLR9.24 These results suggest that ligand binding induces a coordinated redistribution of inter-monomer interactions and uncovers previously unrecognized functional residues (H612, H641, T642, and K690) that participate in stabilizing CpG and IM ligands. Based on our simulations, we hypothesize that CpG binding induces a conformational rearrangement around the Z loop that enables IM engagement, and that the orientation of IM further promotes favorable interactions between the 5′-xCx motif and the TLR9 dimer.

Residue betweenness centrality analysis of apo- and CpG ssDNA and IM-bound TLR9 dimer complexes

To determine the critical residues in TLR9 and its ligand-bound complexes, we analyzed the betweenness centrality (CB) of these residues in the protein, which reflects the contribution of each node to intra-protein signal flow (Figure S4). Subsequently, to quantify changes in signal propagation between apo and ligand-bound states, we calculated the absolute differences |CB (TLR9_apo) _ CB (TLR9_CpG) |, |CB (TLR9_apo) _ CB (TLR9_CpG_IM) |, |CB (TLR9_CpG) _ CB (TLR9_CpG_IM) | using a cutoff of ≥0.3 to identify residues with the largest shifts. These residues represent key allosteric nodes and are shown in Figure 9.

Figure 9: Residue CB variations calculated between apo and ligand bound states of TLR9The plot shows the difference in CB among the three complexes with cutoff ≥0.3 (absolute values) for both monomers of TLR9. Locations of the residues observed ≥0.3 are mapped using dot representation on the structures and also listed on the tables.

We examined the residues shared by the TLR9 and TLR9∗ monomers, and also reported residues in the dimer interface shared by TLR9∗ as they showed high CB variation. For |CB(TLR9_apo) – CB(TLR9_CpG)|, residues P105 (LRR2), H260 (LRR8), T642 (LRR21) and L665 (LRR22). Y180∗ (LRR5), K532∗ (LRR17), N555∗ (LRR18), S589∗ (LRR19), T642∗ (LRR21), L665∗ (LRR22), P785∗ (C-terminal LRR), and S786∗ (C-terminal LRR) exceeded the cutoff in TLR9∗ monomer. Among these, Y180 and T642 lie at the dimer interface and mediate inter-monomer interactions as seen in interaction analysis.

In contrast, for |CB(TLR9_apo) – CB(TLR9_CpG_IM)|, residues N-terminal LRR E36, D424 (LRR14), C507 (LRR16), M561 (LRR18), I587 (LRR19), N608 (LRR20) and G783(C-terminal LRR) were common residues in both monomers with high CB values. C507∗ (LRR16), F559∗ (LRR18), M561∗ (LRR18), T642∗ (LRR21), and K738∗ (LRR25) showed the largest CB differences in TLR9∗. Similar to the CpG-bound comparison, these residues cluster near the dimer interface and the Z loop, suggesting a role in stabilizing ligand-dependent conformational rearrangements. Although D424∗ (LRR14), C507∗ (LRR16), and F559∗ (LRR18) displayed high CB values, they did not form direct interactions; however, their side chains project toward the interface region. Additional residues: I587, N608, L610, and W614 (LRR19–20), also exceeded the cutoff but are located away from the interface and do not orient toward the opposing monomer.

Comparison of |CB(TLR9_CpG) – CB(TLR9_CpG_IM)| revealed further shifts in the Z loop and LRRs 15–21, indicating altered communication pathways upon IM binding. All residues with CB ≥ 0.3 for TLR9 monomers are listed in Figure 9. Notably, residues C507 and T642 have not been previously implicated in TLR9 function. In our analysis, however, both exhibit substantial increases in CB upon ligand binding, suggesting that they may serve as previously unrecognized allosteric mediators linking the CpG and IM binding sites. Their pronounced network-level response indicates a potential role in long-range signal propagation within the receptor. Together, these shifts in CB reveal a coordinated network of allosteric nodes spanning the dimer interface and Z-loop region, defining the pathways through which ligand binding reshapes signal propagation in TLR9.