Insights into the species-specific TLR4 signaling mechanism in response to Rhodobacter sphaeroides lipid A detection.

TLR4 in complex with MD2 senses the presence of lipid A (LA) and initiates a signaling cascade that curb the infection. This complex is evolutionarily conserved and can initiate the immune system in response to a variety of LAs. In this study, molecular dynamics simulation (25 ns) was performed to elucidate the differential behavior of TLR4/MD2 complex in response to Rhodobacter sphaeroides lipid A (RsLA). Penta-acyl chain-containing RsLA is at the verge of agonist (6 acyl-chains) and antagonist (4 acyl-chains) structure, and activates the TLR4 pathway in horses and hamsters, while inhibiting in humans and murine. In the time-evolved coordinates, the promising factors that dictated the differential response included the local and global mobility pattern of complexes, solvent-accessible surface area of ligand, and surface charge distributions of TLR4 and MD2. We showed that the GlcN1-GlcN2 backbone acquires agonist (3FXI)-like configurations in horses and hamsters, while acquiring antagonist (2E59)-like configurations in humans and murine systems. Moreover, analysis of F126 behavior in the MD2 F126 loop (amino acids 123-129) and loop EF (81-89) suggested that certain sequence variations also contribute to species-specific response. This study underlines the TLR4 signaling mechanism and provides new therapeutic opportunities.

class="Gene">Toll-like receptor 4 (class="Chemical">n class="Gene">TLR4), a pathogen recognition receptor (PRR) family member, is principally involved in the sensing of lipid A (LA), a pathogen-associated molecular pattern (PAMP), and mounts an immune response against invading pathogens1. LA detection involves complex formation with TLR4/MD2. When stable, this complex triggers downstream mediators that converge on nuclear factor (NF)-κB, leading to inflammatory responses2. Regulated immune responses can abate bacterial threats; however, uncontrolled responses could lead to sepsis, a life-threatening condition that is difficult to treat3.

LA is composed of a class="Chemical">diglucosamine diphosphate head group, liclass="Chemical">nked iclass="Chemical">n a β(1–6) maclass="Chemical">nclass="Chemical">ner, aclass="Chemical">nd is appeclass="Chemical">nded with 4–8 acyl chaiclass="Chemical">ns depeclass="Chemical">ndiclass="Chemical">ng oclass="Chemical">n the species of origiclass="Chemical">n4. LA iclass="Chemical">nitiates the complex formatioclass="Chemical">n through dimerizatioclass="Chemical">n of [class="Chemical">n class="Gene">TLR4/myeloid differentiation factor (MD)2-LA]256, which has the ability to trigger or suppress the immune system depending on the number of acyl chains. LAs containing 4 to 6 chains exhibit as inhibitor to activator in a chain-number dependent manner respectively7891011. Rhodobacter sphaeroides lipid A (RsLA) is a 5-acyl chain-containing LA that activates the immune systems of horses and hamsters, and inhibits the immune systems of humans and murine1213. It also competes with Escherichia coli's lipopolysaccharide (LPS) for soluble CD14, lipid-binding protein (LBP) and MD2, prevents LPS-induced shock in mice, and also blocks the binding of LPS to cells141516.

The role of class="Gene">TLR4 iclass="Chemical">n ligaclass="Chemical">nd biclass="Chemical">ndiclass="Chemical">ng aclass="Chemical">nd discrimiclass="Chemical">natioclass="Chemical">n is class="Chemical">not fully clear. class="Chemical">n class="Gene">TLR4 is thought to play a secondary role in ligand recognition because the MD2 residues that interact with TLR41718 are located at the edge of the ligand-binding cavity19. Moreover, the TLR4/MD2 complex has higher affinity for LPS than MD2 alone20. The binding region of MD2 to TLR4 is near to N-terminus6; whereas, TLR4 single nucleotide polymorphisms (D299G and T399I)2122 and mutation studies in humans23 that attenuated TLR4 signaling are far away from primary or dimer interface. The plausible interpretation of this unexplained phenomenon could be disrupted cooperative binding to LPS or distorted ligand-induced conformational changes during signal transduction24.

Crystallographic studies revealed that class="Gene">MD2 is a cup-like hydrophobic structure that hosts the acyl part of LA, coclass="Chemical">nclass="Chemical">nected by various loops5619. These loops play crucial roles iclass="Chemical">n class="Chemical">n class="Gene">TLR4/MD2 activation. Specifically, the loop that joins the βG and βH strands harbors the crucial amino acid F126 that, when mutated, abolishes MD2 activation by a variety of TLR4 activators25. This flexible loop dynamically switches between activator and inhibitor modes, i.e. it is projected into the solvent and away from solvent in inhibitor and activator mode respectively619, possibly stabilizing contacts between TLR4* (the second TLR4 in a complex) and MD2 in the [TLR4/MD2-RsLA]2 complex526. Mutating F126A, however, does not influence ligand binding, but abolishes the signal transduction26.

In recent years, many studies have clarified various aspects of class="Gene">TLR4/class="Chemical">n class="Gene">MD2 signaling in response to LA binding; however, molecular simulation approaches can be used to explain the microscopic details of these signaling complexes27282930. To reveal the dynamic behavior and underlying mechanisms of the differential responses of the TLR4/MD2 to RsLA, we performed molecular dynamics simulations (MDS) of RsLA in complex with human, murine, horse, and hamster TLR4/MD2. Moreover, conformational amenability, dynamic motion, and a number of other physical features were focused on in the analysis. Finally, binding energies between LA and TLR4/MD2 were computed, providing a rationale for the existence of polar and apolar forces in the TLR4/MD2-LA complex. These findings can facilitate the design of therapeutic interventions for septic shock.

Results

Sequence analysis and homology modeling

class="Gene">TLR4 aclass="Chemical">nd class="Chemical">n class="Gene">MD2 (human, horse, murine, and hamster) sequences were aligned using CLUSTALW (Figure S1A and S1B). The human and horse proteins were closely related [TLR4 (81.73%) and MD2 (75.62%)], while murine and hamster proteins showed higher similarity [TLR4 (81.26%) and MD2 (83.12%)] (Table S1). However, RsLA acts as antagonist in human and murine systems and as an agonist in horses and hamsters, suggesting that something beyond TLR4 and MD2 sequence similarity accounts for these differential behaviors (Figure 1).

Figure 1

TLR4/MD2-RsLA complex.

Human TLR4/MD2-RsLA along with their van der Waal surface has been given as a representative model to outline the overall complex. TLR4 is shown in magenta, TLR4* in orange, MD2 in cyan while RsLA is shown in heteratomic color with C, O, N and P in yellow, red, blue and magenta respectively. Hydrogen has not been shown to simplify the structure. Only one dimerization interface is given. The loops are labeled in their respective color. (*) indicates the second TLR4 in a complex.

Whether the ligand-specific actclass="Disease">ivatioclass="Chemical">n of the complex is imparted by class="Chemical">n class="Gene">TLR431 or MD232 has yet to be established. Both proteins have been implicated as the principal modulator of ligand-based activation in different studies; however, contributions of both TLR4 and MD2 are imperative for activation. For instance, Walsh et al.23, reported that R384 residue in horse TLR4 is fundamentally important for complex activation, and its mutation to glycine completely abolishes lipid IVa-induced TLR4 activation. Notably, only horse TLR4 has an R residue at this position, the other three species having G or A instead (Table 1). Therefore, the influence of this amino acid might be limited to lipid IVa, while the response to RsLA might account other differences as well33. Moreover, N, L, D, and F residues found at TLR4 position 468 in humans, murine animals, horses, and hamsters, respectively, can potentially affect dimerization with MD2. TLR4 amino acids 290–350, which correspond to the intermediate region of leucine-rich repeat (LRR)8 and LRR9 including LRR9 as well, represent the most heterogeneous region, containing few conserved amino acids. This region is located immediately above of the region in which TLR4 interacts with the LA backbone and influences the behavior of the ligand (Figure S1A).

Table 1

Characteristic sequence specificity in different species. Human and murine behave as antagonist while horse and hamster behave as agonist toward RsLA. Amino acid positions marked with (*) and without asterisks belong to TLR4 and MD2, respectively. Only those residues that have proved influence or expected due to their positions are given. Human proteins were used as references for numbering the residues in both TLR4 and MD2

Species/Position	Human	Murine	Horse	Hamster
42	Y	F	S	F
57	S	T	T	T
61	L	V	L	L
64	F	E	L	K
65	Y	F	F	F
69	R	G	R	R
73	Q	Y	K	K
82	V	V	M	V
85	M	I	L	I
87	L	L	F	V
89	K	K	M	T
122	K	E	R	K
125	K	L	R	L
369*	E	K	E	E
384*	G	A	R	G
388*	K	S	K	R
397*	G	G	K	G
400*	S	S	R	S
436*	Q	R	Q	K
441*	S	S	P	S
468*	N	L	D	F

At class="Gene">MD2 positioclass="Chemical">n 89, iclass="Chemical">n the loop coclass="Chemical">nclass="Chemical">necticlass="Chemical">ng βE aclass="Chemical">nd βF (hereiclass="Chemical">n referred to as loop EF), class="Chemical">n class="Species">human and murine proteins contain a basic residue (K), while hydrophobic/polar residues are present in horse (M) and hamster (T) proteins. Loop EF is near the dimerizing interface, potentially influencing the complex stability (Figure 1). Moreover, Muroi et al.32, reported that few murine MD2 residues are important for the response to lipid IVa. Mutating these residues to corresponding human residues (T57S, V61L and E122K), impeded the response to lipid IVa. However, in the species under study, this did not yield any meaningful correlation, consequently eliminating any possible sequence-activity correlation. MD2 amino acids 57–66 and 82–8932, or 57–79 and 108–13523 are crucial for the response to lipid IVa (Table 1). Swapping these residues for corresponding residues from other mammalian species altered the response pattern. These blocks yielded heterogeneous sequence specificity precluding any sequence-structure relationship.

Structural flexibility and stability of TLR4/MD2 complexes

n class="Gene">TLR4 aclass="Chemical">nd class="Chemical">n class="Gene">MD2 complexes with docked RsLA were assembled in a hexameric form, and simulated for 25 ns. The human complex in the absence of RsLA and RsLA alone were simulated as controls. The global structural stability and compactness of the complexes were measured using the root-mean-square deviation (RMSD) of the backbone atoms (Figure 2A and B). The overlapped images obtained at 100 ps intervals showed that the secondary and tertiary structures were well maintained.

Figure 2

Root Mean Square Deviation (RMSD) plots of ligand bound and unbound complexes.

(a) RMSD plots of the backbone atoms in the TLR4/MD2 complex after least-square fitting to the initial conformation of the backbone atoms over a 25 ns simulation. The human TLR4/MD2 complex without ligand was simulated as the control. (b). RMSD plots of RsLA atoms bound to the TLR4/MD2 complex after least-square fitting to the initial conformation over a 25 ns simulation. RsLA alone was simulated in water as a control.

The F126 in class="Gene">MD2 is of critical importaclass="Chemical">nce iclass="Chemical">n propagaticlass="Chemical">ng the sigclass="Chemical">nal53435. Therefore, we decided to evaluate the dyclass="Chemical">namic behavior of the loop harboriclass="Chemical">ng this amiclass="Chemical">no acid as well as its coclass="Chemical">nformatioclass="Chemical">nal orieclass="Chemical">ntatioclass="Chemical">n. Mobility of the class="Chemical">n class="Gene">MD2 F126 loop was evaluated over the trajectory and was found to be more flexible in the control and human protein compared to horse. Similarly, murine was slightly more flexible than hamster (Figure 3A). Furthermore, conformational orientation the F126 side chain in human was not as stable as that of horse MD2, in which the side chain was rearranged from the outside into inside of protein and was stable when in contact with the acyl chain. The hamster and murine side chain was projected outward in the initial configuration and did not rearrange. In the hamster protein, the side chain was rearranged to an agonistic orientation to some extent, but because of the short time scale, the complete rearrangement phenomenon was not observed. We also examined the RMSD of the loop EF, which could play decisive role. The results showed that loop EF was stable in the horse and hamster proteins, while it had higher fluctuation in human, murine and control (Figure 3B). This plausibly explains the species-specific differences and the TLR4 signaling pathway.

Figure 3

Fluctuation of loops and solvent exposure of ligand.

(a) Root Mean Square Distance (RMSD) of atoms of the F126 loop over time. (b) RMSD of atoms of loop EF over 25 ns. (c) Solvent Accessible Surface Area (SASA) of whole RsLA compared to that of the protein complex over 25 ns.

The distance between the glycosidic class="Chemical">oxygen of class="Chemical">n class="Chemical">RsLA and the center of mass of the F126-loop was evaluated. Human had a higher distance and less stability than horse, while the murine had a greater distance and was less stable than that in the hamster protein. This distance decreased initially in horse and hamster proteins due to RsLA displacement, and then later stabilized, while the average distances in human and murine proteins were not constant (Figure S2). Interestingly, the distances remained stable in proteins with agonist behavior, but fluctuated in those with antagonist behavior, indicating that the stability of this loop is crucial to mount an immune response.

n class="Gene">MD2 is composed of a hydrophobic cavity where the class="Chemical">noclass="Chemical">n-polar portioclass="Chemical">n of class="Chemical">n class="Chemical">RsLA was buried, while the polar portion remained at the mouth of this cavity. Here, the hydrophobic interactions largely affected the outcome. Thus, the ratio of solvent-accessible to non-accessible area became a dominant factor in these protein-ligand interactions. The solvent-accessible surface area (SASA) was calculated, and to our intrigue, the SASA value was the highest and lowest for human and horse TLR4/MD2 complexes, respectively. Murine and hamster complexes had intermediate values, with the murine complex having a slightly higher value (Figure 3C).

Spatial orientation of RsLA in the MD2 hydrophobic core is also a factor that determines signaling behavior

Binding of the ligand to the receptor can reveal its possible actclass="Disease">ivatioclass="Chemical">n mode. LA coclass="Chemical">ntaiclass="Chemical">niclass="Chemical">ng 6 acyl chaiclass="Chemical">ns bouclass="Chemical">nd to class="Chemical">n class="Species">human TLR4/MD2 in an orientation where the 4′ PO4 group was pointing at the dimer interface (3FXI). This orientation was deemed the “normal” pose, while other ligand poses, were referred to as “flipped” poses (Table 2). RsLA orientation in human MD2 hydrophobic core was comparable to that of Eritoran (2Z65) and lipid IVa (2E59), indicating an antagonistic orientation. LA bound to murine MD2 in a similar orientation in all instances, regardless of whether LA contained 4 or 6 chains (3VQ1 and 3VQ2). Because the crystal structures for the horse and hamster complexes had not been resolved until now, we had to rely on docking poses. The ligand orientation in the horse MD2 was similar to that observed in the murine and to the previously reported pose36, while the orientation in the hamster core was similar as in the chicken (3MU3) and human proteins.

Table 2

Observed ligand binding patterns of TLR4/MD2 complexes and their resulting activities. Prefixes correspond are as follows: h (human), m (murine), e (horse), and ham (hamster)

Effecter activity	Complex	Orientation	Volume of ligand (Å3)	Reference
Agonist	hTLR4/MD2/LPS	Normal	2016	3FXI
Antagonist	hTLR4/MD2/RsLA	Flipped	1582	Modeled
Antagonist	hTLR4/MD2/Eritoran	Flipped	1412	2Z65
Antagonist	hTLR4/MD2/Lipid IVa	Flipped	1393	2E59
Agonist	mTLR4/MD2/LPS	Normal	2016	3VQ2
Agonist	mTLR4/MD2/Lipid IVa	Normal	1393	3VQ1
Antagonist	mTLR4/MD2/RsLA	Normal	1582	Modeled
Agonist	eTLR4/MD2/RsLA	Normal	1582	Modeled
Agonist	hamTLR4/MD2/RsLA	Flipped	1582	Modeled

The differential behavior of 5-acyl chains class="Chemical">RsLA is solely attributed to the structural features of the receptor iclass="Chemical">n the species examiclass="Chemical">ned iclass="Chemical">n this study37. Wheclass="Chemical">n this ligaclass="Chemical">nd was simulated iclass="Chemical">n the class="Chemical">n class="Gene">TLR4/MD2 complexes, behavior-dependent conformational changes were observed. For example, the displacement of RsLA was approximately 1.4 Å, 5.5 Å, 3.7 Å, and 6.4 Å into the human, horse, murine, and hamster MD2 cavities, respectively. The increased shift of RsLA in the horse and hamster cavities pushed their acyl chains, causing them to fold back onto themselves to a greater extent. Furthermore, the spatial orientation of GlcN1 → GlcN2 in RsLA disaccharide acquired a conformation resembling that of lipid IVa in human and murine complexes, while it was closely related to the confirmation of ReLA in horse and hamster complexes (Figure 4). This structural orientation has a strong influence on ligand behavior and signaling mechanisms38.

Figure 4

Representative configuration of GlcN rings.

A representative image of GlcN rings of RsLA bound to the (a) human (b) murine (c) horse, and (d) hamster TLR4/MD2 complex. 3D orientation of GlcN rings acquired the co-planer and twisted configuration as reported by crystal structures.

Moreover, the flexibility of dihedral angle around C1 → C6 linkage in class="Chemical">RsLA backboclass="Chemical">ne facilitates to rearraclass="Chemical">nge the acyl chaiclass="Chemical">ns iclass="Chemical">n class="Chemical">n class="Gene">MD2 such that MD2 acquired a favorable configuration for activation, without substantial entropic loss38. This dihedral was also shown to differ significantly in solution- and membrane-incorporated LA structures3940. These angles were found to be more stable in human complexes than horse. Similar patterns were observed in murine complexes, where the angles were more stable compared to those in hamster complexes (data not shown). In human and murine complexes, the dihedral angels in RsLA were stable around −120°, while in horse and hamster complexes, angels were adopted a conformation of around −150° or 150°, respectively.

RsLA agonizes and antagonizes downstream signaling from horse and human TLR4/MD2, respectively, despite their sequence similarity

To elucidate this differential behavior, the electrostatic surface potential of class="Gene">TLR4 aclass="Chemical">nd class="Chemical">n class="Gene">MD2 were evaluated and substantial differences among specific localized regions were found. The electrostatic potential of human TLR4 was strongly negative over the entire surface, especially at the primary interface with MD2. The portion of TLR4 that interacts with TLR4* was also populated with negative residues. The surface of horse TLR4 interacting with MD2 was negative, while the region between LRR8-15 was comparatively neutral. In addition, the dimer interface was less negative compared to that in the human complex (Figure 5).

Figure 5

Electrostatic surface potential of TLR4 ectodomains of different species.

TLR4 ectodomains are shown for human (a), murine (b), horse (c), and hamster (d) proteins. The electrostatic surface potentials of the TLR4 ectodomains (a–d, respectively) vary greatly among species. These electrostatic potentials were built on last snapshots of the 25 ns trajectory. PyMol was used to render the images, and potential values were set to −10 and 10 for negative and positive charges, respectively. Contour values were set to −5 and 5, respectively.

In n class="Species">human class="Chemical">n class="Gene">MD2, the dimerizing interface was more electrostatically positive and slightly positive in the core region, while the primary interface was predominantly negative (Figure 6). In horse MD2, the dimerizing interface was less positive than that of human MD2; the core was neutral while the primary interface was slightly more negative.

Figure 6

Electrostatic surface potential of MD2 of different species.

The surface potentials are shown for human (a), murine (b), horse (c), and hamster (d) MD2 proteins. These electrostatic potentials were built on last snapshots of the 25 ns trajectory. Potential values were set to −10 and 10 for negative and positive charges, respectively. Contour values were set to −5 and 5, respectively.

Structural flexibility of protein has been corclass="Gene">related to differeclass="Chemical">nt biological fuclass="Chemical">nctioclass="Chemical">ns iclass="Chemical">ncludiclass="Chemical">ng ligaclass="Chemical">nd recogclass="Chemical">nitioclass="Chemical">n aclass="Chemical">nd eclass="Chemical">nzymatic activity41. These structural flexibilities caclass="Chemical">n be studied by esseclass="Chemical">ntial dyclass="Chemical">namics or PCA42. This aclass="Chemical">nalysis has beeclass="Chemical">n desigclass="Chemical">ned to capture large-scale motioclass="Chemical">ns, where first few eigeclass="Chemical">nvectors caclass="Chemical">n accouclass="Chemical">nt for majority of proteiclass="Chemical">n motioclass="Chemical">ns. Similarly, the differeclass="Chemical">ntial respoclass="Chemical">nses of class="Chemical">n class="Species">human and horse receptors could be correlated to the localized fragmental motions when analyzed using PCA (Figure 7 and Figure S3). Human receptor was largely stayed in one conformational space indicating lower energy, while transition to any other space was sparse (Figure 7). In contrary, horse receptor visited two prominent conformational spaces separated by small energy barrier.

Figure 7

2D plots of first and second eigenvectors.

(a) Human and horse plots (b) murine and hamster plots have been shown in the 2D graph, where the X-axis represents first eigenvector and the Y-axis represents the second eigenvector.

class="Gene">MD2 movemeclass="Chemical">nt iclass="Chemical">n apoclass="Chemical">n class="Gene">TLR4/MD2 and ligand-bound TLR4/MD2 was surprisingly in opposite directions in human, as revealed by porcupine plots (Figure S3). In the apo complex, the movements of the TLR4s opposed each other at the C-terminal region as well as at the dimer interface. The N-terminus of one TLR4 was directed upward, while its counterpart was the least mobile component in the whole complex. Similarly, MD2 was either the least mobile or pointing away from the dimer interface. Ligand-bound human TLR4 had uniform motions, where the N- and C-termini were facing each other. The intermediate region was also pointing along the axis of TLR4. In MD2, mixed motions were observed, which accounting for as well as away from dimer with intermediate magnitude. Binding of the ligand may have influenced the essential motions of the surrounding atoms, facilitating dimer formation and increasing its stability compared to the unbound state. Furthermore, in horse TLR4, the N-terminal regions moved in opposite directions, while the C-terminal regions moved away from each other. The dimer-forming regions of both TLR4 proteins moved unidirectionally. Horse MD2 was largely moving towards the dimer interface.

In RMSF, the terminal regions were highly flexible in both class="Species">human aclass="Chemical">nd class="Chemical">n class="Species">horse TLR4, whereas the centre displayed restricted fluctuation except loop regions. MD2 behaved similarly where loops were also highly mobile as it was obvious for loop EF and loop Phe126 (Figure S4). These are in line to the observations as concluded by porcupine plots.

Binding free energy is valuable to consider in complex stability. There are various energy terms that contributes to the complex formation, broadly categorizing into polar and non-polar energies. The protein-ligand binding energies of class="Species">human aclass="Chemical">nd class="Chemical">n class="Species">horse receptors surprisingly varied in both cases. Various non-polar energies terms i.e. van der Waals (VdW), SASA and solvent accessible volume (SAV) energies—favored TLR4/MD2-based downstream signaling, while polar salvation energy and electrostatic energy hindered this behavior (Table 3 and S5). These energy terms delineate that non-polar energy terms are predominant in human and horse receptors.

Table 3

MM-PBSA values (binding free energies) of RsLA-protein complexes. Each value represents the average value calculated for the last 10 ns of each trajectory with 10 snapshots. SASA, solvent accessible surface area; SAV, solvent accessible volume; WCA, Weeks-Chandler-Andersen

	Energies (kJ/mol)	Human	Murine	Horse	Hamster
Polar Energies	Electrostatic	−3798.51 +/− 244.91	1360.2 +/− 129.72	282.43 +/− 48.08	245.63 +/− 32.63
	Polar Solvation	647.50 +/− 89.15	332.89 +/− 19.71	272.13 +/− 28.87	350.46 +/− 14.38
Non-Polar energies	Van der Waal	−873.54 +/− 21.90	−970.53 +/− 30.82	−1069.11 +/− 29.5	−1012.96 +/− 32.0
	SASA	−108.36 +/− 2.25	−115.83 +/− 2.55	−121.33 +/− 2.62	−117.65 +/− 2.51
	SAV	−1188.88 +/− 36.96	−1206.60 +/− 42.9	−1289.9 +/− 48.78	−1254.09 +/− 23.0
	WCA	357.31 +/− 113.07	573.02 +/− 6.83	600.97 +/− 3.99	570.31 +/− 11.94
	Binding Energy	−4964.5 +/− 235.53	−26.85 +/− 146.03	−1324.80 +/− 74.5	−1218.31 +/− 43.1

RsLA triggers downstream signaling through hamster receptors, while hinders signaling through murine receptors despite their high sequence similarity

The intriguing behavior of higher sequence similarity and different responses could be explained by multiple ways. class="Species">Murine class="Chemical">n class="Gene">TLR4 was negatively charged at the N-terminal, while the central region was slightly positive and the C-terminal is negative to neutral. Compared to the other receptors, hamster receptors were slightly unusual in that the N-terminal of TLR4 was negative, the intermediate region was slightly neutral, and the C-terminal was negative (Figure 5). The dimerizing face of murine MD2 was minimally positive, the core was neutral, and primary interface was positive. In hamster MD2, the dimer interface was substantially positive along with the core, while a patch of negative residues was observed at the primary interface (Figure 6).

The essential dynamics revealed that different patterns of motion occurred in the complexes in these species. class="Species">Murine dyclass="Chemical">namics data showed a similar behavior to that of class="Chemical">n class="Species">human, where murine receptor visited multiple conformational spaces almost equally with least energy distinctions or prominent sampling. In contrast, hamster could be sampled in two prominent conformational spaces with clear energy barriers as observed in horse (Figure 7). This dynamics, although slightly different from that of the horse, might be the reason that hamster receptors behave agonistically.

Porcupine motions indicated a highly mobile nature of class="Species">murine class="Chemical">n class="Gene">TLR4. All parts of TLR4 showed coordinated movements in opposing directions (Figure S3). MD2 largely pointed in the upward direction, tilting slightly toward the dimerizing interface, potentially accounting for its antagonistic behavior. In hamster, the motion of the TLR4 N-terminus was directed opposite to that of the C-terminus, while the dimer region displayed mixed motion directions. In MD2, mixed directional motions were observed, accounting for both dimerizing as well as away from dimer interface.

class="Species">Murine aclass="Chemical">nd class="Chemical">n class="Species">hamster has the highest fluctuation as seen by RMSF graphs, both in TLR4 and MD2. In TLR4, murine had the higher fluctuation followed by horse, while in MD2, murine was the one displaying higher fluctuations irrespective of the chain. In case of hamster, the fluctuations were of intermediate intensity in MD2, whereas these were of higher in TLR4 (Figure S4).

The binding energies also varied for various sub-components. The non-polar terms were also dominant in this case. For example, VdW, class="Chemical">SASA, aclass="Chemical">nd class="Chemical">n class="Chemical">SAV energies were lower in hamster complexes, making them more stable than murine complexes (Table 3). In contrast, polar salvation energy was the only factor that opposed complex stability. Together, these factors strongly influence the outcomes of protein-ligand interactions during complex formation. These data suggested that non-polar terms markedly affected complex stability (Figure S5).

Alanine scanning mutagenesis

Hot spots residues at the protein-ligand interfaces in all four complexes of class="Gene">LPS with class="Chemical">n class="Species">human, mouse, horse and hamster MD2 were identified through DrugScorePPI server43. DrugScorePPI is a scoring function based on the knowledge of computational alanine scanning in protein-protein and protein-ligand interfaces. The binding free energy differences between the wild-type residues in hot spots and the Alanine mutants in respective chain in complex were calculated with the following formula.

A higher positive ΔΔG (kcal/mol) value indicates a hot spot with strong binding affinity and vice versa. Upon analysis, we found that the hydrophobic pocket of class="Gene">MD2 iclass="Chemical">n all four species is completely filled by class="Chemical">n class="Chemical">RsLA with slight difference in their binding pattern and residues involved in binding. MD2-RsLA complexes were held together primarily by hydrophobic interaction between the hydrophobic tails of RsLA and binding cavity of MD2. The important residues were labeled. The conserved residues are all having higher energy penalty such as F76, I80, I94, and F121. Particularly, Y102, which is one of the key hot spot showing higher energy perturbation, may account for antagonistic characteristics of RsLA in human and murine by stabilizing these complexes in one conformation rather to explore other conformations (Figure 8).

Figure 8

Alanine scanning mutagenesis and Hotspot identification.

Alanine scanning was performed through DrugScorePPI for MD2-LPS of human (a), murine (b), horse (c), and Hamster (d). The binding free energy of key hotspots has been given along Y-axis, residues having >600 calc/mol are labeled and considered as most important hot spots. The positive ΔΔG values represent the potential hot spot residues contribution in the ligand-receptor complexes.

Discussion

In this study, we described the class="Chemical">RsLA iclass="Chemical">nduced class="Chemical">n class="Gene">TLR4/MD2 signaling mechanism in human, murine, horse and hamster. Based on these data, various factors have been identified that must be considered when designing ligands to manipulate these receptors. Importantly, when considering the species-specific differences in these cases, a few factors, such as protein flexibility and non-polar character, play crucial roles than the conserved regions of TLR4 and MD2.

Sequence similarities among these species were strikingly different and deviated from the assumption that proteins with greater sequence similarity behave similarly. Therefore, we can assume that more than just sequence similarity accounts for the differences in behavior. Moreover, specific foods and habitat patterns might evolve the class="Species">murine immuclass="Chemical">ne system that specifically recogclass="Chemical">nizes eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">ntal bacteria36, but class="Chemical">n class="Species">hamster response lacks any comprehensible explanation. Similar explanation can be implied to the observed TLR4 residue R384 (Figure S1 and Table 1). The sequence substitutions may influence the local conformation alterations in 3D structure of TLR4 as observed by X-ray crystallography2344. The one residue change (T399I), though it is not present on or near active site, can hamper the signaling by altering the local structure, this behavior can be implicated to account the influence of other residues substitution (Table 1). So far, the sequence-structure and phylogenetic relationship are the most plausible and strongly correlated factor in these analyses4546. In our study, we selected two species each for the agonistic (horse and hamster) and antagonistic (human and murine) groups to strengthen the in silico analysis23. Furthermore, we evaluated the polar/non-polar characteristics of these complexes based on the rules described by Kyte and Doolittle47. When whole ectodomain sequences were evaluated, no obvious relationship was observed between their behavior and hydrophobic character (data not shown). When only non-conserved residues were analyzed, horse and hamster TLR4 receptors were less hydrophobic than the TLR4 proteins of the antagonistic group. Similarly, when the relative hydrophobicity of RsLA to protein was analyzed over the trajectory, horse and hamster were lower than human and murine (Figure 3C). The electrostatic surfaces of TLR4 and MD2 of horse and hamster are more neutral than human and murine (Figures 5 and 6). Because of these differences in surface potential, it can be argued that protein-ligand binding energies vary significantly. Electrostatic and solvation free energies are markedly different in these complexes, while the non-polar energy terms are in line with the experimental observations (Table 3). The lower non-polar terms signified the formation of stable complexes and highlighted the importance of non-polar energies in TLR4/MD2 signaling. Cumulatively, the hydrophobicity of the receptor molecules could be correlated to the signaling pattern.

Most studies have compared class="Species">human aclass="Chemical">nd class="Chemical">n class="Species">murine receptors with lipid IVa that possesses a characteristics agonist and antagonist behaviors. Therefore, in these cases, mutation of relevant amino acids would profoundly affect the outcome. In the current study, 5-acyl chain LA may not be able to impart much of an effect if similar amino acids were replaced. Moreover, this ligand behaves similarly in human and murine receptors, which may require complementary explanations alongside mutation-based conclusions. Recently, a double mutation of R384G/P441S completely abolished the horse TLR4 activation33. Intriguingly, human and hamster possess the identical residues at these positions, which suggests that these mutations may be limited to horse, and to explain hamster response other factors may also be needed. Moreover, a few caveats should be borne in mind, such as, horse and hamster proteins were modeled using human and mouse proteins as templates. This would generate structures overly similar to their templates, which might bias the results and preclude the vivid structural changes. Moreover, human and horse can be justified with similar explanations, whereas, being rodents, murine and hamster lack any similar justifications.

n class="Chemical">Alanine scaclass="Chemical">nclass="Chemical">niclass="Chemical">ng aclass="Chemical">nalysis also revealed hot spots that are mostly coclass="Chemical">nserved amoclass="Chemical">ng the sequeclass="Chemical">nces (Figure 8 aclass="Chemical">nd Figure S1). While coclass="Chemical">nsidericlass="Chemical">ng their side chaiclass="Chemical">n property, almost all are of hydrophobic iclass="Chemical">n class="Chemical">nature with aclass="Chemical">n exceptioclass="Chemical">n as iclass="Chemical">n case of class="Chemical">n class="Species">human, R90, E92, D101 are polar. This analysis also confirmed the role of hydrophobic residues in TLR4 signaling.

Proteins are inherently flexible that, upon binding to its cognate molecule, show transition from one ensemble to energetically favorable ensemble to perform its functions. These local or global transitions are vital to understand the complex formational modulations in protein receptors. Moreover, such intrinsic dynamics of biological macromolecules such as class="Gene">TLR4/class="Chemical">n class="Gene">MD2 complex is the crucial feature in structure-activity relationships4849. PCA was performed to elucidate such intrinsic flexibility. In PCA analysis, 2D plots of the first two eigenvectors revealed a conformationally restrained space for human and murine receptors that were either stuck in one space or found that conformational ensemble energetically favorable with RsLA. Intriguingly, horse and hamster proteins visited two distinct conformational spaces that indicate their favorable energy ensembles and plausibly explain their deviated response (Figure 7). The observing of two distinct energetically favorable conformations in horse and hamster, while human and murine are lacking such behavior, might be correlated to their observed responses. Similarly, the movements in ligand-bound and unbound MD2 complexes implied that this protein had inherent mobility that might be independent of ligand (Figure S3). This flexible behavior of MD2 favored the binding of endotoxins with different characteristics and would be explored in further studies.

In conclusion, the host response to class="Chemical">RsLA could be attributed to local aclass="Chemical">nd global fluctuatioclass="Chemical">n aclass="Chemical">nd class="Chemical">noclass="Chemical">n-polar eclass="Chemical">nergy coclass="Chemical">ntributioclass="Chemical">ns to the dimerizatioclass="Chemical">n. These results revealed uclass="Chemical">ndetermiclass="Chemical">ned mechaclass="Chemical">nisms of proteiclass="Chemical">n iclass="Chemical">nteractioclass="Chemical">n aclass="Chemical">nd provided a plausible explaclass="Chemical">natioclass="Chemical">n for species-specific sigclass="Chemical">naliclass="Chemical">ng respoclass="Chemical">nses. Thus, class="Chemical">n class="Disease">MDS is a valuable tool for designing and testing new drugs. The timescale used in this study was not sufficient to decipher complete structural details; however, these data will contribute to existing knowledge and assist scientists in the exploration of other unidentified factors that are involved in TLR4 signaling, and guide the development of new molecules for pharmacological manipulation of TLRs.

Methods

Homology modeling and ligand docking

class="Species">Human aclass="Chemical">nd class="Chemical">n class="Species">murine X-ray crystal structures of TLR4/MD2 were retrieved from the Protein Data Bank (3FXI and 2Z64) and were used as templates to construct homology models for horse and hamster TLR4/MD2 by using MODELLER v.9.1150 respectively. The generated models were selected owing to their maximum score and minimum violation. These models were then evaluated by using ProSA51, Verify3D52, and Ramachandran Plot53, and were found to be in an acceptable range. For further refinement, all fragments were simulated independently, improving the quality of these models, and the last frames were used as the starting structures for further experiments.

The class="Chemical">RsLA structure was derived from the LA structure of class="Chemical">n class="Species">E. coli, and was modified using ChemBioDraw 13.0 (Trial Version). For ligand topology, hydrogens were first added and then the partial charges were calculated according to the AM1-BCC model54, and using the online version of ACPYPE55 topology was created. The created topology was edited in order to be compatible with the AMBER force field in GROMACS. Protein-ligand docking was performed using AUTODOCK VINA56, which is superior in finding the correct ligand pose to AutoDock 4.0. The implemented algorithm in AUTODOCK VINA (Iterated Local Search global optimizer) follows an iteration process involving mutations and local optimization allowing the acceptance of each iteration based on Metropolis criterion. For receptor, TLR4/MD2 heterodimer that was kept rigid, was docked against RsLA, which was partially flexible with the maximum allowed (i.e. 32) torsion angles during docking. A minimum grid of 30 × 30 × 28 with a spacing of 1.0 Å was selected around the MD2 hydrophobic core as the docking site. Unless otherwise stated, parameters were kept at their defaults. Moreover, increasing the exhaustiveness did not yield any better results. Multiple docking cycles were carried out and a suitable ligand-bound pose was selected based on least energy and conformational acceptability in each complex. The final complexes with ligand-bound conformations were generated based on the available crystal structures of LA (3FXI)5, lipid IVa (2E59, 3VQ1, and 3VQ2)24, and Eritoran (12) (2Z65)6, by performing pair-wise structural alignments of the TLR4/MD2 conformations with PyMol (www.pymol.org). The control system was established by removing the ligand from the TLR4/MD2 complex crystal structure (3FXI).

Molecular dynamics simulation protocol

All simulations were performed in GROMACS v.4.6.257 with the AMBER99SB-ILDclass="Chemical">N58 parameter sets aclass="Chemical">nd the TIP3P class="Chemical">n class="Chemical">water model59. For long-range electrostatics, the particle mesh Ewald method was employed60, using a 10 Å cut-off distance for real-space Ewald interactions and VdW interactions. To accommodate energy and pressure terms due to VdW truncation, dispersion correction was applied. To mimic the infinite system, periodic boundary conditions were applied in all directions. LINC algorithm61 was used to constrain bond lengths in all atoms. All production simulations were performed with a 2 fs time step under 1 bar of pressure and 300 K, without restraining any position.

The complexes were placed into cubic boxes filled with TIP3P explicit class="Chemical">water represeclass="Chemical">ntatioclass="Chemical">n aclass="Chemical">nd class="Chemical">neutralized with couclass="Chemical">nter ioclass="Chemical">ns before addiclass="Chemical">ng 100 mM class="Chemical">n class="Chemical">NaCl, in total ~400,000 atoms in each box. Steepest descents minimization was performed to remove any unfavorable interactions, and then two-step equilibrium was used to generate the starting structures for the production simulations. During equilibration, position restraints were applied to all atoms to avoid any configuration changes. In the first phase, system was simulated for 100 ps under a constant volume (NVT) ensemble to achieve 300 K by V-rescal method62. For temperature-coupling, the protein and ligand was treated as a single group and ions and water were considered as a second group. The equilibrated structures from NVT ensemble were subjected to 100 ps of constant pressure (NPT) equilibration, under an isotropic pressure of 1.0 bar, using the Parrinello-Rahman barostat63. Production MDSs were performed for 25 ns in the absence of any restraints. During this data collection period, the V-rescale thermostat and Parrinello-Rahman barostat were used to maintain temperature and pressure at 300 K and 1 bar respectively. These settings have been shown to yield a true NPT ensemble636465.

To evaluate the vital role played by the hotspot residues of class="Gene">MD2 iclass="Chemical">nteracticlass="Chemical">ng with class="Chemical">n class="Gene">LPS, the last snapshot of each complex was subjected to Alanine scanning mutagenesis using DrugScorePPI web server43. This server, using knowledge based function, search for interface residues, calculates ΔGWT (wild type) for all interface residues, mutates them to Alanine sequentially and then calculates the ΔGMUT (mutant type), which allows subsequent calculation of ΔΔG (change in binding free energy). This procedure iterates until ΔΔG for all interface residues have been calculated. The calculated ΔΔG points to the hot spot residues.

Data analysis

GROMACS was used for trajectory analyses, and PyMol (www.pymol.org) and UCSF Chimera66 were used to create images. All electrostatic surfaces were calculated using the online class="Gene">PDB2PQR server (http://class="Chemical">nbcr-222.ucsd.edu/class="Chemical">n class="Gene">pdb2pqr_1.8/)67 and images were rendered using PyMol. Multiple sequence alignments were performed using CLUSTALW and further analyses were performed using BioEdit. Numbers were assigned based on human TLR4 and MD2 as reference. For Principal component analysis (PCA)42, eigenvectors of protein backbone have been calculated, while, coordinates of every 900 ps were used for porcupine plots. Binding energies were calculated for the last 10 ns by using the MM-PBSA tool as described by Kumari et al.68, utilizing APBS69 to solve Possion-Boltzman equations. The extensive details have been described in the paper and on the web (http://rashmikumari.github.io/g_mmpbsa/ accessed 21040920). For a brief overview, periodicity was removed in each complex before calculating the binding energies. The last 10 ns of each trajectory saved at 2 ps frame subjected for the binding free energies that have been calculated for electrostatic, polar salvation, SASA, SAV, and WCA energy terms as given in Table 3. The default parameters were used in all instances.

Author Contributions

M.A.A. and S.C. planned experiments. M.A.A. and S.P. performed experiments. M.A.A., S.P. and M.S. analyzed data. S.C. contributed for material. M.A.A., S.P., M.S. and S.C. wrote the paper.