Subscriber access provided by CMU Libraries - http://library.cmich.edu
Article
T-cell Receptor Binding Effects the Dynamics of the Peptide/MHC-I Complex Bernhard Knapp, and Charlotte M Deane J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.5b00511 • Publication Date (Web): 03 Dec 2015 Downloaded from http://pubs.acs.org on December 8, 2015
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
T-cell Receptor Binding Effects the Dynamics of the Peptide/MHC-I Complex
Bernhard Knapp1*, Charlotte M. Deane1
1
Department of Statistics, Protein Informatics Group, University of Oxford
*
[email protected] Abstract 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, so far 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 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 two orders of magnitude larger than the average T-cell receptor / peptide / MHC Molecular Dynamics simulation study. This dataset provides reliable insights into alteration in peptide/MHC-I dynamics caused by T-cell receptor presence. We found that simulations in the presence of T-cell receptors have more H-bonds between peptide/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 T-cell receptor presence 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 (MHC) / β2 microglobulin (β2M) is one of the most important events in the adaptive immune system 1. Intracellular proteins are degraded by
ACS Paragon Plus Environment
1
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 21
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/MHCs (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 theories2,3,4,5 suggest subtle structural alterations or differences in dynamics. Involvement of structural dynamics in TCR triggering is supported by several lines of evidence (reviewed in 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; mechanical forces acting on the TCR could be directly transmitted into the cell; iii) T-cell activation is optimal if 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; and v) mechanical forces on the TCR enhance TCR triggering 10. Dynamic changes at atomistic level cannot be measured using current experimental techniques such as X-ray crystallography or NMR. Instead computational Molecular Dynamics (MD) studies have been carried out to investigate the dynamics of pMHC and TCRpMHC structures (reviewed in 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 IAu 16. Other studies have analysed the effect of different MHCs for the same set of peptides 5,13, the effect of different TCRs on pMHC 6,17 or TCR pulling 6. Only limited attention has been paid to the comparison of the dynamics of TCR bound and unbound pMHC. Experimental studies have however investigated these effects. Examples include the X-ray structure comparison of FLRGRAYGL/HLAB*08:01 with the LC13 TCR 18 and without TCR 19 as well as solution mapping of the 2C TCR docking footprints on the pMHC QL9/H–2Ld 20. Only two MD studies have investigated changes in pMHC in reaction to TCR presence. 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 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 signalling cascade.
ACS Paragon Plus Environment
2
Page 3 of 21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
In this paper we examine the differences in the dynamics between 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 of 100 ns length each giving a total simulation time of 34 400 ns. This dataset is about 250 times larger than 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 suggests that pMHC alone can undergo breathing-like movements of the MHC groove while the TCRpMHC complex cannot.
Methods Modelling 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 few simulations are likely to produce misleading results 14,22,23,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 TCR presence/absence we performed all possible single point mutations to the wild-type peptide (FLRGRAYGL) of the X-ray structure. This leads to a total of 172 structures (9 positions by 19 mutations and the wild-type). This procedure was carried out twice, once with TCR present and once without, yielding 344 complexes in total. All mutations were carried out using SCWRL 25 via the peptX frame work 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 to 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 parameterisation 30. Each complex was immersed into 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 as 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/litre and a neutral charge. Each system was energetically minimised to avoid initial energies at s too high level and then
ACS Paragon Plus Environment
3
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 21
warmed up 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 cut-off for coulomb interactions was set to 1.4 nm. Long range electrostatic interactions were calculated using the Particle-Mesh Ewald method. Short range neighbourhood lists were updated every 10 steps. The distance for the Lennard-Jones cut-off was set to 1.4 nm. Initial velocities were generated according to a Maxwell distribution at temperature 310K. The total simulation time was 34 400 ns (34.4 µs).
Methods of trajectory evaluation The GROMACS standard tools g_rmsf was used 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 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 distributions of parameters deduced from TCR bound and unbound pMHC MD simulations differ. The tvd is defined as:
tvd ( f1 , f 2 ) =
1 | f1 ( x ) − f 2 ( x ) |dx 2∫
where f1(x) is the first distribution normalized and f2(x) the second distribution normalized. 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 referred to as d/r and defined as:
d /r =
| X1 − X 2 | range( X 1; 2.5−97.5% ∪ X 2; 2.5−97.5% )
where ܺതଵ and ܺതଶ are the mean value 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.
ACS Paragon Plus Environment
4
Page 5 of 21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
Results In total we analysed 344 simulations (172 pMHC and 172 TCRpMHC) of 100 ns length each. All simulations were performed with the MD package GROMACS 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 peptide and MHC in reaction to TCR binding Hydrogen bonds play an important role in stabilising peptides inside the MHC binding groove (Figure 1A and Janeway et al. 1). We analysed the number of hydrogen bonds between 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 peptide and MHC than simulations without TCR. TCR presence appears to be aiding the occurrence of H-bonds between peptide and MHC. To investigate which peptide residues contribute most to this change in the number of H-bonds we analysed each peptide residue separately (Figure 1C). Peptide positions five and nine show the largest number of H-bonds to the MHC. At the same time they also have the largest difference in Hbonds between trajectories with and without TCR (tvd 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 of TCRpMHC simulations indicates that TCR presence increases peptide/MHC binding.
ACS Paragon Plus Environment
5
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 21
Figure 1: Hydrogen bonds between peptide and MHC. (A) Illustration of H-bonds on the basis of the X-ray structure of PDB accession code 1mi5. Blue licorice: peptide; White licorice: MHC residues within 3 Å of the peptide; White and transparent cartoon: MHC; Dashed red: H-Bonds. The MHC α1-helix (R45-90) is not shown as cartoon to allow insight into the MHC binding groove from the side. (B) Distribution of H-bonds between peptide and MHC over 344 MD simulations split by the presence and absence of the TCR. (C) Distribution of the H-bonds in the 344 MD simulations over the nine peptide residue positions shown in (A).
Flexibility of pMHC in reaction to TCR binding Previous studies have suggested that flexibility in the pMHC is important for TCR recognition 5,35. We investigated if and how TCR presence affects pMHC flexibility 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. TCR presence further lowers the RMSF for large areas of the two α-helices that flank the peptide binding groove. In contrast TCR presence increases the flexibility in the linker-region between MHC β-strand one and two as well as in the area between β-
ACS Paragon Plus Environment
6
Page 7 of 21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
strand three and the α1-helix of the MHC (Figure 2A). Residues not forming parts of the MHC head show almost identical RMSF patterns for simulations with and without TCR. Taking all atoms instead of only Cα into account shows similar, but slightly noisier, patterns (compare Figure 2A and Figure 2B). TCR presence lowers the flexibility of the peptide backbone (Figure 2C) as well as side-chain atoms (Figure 2D). The high side-chain flexibility of peptide residue seven is curbed by the presence of the TCR (Figure 2D). Peptide position seven is known to protrude deeply into the TCR in the X-ray structure 18. Therefore it is unsurprising that TCR presence curbs its flexibility. Figure 2E illustrates regions with increased and decreased residue flexibility on a three dimensional representation of pMHC. From the viewpoint of the TCR it can be seen that most regions that are oriented towards the TCR have a decreased flexibility if a TCR is present (Figure 2E in blue). This is expected as TCR presence reduces the conformational freedom of these residues. It is, however, surprising that some regions acting as “suspenders” for the MHC head have an increased flexibility if the TCR is present (Figure 2E in red). It appears that these “suspenders” compensate for the loss of flexibility in the TCR interface regions.
ACS Paragon Plus Environment
7
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 21
Figure 2: Root mean square fluctuation of pMHC over the 344 MD simulations split by the presence and absence of the TCR. (A) Backbone RMSF of the MHC. The vertical dashed lines separate secondary structure elements: The “s” labels mark β-strands. The first 8 β-strands form the MHC floor (B) As in (A) but all-atom RMSF of the MHC. (C) Backbone RMSF of the peptide. (D) All-atom RMSF of the peptide. The dotted lines in A-D represent the standard error of mean. (E) Visualisation of RMSF differences of (A) on a 3D representation of pMHC. For clarity only the MHC head (residues 1 to 180) is shown. Blue: regions with increased RMSF for pMHC simulations; Red: regions with increased RMSF for TCRpMHC simulations; White: regions with roughly equal RMSF between pMHC and TCRpMHC.
Width of the MHC binding groove in reaction to TCR binding H-bonds between peptide/MHC as well as the flexibility in pMHC are affected by the TCR presence. To investigate how loss of 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 TCR presence (Figure 3B). This agrees with Figure 2E which shows that the RMSF of this region does not change in reaction to TCR presence. In both simulations with TCR and simulations without 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 almost normally distributed distances (Figure 3C). However, simulations without TCR show a larger range of ACS Paragon Plus Environment
8
Page 9 of 21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
values than simulations with 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 Å while the open state (right peak) has a mean value of 15.2 Å. If we define a change in distance from ≤ 14 Å to ≥ 15.2Å or vice versa during simulation time 10-100 ns as a state transition then pMHC simulations undergo on average 4.8 state transitions and spend 18.1% of their simulation time in the closed state, 42.5% in the open state, and 39.5% in between. TCRpMHC simulations under go on average 3.4 state transitions and spend 37.8% in the closed state, 18.3% in the open state, and 43.9% in between (Figure S1). Overall this suggests that the MHC binding groove can open and close both in the middle and at the end which binds to the peptide’s C-terminal if no TCR is present, but this is hampered in the presence of a TCR.
Figure 3: MHC binding groove width over 344 MD simulations split by the presence and absence of the TCR. (A) Visualisation of the distances shown in (B), (C), and (D) on an MHC molecule. For (B) the distance between the average coordinates of the yellow regions is used. For (C) the distance between the average coordinates of the evolutionary conserved MHC kink points (red) is used. For (D) the distance between the average coordinates of the green regions is used. (B) Distance distribution between the starting part of the α1-helix and the ending part of the α2-helix. (C) Distance distribution between the central kink of the α-1 helix and the central kink of the α-2 helix. (D) Distance distribution between the ending part of the α1-helix and the starting part of the α2-helix.
Solvent accessible surface area of 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: TCR free pMHC/β2M show an average SASA of
ACS Paragon Plus Environment
9
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 21
228 nm2 while the SASA for TCR bound pMHC/β2M complexes drops to 219 nm2. This means that the TCR covers on average only about 9 nm2 or 4% of the pMHC/β2M surface (a small interface area in comparison to other proteins 38). This low change in SASA agrees with experimentally known low affinity binding between TCR and pMHC. The initial assumption would be that the TCR is causing the change in SASA solely by forcing water molecules away from pMHC. Figure 4 shows that this is not the case. In Figure 4A we compare the peptide SASA of simulations without TCR against simulations with TCR. In Figure 4B we compare the peptide SASA of simulations without TCR against simulations with TCR but the SASA is measured as if the TCR is not there (i.e. the simulation trajectory was calculated with TCR while 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 nm2 (pMHC) to 4.582 nm2 (TCRpMHC) (a decrease of 7%; Figure 4B). This indicates that TCR presence stabilises peptide/MHC interactions by hampering water from entering the MHC binding groove.
Figure 4: Solvent accessible surface area distribution of the peptide over 344 MD simulations split by the presence and absence of the TCR. (A) SASA of the peptide. The presence of a TCR reduces the SASA of the peptide on average by 2 2.1 nm . (B) Same as (A) but the SASA of the TCR bound peptides is calculated as if the TCR would not be present. The slight decrease of SASA in the TCRpMHC trajectories is due to structural changes in the peptide induced by the TCR. Orange transparent cartoon: TCR α-chain; Black transparent cartoon: TCR β-chain; White surface: MHC and β2microglobulin; Blue surface: peptide; Overlaying red and transparent surface: approximation of the SASA area measured.
ACS Paragon Plus Environment
10
Page 11 of 21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
All 172 peptides are ligands for HLA-B*08:01 All our simulations start with the modelled peptide inside the MHC binding groove. It is likely that this is a valid assumption for two reasons. Firstly, we have modelled our TCRpMHC complexes on the basis of an experimental X-ray structure. The wild-type peptide (FLRGRAYGL) is known as a strong binder for HLA-B8. The 171 peptide mutants differ from the X-ray structure only by a single point mutation. According to the SYFPEITHI 34 database HLA-B8 has 3 anchor residues (position 3, 5, and 9). All three anchor residues match in the wild-type peptide. When we perform a single mutation to this peptide, 2 out of 3 anchors remain intact. On this basis it is likely that all of the 172 peptides would bind HLA-B8 to a greater or lesser degree. Secondly, 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