T cells use their T-cell receptors (TCRs) to scan other cells for antigenic peptides presented by MHC molecules (pMHC). If a TCR encounters a pMHC, it can trigger a signalling pathway that could lead to the activation of the T cell and the initiation of an immune response. It is currently not clear how the binding of pMHC to the TCR initiates signalling within the T cell. One hypothesis is that conformational changes in the TCR lead to further downstream signalling. Here we investigate four different TCRs in their free state as well as in their pMHC bound state using large scale molecular simulations totalling 26 000 ns. We find that the dynamical features within TCRs differ significantly between unbound TCR and TCR/pMHC simulations. However, apart from expected results such as reduced solvent accessibility and flexibility of the interface residues, these features are not conserved among different TCR types. The presence of a pMHC alone is not sufficient to cause cross-TCR-conserved dynamical features within a TCR. Our results argue against models of TCR triggering involving conserved allosteric conformational changes.
The interaction between T-cell receptors (TCRs) on the surface of T-cells and peptides bound by Major Histocompatibility Complexes (MHCs) on the surface of antigen presenting cells is one of the most important processes of the adaptive immune system [1]. In the case of MHC class I molecules intracellular proteins are degraded by proteasomes into peptides, the peptides are loaded onto MHCs, and subsequently the peptide/MHC (pMHC) structures are presented on the cell surface. The TCRs of T-cells bind to pMHCs with their six hypervariable Complementarity Determining Regions (CDRs) and thereby scan the pMHC (Fig 1A and 1B). Based on this interaction further downstream signalling cascades can be activated and an immune response can be elicited against a particular antigenic peptide. The TCR/pMHC interaction is of relatively low affinity (KD ~0. 1–500 μM) and degenerate: One TCR can recognise multiple pMHC and one pMHC can be recognised by multiple TCRs but not every TCR can recognise every pMHC. A long standing question has been how TCR binding to pMHC results in changes in the cytoplasmic domains of the TCR/CD3 signalling subunits (e. g. phosphorylation), a process termed TCR triggering. Several mechanisms of TCR triggering have been proposed [2], which can be grouped roughly into segregation/redistribution, aggregation, and conformational change models. Conformational change models are of two types [2]: one group, based on the observation that pMHC binding will impose a mechanical pulling force on the TCR [3], proposes that this mechanical force somehow alters the conformation of TCR relative to the membrane and/or CD3 subunits [4–6]. The second group postulates that pMHC binding is accompanied by an allosteric conformation change transduced through the TCRαβ heterodimer [7]. While there is growing support for mechanical models [8–11], evidence for ‘allosteric’ models is more equivocal (reviewed in [2]). Rossjohn and colleagues have reported a change in the AB loop of TCR Cα domain accompanying binding to agonist pMHC [7,12], but structural studies have failed to identify a conformation change conserved in all TCRs upon pMHC binding (reviewed in [2]). The experimental studies described above relied on X-ray crystallography or NMR spectroscopy, which cannot measure dynamic changes in the TCR structure at atomistic resolution, leaving open the possibly that pMHC binding results in conserved changes in TCR dynamics. Fortunately, computational methods such as Molecular Dynamics (MD) or Monte Carlo simulations can be used to explore this possibility (reviewed in [13]). Recent approaches include TCR/pMHC interface H-bond network analysis [14,15], binding free energy [16] and detachment [17,18] simulations, effects on the MHC [19], peptide [20], and on the TCR [21] or CDR loop characterisation [22]. In a previous large scale study we analysed the same TCR/MHC in combination with 172 different peptides of known experimental immunogenicity [23]. The aim was to use simulations to seek distinct dynamical TCR behaviours for encounters with immunogenic versus non-immunogenic peptides. While no such differences were identified to be involved in discrimination, this does not rule out a role for conformational dynamics in triggering itself. In the current study we explore this by simulating four different TCRs in their agonist pMHC bound and unbound state and then look for conserved conformational features that distinguish bound from unbound TCRs. Using a total of 26 000 ns of simulation time we show that no such conserved and non-obvious features distinguish bound and unbound TCR.
We extracted the structures of the LC13-TCR/FLRGRAYGL/HLA-B*08: 01 (accession code 1MI5), A6-TCR/LLFGYPVYV/HLA-A*02: 01 (accession code 1AO7), JM22-TCR/GILGFVFTL/HLA-A*02: 01 (accession code 1OGA), and 1G4-TCR/SLLMWITQC/HLA-A*02: 01 (accession code 2BNR) from the Protein Data Bank (PDB) [24]. Constant TCR regions were included in all simulations as we have previously shown that the constant regions are important for reliable conclusions from molecular TCR and antibody simulations [25]. The LC13 was chosen as a model system as this system has been the target of extensive MD simulations (e. g. [23,26]) due to the availability of a large number of experimentally tested mutations. Additionally we chose the A6, JM22, and 1G4 systems as they have been investigated computationally before (e. g. [14,15,18,27]) and all three TCRs bind to HLA-A*02: 01 while the LC13 TCR binds HLA-B*08: 01. This allows us to investigate if there are conserved TCR reaction features across MHC types as well as within HLA-A*02: 01. Eight different structures (LC13-TCR, A6-TCR, JM22-TCR, and 1G4-TCR with and without pMHC) were simulated. Each structure was immerged into a dodecahedronic simulation box filled with explicit SPC water allowing for a minimum distance of 1. 2 nm between protein and box boundary. Na+ and Cl- ions were added to achieve a neutral charge and a salt concentration of 0. 15 mol/litre. Protonation states were determined automatically by Gromacs [28]. Energy minimisation using the steepest descent method was applied. The systems were warmed up to 310 K using position restraints. Hydrogen atoms were replaced by virtual sites to allow for an integration step of 5 fs for the production runs [29]. Final production runs were carried out using Gromacs 4 [28] and the GROMOS 53a6 force field [30]. Parts of the LC13 simulations were taken from our previous work in [21] and [23] while parts of the A6, JM22, and 1G4 simulations were taken from [14] and [15]. Multiple replicas (identical parameter but different seeds) per simulation are important for reproducible conclusions as the results of several studies [15–17,25,31] and a systematic evaluation have shown [32]. Therefore we use 100 LC13 TCR simulation replicas at 100ns each. A boot strapping analysis of these 100 replicas (S1 Fig) shows that 10 replicas reduce the variance between replicas by about 70% while more replicas reduce the variance only slowly further (25 replicas: 80% and 50 replicas 87%). Therefore we simulated the A6, JM22, and 1G4 TCR with and without pMHC for 100 ns using 10 replicas totalling 26 μs (Table 1). Trajectories were manually inspected using VMD [33] and the vmdICE-plugin [34]. Solvent accessible surface area (SASA), root mean square fluctuation (RMSF), radius of gyration (RG), hydrogen bonds (H-bonds), and distances were calculated by the GROMACS [28] modules gmx sasa, gmx rmsf, gmx gyrate, gmx hbond, and gmx distance respectively and imported into pymol/Matlab using gro2mat [35]. The H-bond networks were visualized using pyHVis3D [14]. We used three different types of measurements to quantify the magnitude of difference between descriptors of TCRpMHC and TCR simulations (modified from our previous study [19]): Firstly, the simple difference in the mean values that is referred to as: d=X¯TCRpMHC−X¯TCR Where X-TCRpMHC and X-TCR are the mean values over all frames and replicas of descriptor X (e. g. SASA or H-bonds of a region). The value d helps to quantify the actual magnitude of difference e. g. TCRpMHC simulations have on average 0. 5 H-bonds less between their TCR chains than TCR simulations. Secondly, we normalise d by the range of the combined distributions excluding the highest and lowest 2. 5% of the values: d/r=X¯TCRpMHC−X¯TCRrange (X¯TCRpMHC; 2. 5−97. 5%∪X¯TCR; 2. 5−97. 5%) The value d/r helps to quantify the magnitude of difference related to the width of the combined distributions. Thirdly, we calculate the total variation difference (tvd) to quantify the difference in the probability distributions: TVD (f1, f2) =12∫|f1 (XTCRpMHC) −f2 (XTCR) |dx Where f1XTCRpMHC is the normalized distribution of all TCRpMHC simulation frames and replicas and f2 (XTCR) the normalized distribution of all TCR simulation frames and replicas. The tvd is normed between 0 and 1 where a tvd value of 0 represents perfect overlap of the distributions and a tvd value of 1 represents no overlap. In contrast to d and d/r the tvd does not have a sign i. e. is always positive. We investigated the six CDR loops according to the IMGT [36] definition as extracted from the STCRDAB database [37] as well as the AB loop of the TCR α-chain which was previously hypothesised as influenced by antigen recognition of the LC13 TCR [7], the loops positioned between the variable and constant TCR domains as the linker of the β-chain (“CβFG”) was hypothesised to be important for TCR mechanosensing [10,11], αA and αΒ helices of Cβ as they were reported to be important for CD3 interaction [38], the F and C-strand of Cα as they were reported to be involved in a possible allosteric TCR signalling mechanism [39], and Cα DE and Cβ CC’ loops [40]. These regions and their sequence positon in our four different TCRs are summarised in S1 Table. In order to evaluate how likely an observed descriptor difference (in d, d/r and tvd) would be seen by chance we performed permutation tests. We merged the n replicas of TCRpMHC simulations with the n replicas of the TCR simulations into one group of size 2n. From this 2n group we picked randomly and with repetition n members for group 1 and n members for group 2. Then we calculated d, d/r and tvd between group 1 and group 2. We repeated this procedure 1000 times and obtained a distribution of d, d/r and tvd values. Finally we determine the quartile (q) of the observed d, d/r and tvd between TCRpMHC and TCR within these distributions of random boot strapping samples as a indicator of significance. E. g. q = 0. 98 for d means that 98% of all randomly created group pairs have a smaller d value between them than the d value between TCRpMHC and TCR simulations. An example for a difference in CDR1 loop distance between TCRpMHC and TCR simulations that is likely to be true for the LC13 TCR (d = 0. 9 Å and q = 1. 0 (i. e. none of the 1000 permutations had a larger difference) ) and unlikely to be true the JM22 TCR (d = 0. 3 Å and q = 0. 7 (i. e. about 300 of the 1000 permutations had a larger difference) is given in S2 Fig.
The distances (DIST) between the CDR loops of a TCR can be a descriptor for an opening or closing of the TCR binding interface. Especially the CDR3α and CDR3β loops that are positioned centrally over the scanned peptide (Fig 1B) and could be further apart as the pMHC presses between them or could be closer together as a result of binding interface rigidification. Our results show that the first is true for the LC13 and JM22 TCR while the A6 TCR shows a wider distribution and higher distance for TCR simulations (Fig 2). The 1G4 TCR does not show significant changes. The CDR1 and CDR2 distances also do not show patterns conserved across TCRs (Table 2). The radius of gyration (RG) is a proxy for the compactness of a structure. A high RG indicates a more extended structure while a low RG indicates a more compact structure. We measure the RG of all atoms of the six CDR regions as well as the ABloop, the variable/constant domain linkers and 6 further regions within the TCR constant domains to investigate if pMHC presence has “cramping” effect on any of these regions that could be involved in signalling. Several regions and TCRs show significant differences between their TCRpMHC and TCR states (Table 2); however, these are not conserved. The LC13 and JM22 TCR tend to have lower CDR RGs in their TCRpMHC simulations than in their TCR simulations while the opposite tends to be the case for A6 and 1G4. As an example we show the RG of CDR3α in Fig 3: All differences found are highly significant but they have opposing signs. In contrast the RG of the ABloop and variable/constant linkers are almost unaffected (Table 2). The solvent accessible surface area (SASA) quantifies the extent that a region is exposed to solvent (here measured using the gmx sasa method of Gromacs). When two proteins bind the solvent accessible area of the binding interface will be reduced. This is the case for all six CDRs of all four TCRs upon pMHC binding (Table 2). The more interesting question is if there is also a change in the SASA if the SASA is measured as if no pMHC would be present for TCRpMHC simulations. i. e. is the protruding of the CDRs altered by pMHC presence? Here we obtain a picture that is partly similar to the RG-analysis. The A6 and 1G4 TCR which tend to have higher RGs in TCRpMHC simulations tend also to have increased SASAs. However, for the LC13 and JM22 TCR the reduced RG in TCRpMHC simulations seems not to lead to a decreased SASA. Fig 4 shows this effect for CDR3α. The root mean square fluctuation (RMSF) is an indication of how stable areas of a structure are. A potential signal transduction could be the increased or decreased flexibility of an area. We have investigated the RMSF of all TCR residues in Fig 5. Areas of special interest are marked with dashed lines. For these areas permutation tests are given in Table 2. As expected the RMSFs of the CDRs are lower in TCRpMHC simulations than in TCR simulations due to the restricted degrees of freedom in the binding interface. The changes in CDR RMSF between TCR and TCRpMHC simulations for the LC13, JM22 and 1G4 TCRs are high and mostly significant while the changes for the A6 TCR are lower and mostly not significant (compare Fig 5 and Table 2). In contrast by far the largest observed difference in the RMSF is found in the ABloop of the A6 TCR. In TCRpMHC simulations this loop is about 50% more flexible than in TCR simulations. This is not the case for the other three TCRs. Note that the difference in ABloop arrangement between TCR and TCRpMHC was originally described for the LC13 TCR [7] and not for the A6 TCR as observed here. We also investigated the number of H-bonds between the TCR chains as a proxy for spatial re-arrangement between the TCR chains (Fig 6). For the JM22 TCR there are on average 1. 37 H-bonds less between the TCR chains for TCRpMHC simulations than for TCR simulations. For the LC13 TCR this number is also slightly reduced by 0. 51 H-bonds. For the A6 TCR there is a change of 0. 61 H-bonds in the opposite direction but due to wider variability this number is not significant based on permutation tests (quartile 0. 76). Also for the 1G4 TCR the change is insignificant. To investigate the overall H-bond network in TCRpMHC and TCR simulations we used pyHVis3D [14] which creates a three dimensional graphical representation of the H-bond distributions (Fig 7). This analysis shows different pictures for the four different TCRs. LC13 TCRpMHC simulations show an increased H-bond presence in the area around CDR3α and CDR3β as compared to LC13 TCR simulations which might indicate an interface rigidification. This effect is also present in the JM22 TCR but it is much less pronounced than in the LC13 TCR. In the JM22 TCR the area around the linker between the variable and constant region of the TCRβ chain is mainly affected by pMHC presence causing an increased presence of H-bonds which is not observeable in the LC13 TCR. A6 TCR simulations without pMHC show a higher H-bond frequency in the CDRs of the TCR α-chain while intra TCR chain H-bonds in this area are increased for TCRpMHC simulations. In the constant area of the TCR α-chain H-bonds dominate for TCRpMHC simulations. Similar to JM22 also 1G4 simulations show a higher H-bonds frequency around the linker between the variable and constant region of the TCRβ chain for TCRpMHC simulations. In contrast to all other TCRs the H-bond patterns of the 1G4 TCR around the CDR loops show a very mixed picture. Taken together this shows that there are no conserved differences in the H-bond patterns between TCR and TCRpMHC simulations.
We have presented an MD study of four different TCRs (LC13, JM22, A6, and 1G4) in their pMHC bound and unbound form using a total of 26 000 ns. To our knowledge this is the first study that investigates four different TCRs on such a large scale. The most similar study was published by Hawse et al. [41] who investigated the A6 and DMF5 TCR and found a global rigidification and dampened coupling in the linker between variable and constant TCR domains upon pMHC binding using computational mutagenesis with gradient-based minimization [27] and experimental hydrogen/deuterium exchange. Also Cuendet et al. [18] investigated the detachment of the A6 and B7 TCR from HLA-A2 using about 4 000 ns of steered MD code. This study found interesting binding interface characteristics but did not address the different dynamics within the TCR in the bound and unbound form. In two studies we have previously investigated the LC13 TCR [21,23]. In one of these studies [21] we found multiple significant differences between the pMHC bound and unbound LC13 TCR. These included CDR distance distributions, CDR compactness as well as differences in the TCR hydrogen bond network. In the current study we investigated if these results hold true for other TCR/MHC combinations. Surprisingly, we obtained very different results for JM22 TCR/HLA-A*02: 01 compared with LC13 TCR/ HLA-B*08: 01. Thus the linker between the C and V regions of the TCR β-chain had a lower RG and RMSF for JM22 TCRpMHC simulations. Furthermore a reduced RG was seen in CDR3α and CDR3β while it was increased in CDR2β. Because LC13 and JM22 bind different MHCs, we investigated two more TCRs (A6 and 1G4) which bind the same MHC as JM22. Even when binding the same MHC type the dynamics within the TCRs vary significantly in, for example, RGs, SASAs, RMSFs (Table 2). These differences between TCRs upon pMHC binding are consistent with our experimentally measured finding that the JM22, A6 and 1G4 TCRs have very different energetic footprints on HLA-A*02: 01 [15]. Other experimental studies support the conclusions drawn in our study. With regard to changes at the binding interface, several structural studies have demonstrated local conformational changes upon TCR binding to pMHC [1], while thermodynamic and kinetic studies of several TCRs, including JM22 [42], are consistent with a reduction in conformational flexibility upon binding. With respect to conformational changes distal from the interface, while the AB-linker of the LC13 TCR has been described to differ between the X-ray structure of the HLA-B*08: 01 bound and unbound structure [7], for other TCRs this finding could not be replicated. Similarly the B4. 2. 3 TCR was reported to be effected in its H3 loop by the binding of H2-Dd presenting a 10-mer HIV-env peptide [39]. But no further support for this being a conserved mechanism could be found in any of the other 10 TCRpMHC complex structures for which a separate TCR structure exists [39]. While our finding argue against pMHC binding inducing conformation changes through allosteric mechanism in the TCR, they remain consistent with models proposing that conformational changes are introduced by mechanical mechanisms [4,5, 8–11]. For example, pulling [8,9] and shearing [10] forces have been shown to enhance TCR triggering; with the CβFG loop seemingly be affected by side-wards pulling on the TCR [11]. We conclude that TCR structural dynamics do not differ between TCR/pMHC and TCR simulations in conserved and non-obvious ways. Taken together with previous studies our findings argue against a role for allosteric conformation change models in TCR triggering. Our results are consistent with mechanical models of conformational change, as well as aggregation and kinetic-segregation models.
The interaction between T-cells and other cells is one of the most important interactions in the human immune system. If T-cells are not triggered major parts of the immune system cannot be activated or are not working effectively. Despite many years of research the exact mechanism of how a T-cell is initially triggered is not clear. One hypothesis is that conformational changes within the T-cell receptor (TCR) can cause further downstream signalling within the T-cell. In this study we computationally investigate the dynamics of four different TCRs in their free and bound configuration. Our large scale simulations show that all four TCRs react to binding in different ways. In some TCRs mainly the areas close to the binding region are affected while in other TCRs areas further apart from the binding region are also affected. Our results argue against a conserved structural activation mechanism across different types of TCRs.
blood cells
medicine and health sciences
immune cells
immunology
simulation and modeling
clinical medicine
mathematics
protein structure
tcr signaling cascade
discrete mathematics
combinatorics
research and analysis methods
immune system proteins
white blood cells
major histocompatibility complex
animal cells
proteins
t cells
molecular biology
biochemistry
signal transduction
t cell receptors
biochemical simulations
macromolecular structure analysis
permutation
cell biology
clinical immunology
biology and life sciences
immune receptors
cellular types
physical sciences
cell signaling
computational biology
signaling cascades
MHC binding affects the dynamics of different T-cell receptors in different ways