T-Cell Receptor Binding Affects the Dynamics of the Peptide/MHC-I Complex.

The recognition of peptide/MHC by T-cell receptors is one of the most important interactions in the adaptive immune system. A large number of computational studies have investigated the structural dynamics of this interaction. However, to date only limited attention has been paid to differences between the dynamics of peptide/MHC with the T-cell receptor bound and unbound. Here we present the first large-scale molecular dynamics simulation study of this type investigating HLA-B*08:01 in complex with the Epstein-Barr virus peptide FLRGRAYGL and all possible single-point mutations (n = 172). All of the simulations were performed with and without the LC 13 T-cell receptor for a simulation time of 100 ns, yielding 344 simulations and a total simulation time of 34 400 ns. Our study is 2 orders of magnitude larger than the average T-cell receptor/peptide/MHC molecular dynamics simulation study. This data set provides reliable insights into alterations of the peptide/MHC-I dynamics caused by the presence of the T-cell receptor. We found that simulations in the presence of T-cell receptors have more hydrogen bonds between the peptide and MHC, altered flexibility patterns in the MHC helices and the peptide, a lower MHC groove width range, and altered solvent-accessible surface areas. This indicates that without a T-cell receptor the MHC binding groove can open and close, while the presence of the T-cell receptor inhibits these breathing-like motions.


■ INTRODUCTION
The interaction between T-cell receptors (TCRs) and protein fragments (called peptides or antigens) presented by major histocompatibility complexes (MHCs)/β 2 microglobulin (β 2 M) is one of the most important events in the adaptive immune system. 1 Intracellular proteins are degraded by proteasomes into peptides. These peptides are transported into the endoplasmic reticulum (ER) and bound inside the binding groove of MHC complexes. The peptide/MHC complexes are then presented on the surface of the cell, where they can be recognized by the TCRs of T-cells. Only limited immune response can take place against a peptide if (1) this peptide cannot be presented by MHCs or (2) TCRs are not triggered by peptide/MHC (pMHC) combinations. Although much research on TCR triggering and subsequent signal transduction via CD3 has been carried out, the detailed structural mechanism remains elusive. 2 Some theories 2−5 suggest subtle structural alterations or differences in dynamics. The involvement of structural dynamics in TCR triggering is supported by several lines of evidence (reviewed in ref 2): (i) a mechanical pushing or pulling force on the TCR is inevitable during pMHC binding and unbinding; 6 (ii) CD3 is relatively rigid, and the FG loop of the TCR constant β-domain is in direct contact with the membrane-distal end of CD3, 7 so mechanical forces acting on the TCR could be directly transmitted into the cell; (iii) T-cell activation is optimal if the pMHC is anchored to an artificial surface; 8 (iv) elongation of the pMHC ectodomain reduces TCR triggering 9 because only a reduced mechanical force acts on the TCR; (v) mechanical forces on the TCR enhance TCR triggering. 10 Dynamic changes at the atomistic level cannot be measured using current experimental techniques such as X-ray crystallography or NMR spectroscopy. Instead computational molecular dynamics (MD) studies have been carried out to investigate the dynamics of pMHC and TCRpMHC structures (reviewed in ref 11; median total simulation time of 135 ns). For example, studies have investigated the effects of different peptides bound by the same MHC for HLA-B*27, 12 HLA-DR1, 13 HLA-B*08, 14,15 or the mouse allele I-A u . 16 Other studies have analyzed the effect of different MHCs for the same set of peptides 5,13 and the effect of different TCRs on pMHC 6,17 or TCR pulling. 6 Only limited attention has been paid to a comparison of the dynamics of TCR-bound and -unbound pMHC. However, experimental studies have investigated these effects. Examples include the X-ray structure comparison of FLRGRAYGL/HLA-B*08:01 with the LC13 TCR 18 and without a TCR 19 as well as solution mapping of the 2C TCR docking footprints on the pMHC QL9/H-2L d . 20 Only two MD studies have investigated changes in a pMHC in reaction to the presence of a TCR. One study 21 was based on two 20 ns simulations of HLA-B*44:05 to investigate the cis−trans isomerization of the MHC-presented peptide after DM1 TCR recognition. The second study 4 performed four 100 ns simulations to investigate HLA-B*35:08 with and without the SB27 TCR.
This limited attention is surprising, as the combined surface of the peptide and MHC is the only site where the T-cell physically interacts with the peptide antigen and hence this site is likely to be an initial trigger of the TCR signaling cascade.
In this paper, we examine the differences between the dynamics of pMHCs in their TCR-bound and -unbound forms using a large set of MD simulations of these complexes. We have carried out 172 pMHC MD simulations and 172 TCRpMHC MD simulations, 14 each 100 ns in length, giving a total simulation time of 34 400 ns. This data set is about 250 times larger than is used in the average study in this area. 11 Using these data, we find that TCRpMHC simulations differ significantly from pMHC simulations in terms of hydrogen bonds, residue flexibility, MHC binding groove width, and solvent-accessible surface area. Our data suggest that the pMHC alone can undergo breathing-like movements of the MHC groove while the TCRpMHC complex cannot.

■ METHODS
Modeling of the pMHC and TCRpMHC Structures. The X-ray structure of the LC13 TCR in complex with HLA-B*08:01-bound Epstein−Barr virus peptide FLRGRAYGL (PDB accession code 1mi5) was used as the basis for this study.
It has previously been shown that comparisons between two or a few simulations are likely to produce misleading results. 14,22−24 For reliable conclusions, a larger number of simulations is needed. To provide a reliable analysis of changes in the pMHC interface in reaction to the presence/absence of

Journal of Chemical Information and Modeling
Article the TCR, we performed all possible single-point mutations of the wild-type peptide FLRGRAYGL in the X-ray structure, leading to a total of 172 structures (9 positions times 19 mutations plus the wild type). This procedure was carried out twice, once with the TCR present and once without, yielding 344 complexes in total. All of the mutations were carried out using SCWRL 25 via the peptX framework, 26 as we have previously shown that this method is most appropriate for pMHC mutations. 27,28 MD Simulations. We performed 172 pMHC simulations and compared them with the 172 TCRpMHC simulations of our previous study. 14 All 344 MD simulations were performed with the same parameters using GROMACS 4 29 and the GROMOS force field in the 53a6 parametrization. 30 Each complex was immersed in a dodecahedronic simulation box filled with explicit SPC water with a minimum distance of 1.5 nm between the (TCR)pMHC and the box boundary. A dodecahedronic simulation box was used because it has only about 71% of the volume of a cubic simulation box with the same periodic distance. Periodic boundary conditions were applied to the box. Random water molecules were replaced with Na + and Cl − ions to achieve a salt concentration of 0.15 mol/L and a neutral charge. Each system was energetically minimized to avoid initial energies at a too-high level and then warmed to 310 K using position restraints. MD simulations of 100 ns per complex were carried out using the Oxford Advanced Research Computing (ARC) facility. The real-space cutoff for Coulomb interactions was set to 1.4 nm. Long-range electrostatic interactions were calculated using the particlemesh Ewald method. Short-range neighborhood lists were updated every 10 steps. The distance for the Lennard-Jones cutoff was set to 1.4 nm. Initial velocities were generated according to a Maxwell distribution at a temperature of 310 K. The total simulation time was 34 400 ns (34.4 μs).
Methods of Trajectory Evaluation. The following GROMACS standard tools were used: g_rmsf to compute the root-mean-square fluctuation (RMSF), g_hbond for the hydrogen bonds (H-bonds), g_sas for the solvent-accessible surface area (SASA), and g_dist as well as the gro2mat

Journal of Chemical Information and Modeling
Article package 31 for the distance calculations. Simulations were manually inspected with VMD 32 and the vmdICE plugin. 33 Comparisons of TCR-Bound and -Unbound pMHC Simulations. The total variation distance (tvd) (defined in a previous study 14 ) was used to quantify how the distributions of parameters deduced from TCR-bound and -unbound pMHC MD simulations differ. The tvd is defined as where f 1 (x) is the first normalized distribution and f 2 (x) the second normalized distribution. Thus, the tvd values can range between 0 and 1. A value of 0 represents perfect overlap (0% variation) of the distributions, while a value of 1 represents no overlap (100% variation). As the tvd does not give a quantification of the difference in means between the two distributions, we calculated in addition a normalized distance between the means according to our previous study. 14 This value is denoted as d/r and defined as where X̅ 1 and X̅ 2 are the mean values of the two distributions (TCRpMHC and pMHC) and the denominator is the range of the combined distributions excluding the lowest and highest 2.5%. Both tvd and d/r were used for H-bonds, distance measurements, and the SASA.

■ RESULTS
In total we analyzed 344 simulations (172 pMHC and 172 TCRpMHC), each 100 ns in length each. All of the simulations were performed with the GROMACS MD package using identical parameters. In each simulation the same MHC and TCR are present in complex with a different peptide. As an example, we show the simulation of the wild-type peptide FLRGRAYGL bound by HLA-B*08:01 with the LC13 TCR (movie S1) and the simulation of FLRGRAYGL bound by HLA-B*08:01 without a TCR (movie S2).
Hydrogen Bonds between the Peptide and MHC in Reaction to TCR Binding. Hydrogen bonds play an important role in stabilizing peptides inside the MHC binding groove ( Figure 1A; also see ref 1). We analyzed the number of H-bonds between the peptide and MHC for each of our 344 simulations ( Figure 1B). It can be seen that simulations in the presence of the TCR show a higher number of H-bonds between the peptide and MHC than simulations without the TCR. The presence of the TCR appears to aid the occurrence of H-bonds between the peptide and MHC.
To investigate which peptide residues contribute most to this change in the number of H-bonds, we analyzed each peptide residue separately ( Figure 1C). Peptide positions 5 and 9 show the largest number of H-bonds to the MHC. At the same time, they also have the largest difference between the numbers of Hbonds in the trajectories with and without the TCR (tvd values of 20% and 23%, respectively). These positions are the experimentally known anchor residues for HLA-B*08:01, which are essential for peptide/MHC binding. 34 The increased number of H-bonds between peptide anchor residues and MHC pockets in the TCRpMHC simulations indicates that the presence of the TCR increases peptide/MHC binding.
Flexibility of the pMHC in Reaction to TCR Binding. Previous studies have suggested that flexibility in the pMHC is important for TCR recognition. 5, 35 We investigated whether and how the presence of the TCR affects the flexibility of the pMHC by calculating the RMSF of the peptide and MHC on a per-residue basis. Figure 2A shows that secondary structure elements such as the β-strand floor and the two α-helices are in general more rigid than the loop regions. The presence of the TCR further lowers the RMSF for large areas of the two α-helices that flank the peptide binding groove. In contrast, it increases the flexibility in the linker region between MHC β-strands 1 and 2 as well as in the area between β-strand 3 and the α1-helix of the

Journal of Chemical Information and Modeling
Article MHC (Figure 2A). Residues not forming parts of the MHC head show almost identical RMSF patterns for simulations with and without the TCR. Taking all atoms instead of only Cα into account shows similar but slightly noisier patterns (compare Figure 2A and Figure 2B).
The presence of the TCR lowers the flexibility of the peptide backbone ( Figure 2C) as well as that of side-chain atoms ( Figure 2D). The high side-chain flexibility of peptide residue 7 is curbed by the presence of the TCR ( Figure 2D). Peptide position 7 is known to protrude deeply into the TCR in the Xray structure. 18 Therefore, it is unsurprising that the presence of the TCR curbs its flexibility. Figure 2E illustrates regions with increased and decreased residue flexibility on a three-dimensional (3D) representation of the pMHC. From the viewpoint of the TCR it can be seen that most of the regions that are oriented toward the TCR have decreased flexibility when the TCR is present (blue in Figure  2E). This is expected, as the presence of the TCR reduces the conformational freedom of these residues. It is, however, surprising that some regions acting as "suspenders" for the MHC head have increased flexibility when the TCR is present (red in Figure 2E). It appears that these "suspenders" compensate for the loss of flexibility in the TCR interface regions.
Width of the MHC Binding Groove in Reaction to TCR Binding. Hydrogen bonds between the peptide and MHC as well as the flexibility of the pMHC are affected by the presence of the TCR. To investigate how the loss of H-bonds and increased flexibility influence the width of the MHC binding groove, we measured the distance between the MHC helices for three different regions ( Figure 3A).
The MHC groove width at the peptide's N-terminal end is not affected by the presence of the TCR ( Figure 3B). This agrees with Figure 2E, which shows that the RMSF of this region does not change in reaction to the presence of the TCR. In both simulations with the TCR and simulations without the TCR, the groove width is almost normally distributed (mean = 16.4 Å) with a longer tail to the right.
The distance between the evolutionarily conserved kinks 36 of the two MHC helices shows an almost normal distribution ( Figure 3C). However, simulations without the TCR show a larger range of values than simulations with the TCR. This larger range of values and the increased RMSF of the MHC helices ( Figure 2E) suggest that the groove can open and close slightly when not bound to a TCR. TCR binding appears to hamper this opening and closing of the MHC groove. Figure 3D shows what may be a bimodal distribution of distances between the MHC helices at the peptide's N-terminal end. The closed state (left peak) has a mean value of 14 Figure S1).
Overall, these results suggest that the MHC binding groove can open and close both in the middle and at the end that binds to the peptide's C-terminus if no TCR is present, but this is hampered in the presence of a TCR.
Solvent-Accessible Surface Area of the pMHC in Reaction to TCR Binding. When any protein binds to another, one would expect the protein to have a reduced solvent-accessible surface area. 37 This is the case for pMHCs:

Journal of Chemical Information and Modeling
Article TCR-free pMHC/β 2 M complexes show an average SASA of 228 nm 2 , while the average SASA for TCR-bound pMHC/β 2 M complexes drops to 219 nm 2 . This means that the TCR covers on average only about 9 nm 2 or 4% of the pMHC/β 2 M surface (a small interface area in comparison with other proteins 38 ). This low change in SASA agrees with the experimentally known low-affinity binding between TCR and pMHC.
The initial assumption would be that the TCR causes the change in SASA solely by forcing water molecules away from the pMHC. Figure 4 shows that this is not the case. In Figure  4A we compare the peptide SASA of simulations without the TCR against that of simulations with the TCR. In Figure 4B we compare the peptide SASA of simulations without the TCR against that of simulations with the TCR in which the SASA was measured as if the TCR were not there (i.e., the simulation trajectory was calculated with the TCR but the TCR was removed from the trajectory for the SASA analysis). If the TCR had no influence on the MHC and peptide conformation, these two distributions would be identical, but they are not. The peptide SASA drops from 4.927 nm 2 (pMHC) to 4.582 nm 2 (TCRpMHC) (a decrease of 7%; Figure 4B). This indicates that the presence of the TCR stabilizes peptide/MHC interactions by hampering water from entering the MHC binding groove.
All 172 Peptides Are Ligands for HLA-B*08:01. All of our simulations started with the modeled peptide inside the MHC binding groove. It is likely that this is a valid assumption for two reasons. First, we modeled our TCRpMHC complexes on the basis of an experimental X-ray structure. The wild-type peptide (FLRGRAYGL) is known to be a strong binder of HLA-B8. Each of the 171 peptide mutants differs from the Xray structure by only a single-point mutation. According to the SYFPEITHI 34 database, HLA-B8 has three anchor residues (positions 3, 5, and 9). All three anchor residues match in the wild-type peptide. When we perform a single mutation of this peptide, two of the three anchors remain intact. On this basis it is likely that all 172 peptides would bind HLA-B8 to greater or lesser degrees. Second, to further support this argument, we used the IEDB T-cell epitope prediction tool 39 on our peptides. The results show that all 172 peptides have a percentile rank of ≤6.5 and that 165 peptides have a score of <3 (a score closer to 0 indicates higher binding affinity, and a score closer to 100 indicates lower binding affinity). This supports our assumption that all of our peptides are at least weak if not moderate/strong HLA-B8 binders.

■ DISCUSSION
Here we have presented by far the largest MD study investigating differences between 172 pMHC and 172 TCRpMHC simulations (total simulation time of 34 400 ns). This is more than 250 times larger than the average simulation study in the field (as reviewed in ref 11). These data allowed us to investigate the effects of the presence of the LC13 TCR on HLA-B*08:01 beyond single examples. Often a larger number of simulations is needed for reliable conclusions. This was for example shown for TCRpMHC interactions, 14 peptide/MHC binding predictions, 24 and HIV-1 protease inhibitors 23 and also holds true for the current pMHC/TCRpMHC comparison (compare Figure S2 with Figure 1, Figure 2, and Figure 3).
We have found several conserved effects of the presence of the LC 13 TCR in HLA-B*08:01 simulations. First, there is a higher number of H-bonds between the peptide and MHC in the TCRpMHC simulations than in the pMHC simulations.
The increased number of H-bonds is in accordance with our second finding, that the MHC helix RMSF and the peptide RMSF are lower in the TCRpMHC simulations than in the pMHC simulations. The TCR hampers the peptide from its frequently observed up and down flapping of the N-or Cterminal end. This flapping behavior has been previously observed in MD simulations for HLA-B*27:05, 40 HLA-A*02:01, 22 and HLA-DR. 13,35 In further agreement with the decreased number of H-bonds and the increased RMSF for the pMHC simulations as opposed to the TCRpMHC simulations, we found that the MHC groove can open slightly wider in the pMHC simulations than in the TCRpMHC simulations. We find this behavior of the MHC groove at the middle of the peptide and at the peptide's Cterminal end. We do not find changes in the groove width at the peptide's N-terminal end.
The conserved differences between pMHC and TCRpMHC are in agreement with analyses of available MHC class I X-ray structures (321 pMHC and 52 TCRpMHC), which found decreased torsion of the MHC α1-helix of TCRpMHC structures in contrast to pMHC structures. 41−43 This decreased torsion is likely to be the result of TCR engagement.
Our findings could also be taken as an indication that peptide/MHC binding or unbinding starts C-terminally. This might be MHC-allele-specific, as a peptide/MHC detachment study has recently shown that peptides do not have a preference for N-or C-terminal detachment for HLA-A*02:01. 22

■ CONCLUSION
Although a large number of MD studies have recently investigated the dynamics of (TCR)pMHC interactions, only limited attention has been paid to the effects of the presence of the TCR on the pMHC dynamics. Here we have presented the first large-scale study based on 344 simulations, each 100 ns in length. We find that the number of H-bonds between the peptide and MHC is higher for TCRpMHC than for pMHC, that the presence of the TCR decreases the flexibility in the peptide and large parts of the MHC α-helices but increases the flexibility in two loops of the MHC β-sheet floor, that the MHC binding groove can widen and narrow at the peptide's Cterminal end if no TCR is present, and that the presence of the TCR has structural effects on the peptide that make it less solvent-accessible. Together these results indicate that the TCR-free pMHC underwent a breathing-like movement of the MHC helices that is hampered if a TCR is present.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.5b00511.
Detail analysis of Figure 3D ( Figure S1) and a misleading comparison of one simulation against another simulation ( Figure S2) (PDF) Molecular dynamics simulation of the wild-type peptide FLRGRAYGL bound by HLAB*08:01 and the LC13 TCR (TCRpMHC system) (AVI) Molecular dynamics simulation of the wild-type peptide FLRGRAYGL bound by HLAB*08:01 (pMHC system) (AVI)

Journal of Chemical Information and Modeling
Article