Intramolecular Domain Movements of Free and Bound pMHC and TCR Proteins: A Molecular Dynamics Simulation Study

The interaction of antigenic peptides (p) and major histocompatibility complexes (pMHC) with T-cell receptors (TCR) is one of the most important steps during the immune response. Here we present a molecular dynamics simulation study of bound and unbound TCR and pMHC proteins of the LC13-HLA-B*44:05-pEEYLQAFTY complex to monitor differences in relative orientations and movements of domains between bound and unbound states of TCR-pMHC. We generated local coordinate systems for MHC α1- and MHC α2-helices and the variable T-cell receptor regions TCR Vα and TCR Vβ and monitored changes in the distances and mutual orientations of these domains. In comparison to unbound states, we found decreased inter-domain movements in the simulations of bound states. Moreover, increased conformational flexibility was observed for the MHC α2-helix, the peptide, and for the complementary determining regions of the TCR in TCR-unbound states as compared to TCR-bound states.


Introduction
The interaction between T-cell receptors (TCR) and antigenic peptides (p) presenting major histocompatibility complexes (pMHC) is a crucial step in adaptive immune response [1]. It triggers the generation of cell-mediated immunity to pathogens and other antigens. The response is driven by TCRs specifically recognizing antigenic peptides bound to and presented by MHC (pMHC) of infected or transformed cells. MHC class I molecules are membrane bound proteins on the surface of antigen presenting cells consisting of an α-chain and a β2-microglobulin. Short fragments of antigens are presented in a pocket formed by the α1 and α2 domains. The MHC-bound peptides are presented to the TCR for examination. The TCR itself is a heterodimer on the surface of T-cells and is formed by an α-chain and a β-chain (each chain containing a variable and a constant region) [2]. The binding site of the T-cell receptor is constituted by the two variable regions TCR V α and TCR V β , which specifically recognize the α-helices of the MHC peptide-binding groove as well as parts of the antigen [3]. The interaction of TCR and pMHC is weak but highly specific, therefore being able to distinguish self-antigens from foreign peptides [1]. Additional molecules and co-receptors, e.g., the CD8, a T-cell surface glycoprotein, help establish a stable, specific and sensitive interaction. the LC13 TCR and Ferber et al. [19] analyzed TCR binding orientations over pMHC. In previous papers we used the 3KPS system for geometric analyses of alloreactive HLA α-helices [11,20] as well as for finding semi-rigid domains in MD simulations [21]. In the present study we consider differences in the dynamics of intramolecular domains between bound and unbound states of the 3KPS system. In the second and third set of configurations, simulations were started from isolated, i.e., unbound TCR-and pMHC-molecules as derived from the original 3KPS crystal structure.  [22].

Molecular Dynamics Simulations
MD simulations were run in an explicit water solvent using GROMACS 5.1.1 [23]. We used the amber99sb-ildn force field [24] to simulate the protein placed in a rhombic dodecahedral box with 2 nm minimum distance between the protein atoms and the box boundary. The whole system contains approximately 397,000 atoms, consisting of 826 protein residues and approximately 128,000 solvent molecules. Based on the SPC water model [25] the protein was solvated and neutralized by replacing solvent molecules with Na + and Cl − ions to reach a salt concentration of 0.15 mol/L.
After solvation the whole complex was energetically minimized using the steepest-descent approach with 10,000 steps. For equilibration we set the temperature of the system to 310 K and the equations of motion were followed for 100 ps using position restraint MD and a Berendsenthermostat with a time constant of 0.1 ps. This equilibration under NVT conditions (constant number of particles, volume and temperature) was followed by a 100 ps equilibration in an NPT ensemble (constant number of particles, pressure and temperature) to control the pressure by a Berendsenbarostat set to 1 bar with a time constant of 1.0 ps.
The MD production runs were carried out with a time-step of 4 fs using the LINCS constraint algorithm for all bonds and virtual sites for hydrogen atoms. For the Van der Walls and Coulomb interactions a cut-off radius of 1.4 nm was applied. Long-range Coulomb interactions were computed using the particle-mesh Ewald method with a maximum grid spacing of 0.12 nm for the FFT and an

Molecular Dynamics Simulations
MD simulations were run in an explicit water solvent using GROMACS 5.1.1 [23]. We used the amber99sb-ildn force field [24] to simulate the protein placed in a rhombic dodecahedral box with 2 nm minimum distance between the protein atoms and the box boundary. The whole system contains approximately 397,000 atoms, consisting of 826 protein residues and approximately 128,000 solvent molecules. Based on the SPC water model [25] the protein was solvated and neutralized by replacing solvent molecules with Na + and Cl − ions to reach a salt concentration of 0.15 mol/L. After solvation the whole complex was energetically minimized using the steepest-descent approach with 10,000 steps. For equilibration we set the temperature of the system to 310 K and the equations of motion were followed for 100 ps using position restraint MD and a Berendsen-thermostat with a time constant of 0.1 ps. This equilibration under NVT conditions (constant number of particles, volume and temperature) was followed by a 100 ps equilibration in an NPT ensemble (constant number of particles, pressure and temperature) to control the pressure by a Berendsen-barostat set to 1 bar with a time constant of 1.0 ps.
The MD production runs were carried out with a time-step of 4 fs using the LINCS constraint algorithm for all bonds and virtual sites for hydrogen atoms. For the Van der Walls and Coulomb interactions a cut-off radius of 1.4 nm was applied. Long-range Coulomb interactions were computed using the particle-mesh Ewald method with a maximum grid spacing of 0.12 nm for the FFT and an interpolation order of 4. Velocity-rescaling temperature coupling was set to 310 K and Berendsen isotropic pressure coupling was set to 1 bar. Coordinates were written to the trajectory file on hard disk every 40 ps, resulting in 1000 frames for each run. This procedure was repeated for each run in the set of initial configurations. We performed 30 independent MD simulation runs with different initial velocities and with a length of 40 ns each: 10 runs for the TCR-pMHC complex (bound configurations), 10 runs for a free pMHC molecule, and 10 runs for a free TCR molecule (unbound configurations). For reasons of comparison, two additional runs of 400 ns in length (one for the free pMHC-and the other for the free TCR-molecule) and one 144 ns run for the TCR-pMHC complex were performed. Prior to the analysis of the trajectories, translational and rotational motions of the protein relative to the equilibrated structure were removed.

Distance between Domains
One commonly used method among many others [26][27][28] to characterize relative movements of intramolecular domains is the calculation of distances between groups of atoms (domains) as a function of simulation time. We evaluated the distances d between the geometric centers of domains, based on the C α atoms in the protein backbone.

Relative Orientation of Two Intramolecular Domains
Based on the methods described previously [29] we monitored the relative orientations of intramolecular domains during the time-course of the simulations. For the definition of domains we used the C α atoms of the protein backbone. As a first pair of domains we chose the two epitope binding domains TCR V α (104 residues) and TCR V β (117 residues), both involved in the formation of the interface with pMHC. As a second pair we considered the antigen presenting sites of MHC, namely MHC α1 (25 residues) and MHC α2 (31 residues).
During an MD simulation the atoms within each domain move according to the equations of motion and thus their relative orientation changes from time step to time step. Principle component analysis (PCA) is commonly used to analyze the motions of domains within the conformations obtained during MD simulations [26]. PCA as such yields eigenvectors unique in directions but ambiguous in orientation. To arrive at standardized eigenvectors allowing for comparison between different MD-runs we chose the crystal structure as a reference frame and re-oriented (if necessary) the resulting eigenvectors after performing the PCA [29]. The three orthonormal eigenvectors from the PCA were used to define an eigenvector matrix to serve as a local coordinate system, moving together with the domain, see Figure 2. For each frame of an MD simulation we monitored the relative orientation between local coordinate systems of the domains.
To actually characterize the relative orientations of domains and their changes over time, we computed the cosines between the corresponding eigenvectors pairs (v 1 , w 1 ), (v 2 , w 2 ) and (v 3 , w 3 ) of the respective local coordinate systems.

Domain Deformations and Fluctuations
To quantify intra-domain deformations we computed the RMSD (root-mean-square deviation) of each frame with respect to the first frame of each trajectory, considering only the Cα atoms within the respective domain. The root mean square deviation RMSD(t) of a group of atoms in a molecule at time t with respect to a given reference structure is a measure of deformations (deviations in shape), i.e., the square root of averaged squared distances between the atom positions at time t and the positions in the reference structure. Here, the starting structure of the production runs (i.e., the first frame of the trajectory) was chosen as the reference structure. As such, RMSD(t) is a statistical measure of the distance of a group of atoms at a particular time t with respect to the same group in the reference structure and was calculated as: where ri(t) is the position of atom i at time t after least square fitting to the reference structure, N is the total number of atoms in the group of atoms considered, and ri ref is the position of atom i in the reference structure. On the other hand, the root mean square fluctuation RMSF(i) is a statistical measure of the deviation between the position of atom i (or a group of atoms, e.g., a residue) and some reference position ri ref : where T is the time interval over which the average is taken and tj is the time of frame j during the simulation of the respective trajectory. Here, we have chosen the time-averaged position of atom i as

Domain Deformations and Fluctuations
To quantify intra-domain deformations we computed the RMSD (root-mean-square deviation) of each frame with respect to the first frame of each trajectory, considering only the C α atoms within the respective domain. The root mean square deviation RMSD(t) of a group of atoms in a molecule at time t with respect to a given reference structure is a measure of deformations (deviations in shape), i.e., the square root of averaged squared distances between the atom positions at time t and the positions in the reference structure. Here, the starting structure of the production runs (i.e., the first frame of the trajectory) was chosen as the reference structure. As such, RMSD(t) is a statistical measure of the distance of a group of atoms at a particular time t with respect to the same group in the reference structure and was calculated as: where r i (t) is the position of atom i at time t after least square fitting to the reference structure, N is the total number of atoms in the group of atoms considered, and r i ref is the position of atom i in the reference structure. On the other hand, the root mean square fluctuation RMSF(i) is a statistical measure of the deviation between the position of atom i (or a group of atoms, e.g., a residue) and some reference position r i ref : where T is the time interval over which the average is taken and t j is the time of frame j during the simulation of the respective trajectory. Here, we have chosen the time-averaged position of atom i as the reference position, i.e., r i ref = <r i >. Whereas the RMSD(t) is an average taken over all atoms within a group for a specific time t, the RMSF(i) is an average over time for a specific atom i.

Relative Movements between TCR V α and TCR V β
The variable regions TCR V α and TCR V β interact with the pMHC surface and are therefore important during TCR-pMHC recognition. Interestingly, if TCR and pMHC are free (i.e., unbound), the distance d between the domains TCR V α and TCR V β is not distinctively different from the bound configurations (results from ten 40 ns runs): < d free > = 2.523 nm, < d bound > = 2.537 nm, ∆< d > = 0.014 nm (difference of mean distances < d > between bound and unbound configurations), σ free = 0.005 nm (standard deviation of distances for unbound configurations), and σ bound = 0.002 nm (standard deviation of distances for bound configurations), see Figure 3. the reference position, i.e., ri ref = <ri>. Whereas the RMSD(t) is an average taken over all atoms within a group for a specific time t, the RMSF(i) is an average over time for a specific atom i.

Relative Movements between TCR Vα and TCR Vβ
The variable regions TCR Vα and TCR Vβ interact with the pMHC surface and are therefore important during TCR-pMHC recognition. Interestingly, if TCR and pMHC are free (i.e., unbound), the distance d between the domains TCR Vα and TCR Vβ is not distinctively different from the bound configurations (results from ten 40 ns runs): < dfree > = 2.523 nm, < dbound > = 2.537 nm, ∆< d > = 0.014 nm (difference of mean distances < d > between bound and unbound configurations), σfree = 0.005 nm (standard deviation of distances for unbound configurations), and σbound = 0.002 nm (standard deviation of distances for bound configurations), see Figure 3. Deformations within each domain can be characterized by RMSD-values. Figure 4 shows boxplots of RMSD-values computed with respect to the first frame of each trajectory for the unbound (runs 1 to 10) and bound states (runs 11 to 20). Deformations within TCR Vα and TCR Vβ. are virtually unaffected by the binding of the pMHC und TCR molecules (see Figure 4, Tables 1 and 2): For TCR Vα the overall mean ± SD of RMSD-values from ten 40 ns runs were 0.115 ± 0.002 nm for the unbound states and 0.101 ± 0.002 nm for the bound states. For TCR Vβ the respective values were 0.096 ± 0.003 nm (unbound) and 0.101 ± 0.002 nm (bound). Deformations within each domain can be characterized by RMSD-values. Figure 4 shows boxplots of RMSD-values computed with respect to the first frame of each trajectory for the unbound (runs 1 to 10) and bound states (runs 11 to 20). Deformations within TCR V α and TCR V β are virtually unaffected by the binding of the pMHC und TCR molecules (see Figure 4, Tables 1 and 2): For TCR V α the overall mean ± SD of RMSD-values from ten 40 ns runs were 0.115 ± 0.002 nm for the unbound states and 0.101 ± 0.002 nm for the bound states. For TCR V β the respective values were 0.096 ± 0.003 nm (unbound) and 0.101 ± 0.002 nm (bound).      Table 2. RMSD (mean ± SD) in nm for various domains of the bound complex TCR-pMHC as calculated from ten 40 ns MD runs (1-10) and from one 144 ns MD run (11). We characterized the relative orientation of the domains TCR V α and TCR V β by the directional cosines between eigenvectors attached to each domain and sampled over time for all trajectories as outlined in the methods section and in our previous paper [29], see Figure 5. The eigenvectors v 1 and w 1 correspond to the main extensions of the respective domain. To quantify the differences between bound and unbound states in the frequency distributions of the respective angles α between v 1 and w 1 we used the test statistic D of the Kolmogorov Smirnov test. For v 1 and w 1 the results from ten 40 ns runs were Larger differences were observed for the angles β between v 2 and w 2 : These results may be summarized as follows: The domains TCR V α and TCR V β remain relatively unchanged in their main orientation after complexation of TCR and pMHC, as seen by the small changes of the angles between v 1 and w 1 (58.96 • on average for the unbound state versus 59.88 • in the bound state), see Figure 5A. An additional change in orientation is due to a rotation around these major axes as reflected in the variation of the angles β between v 2 and w 2 . As a result of complexation their mean values change from 57.38 • for the unbound states to 62.78 • for the bound states, see Figure 5B. w1 correspond to the main extensions of the respective domain. To quantify the differences between bound and unbound states in the frequency distributions of the respective angles α between v1 and w1 we used the test statistic D of the Kolmogorov Smirnov test. For v1 and w1 the results from ten 40 ns runs were D = 0.215, < α(v1,w1) >free = 58.96°, σfree = 2.15°, < α(v1,w1) >bound = 59.88°, σbound = 1.75°. Larger differences were observed for the angles β between v2 and w2: D = 0.643, < β(v2,w2) >free = 57.38°, σfree = 2.36°, < β(v2,w2) >bound = 62.78°, σbound = 3.78°. These results may be summarized as follows: The domains TCR Vα and TCR Vβ remain relatively unchanged in their main orientation after complexation of TCR and pMHC, as seen by the small changes of the angles between v1 and w1 (58.96° on average for the unbound state versus 59.88° in the bound state), see Figure 5A. An additional change in orientation is due to a rotation around these major axes as reflected in the variation of the angles β between v2 and w2. As a result of complexation their mean values change from 57.38° for the unbound states to 62.78° for the bound states, see Figure  5B.

Relative Movements between MHC α1 and MHC α2
The two α-helices of the MHC, together with a β-floor, form the peptide binding cleft that presents antigen fragments to be examined by the TCR. Together with the variable domains of the TCR, TCR Vα and TCR Vβ, they are most important for the formation of the TCR-pMHC complex and for a successful initiation of the immune response. We thus also monitored intramolecular movements of the two α-helices in the free (i.e., unbound) and bound states. If TCR and pMHC are unbound, the distance d between domains MHC α1 and MHC α2 practically does not change with

Relative Movements between MHC α1 and MHC α2
The two α-helices of the MHC, together with a β-floor, form the peptide binding cleft that presents antigen fragments to be examined by the TCR. Together with the variable domains of the TCR, TCR V α and TCR V β , they are most important for the formation of the TCR-pMHC complex and for a successful initiation of the immune response. We thus also monitored intramolecular movements of the two α-helices in the free (i.e., unbound) and bound states. If TCR and pMHC are unbound, the distance d between domains MHC α1 and MHC α2 practically does not change with respect to the bound state (results from ten 40 ns runs): < d free > = 1.676 nm, < d bound > = 1.649 nm, ∆< d > = 0.028 nm (difference of mean distances < d > between unbound and bound configurations), σ free = 0.004 nm (standard deviation of distances for unbound configurations), and σ bound = 0.004 nm (standard deviation of distances for bound configurations), see Figure 6.  To monitor the relative orientation between the helical domains MHC α1 and MHC α2 we calculated directional cosines between corresponding eigenvectors and examined their distributions over time, see Figure 8. The angles α between the eigenvectors v1 and w1 (corresponding to the main extensions of the domains) were slightly different between unbound and bound states (results from ten 40 ns runs): D = 0.255, < α(v1,w1) >free = 20.75°, σfree = 2.74°, < α(v1,w1) >bound = 22.50°, σbound = 1.85°. More pronounced differences and larger fluctuations were observed for the angles β between v2 and w2: D = 0.290, < β(v2,w2) >free = 86.57°, σfree = 12.26°, < β(v2,w2) >bound = 81.16°, σbound = 11.06°. These eigenvectors v2 and w2 point away from the axis of the helices at right angles and indicate a larger variation in the relative orientation between MHC α1 and MHC α2 in the unbound state as compared to the bound state. To monitor the relative orientation between the helical domains MHC α1 and MHC α2 we calculated directional cosines between corresponding eigenvectors and examined their distributions over time, see Figure 8. The angles α between the eigenvectors v 1 and w 1 (corresponding to the main extensions of the domains) were slightly different between unbound and bound states (results from ten 40 ns runs): D = 0.255, < α(v 1 ,w 1 ) > free = 20.75 • , σ free = 2.74 • , < α(v 1 ,w 1 ) > bound = 22.50 • , σ bound = 1.85 • . More pronounced differences and larger fluctuations were observed for the angles β between v 2 and w 2 : D = 0.290, < β(v 2 ,w 2 ) > free = 86.57 • , σ free = 12.26 • , < β(v 2 ,w 2 ) > bound = 81.16 • , σ bound = 11.06 • . These eigenvectors v 2 and w 2 point away from the axis of the helices at right angles and indicate a larger variation in the relative orientation between MHC α1 and MHC α2 in the unbound state as compared to the bound state. To monitor the relative orientation between the helical domains MHC α1 and MHC α2 we calculated directional cosines between corresponding eigenvectors and examined their distributions over time, see Figure 8. The angles α between the eigenvectors v1 and w1 (corresponding to the main extensions of the domains) were slightly different between unbound and bound states (results from ten 40 ns runs): D = 0.255, < α(v1,w1) >free = 20.75°, σfree = 2.74°, < α(v1,w1) >bound = 22.50°, σbound = 1.85°. More pronounced differences and larger fluctuations were observed for the angles β between v2 and w2: D = 0.290, < β(v2,w2) >free = 86.57°, σfree = 12.26°, < β(v2,w2) >bound = 81.16°, σbound = 11.06°. These eigenvectors v2 and w2 point away from the axis of the helices at right angles and indicate a larger variation in the relative orientation between MHC α1 and MHC α2 in the unbound state as compared to the bound state. (A)

Discussion
Crystallographic studies of TCR-pMHC complexes have delivered important insights regarding the structural components involved in stabilizing the TCR-pMHC complex. Despite the huge success and the valuable contributions of these studies, they provide only a limited static view. On the other hand, molecular dynamics simulations allow for following the dynamics of single atoms or groups of atoms with a time-resolution in the order of pico-seconds. Although the interaction between a TCR and a pMHC is weak (10 3 -10 5 M −1 s −1 [30]), it represents an important initial step in T-cell activation. Detailed studies on protein-protein complexes suggest that the binding process of protein molecules can be divided into three phases [31]: A first and short phase of diffusion (initial contacts form within 2 ns) is followed by an intermediate phase where most native contacts are established (between 50-200 ns). In the last phase the final stereospecific complex is formed (this takes hundreds of nanoseconds). Various studies have reported conformational changes at the TCR-pMHC binding interface in the course of complex formation, in particular at the CDRs and peptides and less at the MHC [32][33][34]. In the current work, we present molecular dynamics simulations to study the differences in intramolecular domain movements between bound and unbound states of a specific TCR-pMHC complex.
For the analyses of domain orientations, we attached local coordinate systems (PCAeigenvectors) to the domains and calculated the cosines between corresponding eigenvectors. Based on three sets of independent simulations with different initial configurations (TCR-pMHC bound, TCR unbound, pMHC unbound) we could show that binding of TCR to pMHC decreases the flexibility in the inter-domain movements, in particular the dynamics of TCR Vα relative to TCR Vβ and of MHC α1 relative to MHC α2, see the distribution of angles between corresponding principal component eigenvectors in Figure 5 and Figure 8.
Further examination of domain RMSD-values revealed that MHC α2-dynamics clearly depends on the binding states (i.e., bound versus unbound) of TCR and pMHC molecules, see Figure 9. Conformational flexibility, as characterized by the RMSD-values of MHC α2, was distinctively larger in the unbound states as compared to the bound states. On the contrary, differences between bound and unbound states in RMSD-values of the TCR domains Vα and Vβ were considerably less pronounced, see Figure 10. The decrease of domain RMSD-values of MHC α2 is consistent with the

Discussion
Crystallographic studies of TCR-pMHC complexes have delivered important insights regarding the structural components involved in stabilizing the TCR-pMHC complex. Despite the huge success and the valuable contributions of these studies, they provide only a limited static view. On the other hand, molecular dynamics simulations allow for following the dynamics of single atoms or groups of atoms with a time-resolution in the order of pico-seconds. Although the interaction between a TCR and a pMHC is weak (10 3 -10 5 M −1 s −1 [30]), it represents an important initial step in T-cell activation. Detailed studies on protein-protein complexes suggest that the binding process of protein molecules can be divided into three phases [31]: A first and short phase of diffusion (initial contacts form within 2 ns) is followed by an intermediate phase where most native contacts are established (between 50-200 ns). In the last phase the final stereospecific complex is formed (this takes hundreds of nanoseconds). Various studies have reported conformational changes at the TCR-pMHC binding interface in the course of complex formation, in particular at the CDRs and peptides and less at the MHC [32][33][34]. In the current work, we present molecular dynamics simulations to study the differences in intramolecular domain movements between bound and unbound states of a specific TCR-pMHC complex.
For the analyses of domain orientations, we attached local coordinate systems (PCA-eigenvectors) to the domains and calculated the cosines between corresponding eigenvectors. Based on three sets of independent simulations with different initial configurations (TCR-pMHC bound, TCR unbound, pMHC unbound) we could show that binding of TCR to pMHC decreases the flexibility in the inter-domain movements, in particular the dynamics of TCR V α relative to TCR V β and of MHC α1 relative to MHC α2, see the distribution of angles between corresponding principal component eigenvectors in Figures 5 and 8.
Further examination of domain RMSD-values revealed that MHC α2-dynamics clearly depends on the binding states (i.e., bound versus unbound) of TCR and pMHC molecules, see Figure 9. Conformational flexibility, as characterized by the RMSD-values of MHC α2, was distinctively larger in the unbound states as compared to the bound states. On the contrary, differences between bound and unbound states in RMSD-values of the TCR domains V α and V β were considerably less pronounced, see Figure 10. The decrease of domain RMSD-values of MHC α2 is consistent with the dynamics of the peptide placed in the binding groove between the two MHC α-helices: Root mean square fluctuations (RMSF) of individual residues of the peptide during the course of the simulations show distinctively larger values in the unbound as compared to the bound states ( Figure 11). In particular, central amino acid residues at positions 4-7 within the peptide (Leu, Gln, Ala, Phe) display reduced dynamic flexibility upon binding of the pMHC with the TCR, in agreement with previous studies of various peptide/MHC-I complexes and TCRs [35,36], see [6,[37][38][39] for reviews. On the other hand, the N-and C-termini of positions 1 and 9 in the peptide are much less influenced by TCR-binding due to the network of hydrogen bonds from the β-sheet and from the lateral α-helices which stabilize the peptide binding domain [12]. Macdonald et al. [18] emphasized the importance of the C-terminal region (p6-p8) of the peptide for TCR recognition: In particular, the small p6 (Ala) and p8 (Thr) residues enabled the aromatic p7 (Phe) residue to protrude within a central pocket of the TCR, sandwiched between Tyr31α and Tyr100β; the p6-Ala made interactions with the CDR3α and the CDR3β loop, while p8-Thr and p7-Phe formed contacts with the CDR3β loop.
Residues in the variable domains V α and V β of the T-cell receptor show less pronounced differences in RMSF-values between free and bound states, with the exception of residues 30 (Thr), 31 (Tyr), 58 (Arg), 96-98 (Gly, Gly, Thr), 100 (Tyr), and 102 (Gly) in the TCR V α domain and residues 50-52 (Gln, Asn, Glu), 71-73 (Glu, Gly, Ser), 98-100 (Gln, Ala, Tyr), and 105 (Glu) in the TCR V β domain ( Figure 12). Apart from residue 58 in the V α domain and residues 71-73 in the V β domain, all the other residues with a distinctively larger RMSF-value in the unbound states are part of the CDR-loops CDR1α (30,31), CDR3α (96-98, 100, 102), CDR2β (50-52), and CDR3β (98-100, 105). These results are in agreement with several studies that have found a reduced mobility of the TCR [6,40] and conformational changes in the CDR-loops [4,7,41] upon formation of the TCR-pMHC complex. In particular, we observed a structural change between the free and bound states in the CDR3α loop of residues 95-97 (Ala, Gly, Gly) from a 3 10 -helix (free) to a coil (bound), where Gly96α and Gly97α form interactions with Arg62 and Ile66 of the MHC α1 helix, respectively [18] (see supplementary Figure  S1). As to the type of interactions between CDR-loops and pMHC, Macdonald et al. [18] reported predominantly van der Waals interactions and H-bonds, including one salt bridge between Glu52β of CDRβ2 and Arg79 of the α1 helix. These crystallographic studies have also shown that CDR1α and CDR2α loops make contacts with the α2 helix of pMHC, whereas contacts to the α1 helix are established by the CDR2β and CDR3α loops (CDR3α also makes interactions with the α2 helix and the peptide). While the CDR1β loop participates only minimally in the TCR-pMHC interactions, the contacts with the peptide are dominated by the CDR3β loop, which sits centrally above the peptide-binding cleft, adjoining the CDR3α loop (see also Figure 12C).  Figure 11). In particular, central amino acid residues at positions 4-7 within the peptide (Leu, Gln, Ala, Phe) display reduced dynamic flexibility upon binding of the pMHC with the TCR, in agreement with previous studies of various peptide/MHC-I complexes and TCRs [35,36], see [6,[37][38][39] for reviews. On the other hand, the N-and C-termini of positions 1 and 9 in the peptide are much less influenced by TCRbinding due to the network of hydrogen bonds from the β-sheet and from the lateral α-helices which stabilize the peptide binding domain [12]. Macdonald et al. [18] emphasized the importance of the Cterminal region (p6-p8) of the peptide for TCR recognition: In particular, the small p6 (Ala) and p8 (Thr) residues enabled the aromatic p7 ( (98-100, 105). These results are in agreement with several studies that have found a reduced mobility of the TCR [6,40] and conformational changes in the CDR-loops [4,7,41] upon formation of the TCR-pMHC complex. In particular, we observed a structural change between the free and bound states in the CDR3α loop of residues 95-97 (Ala, Gly, Gly) from a 310-helix (free) to a coil (bound), where Gly96α and Gly97α form interactions with Arg62 and Ile66 of the MHC α1 helix, respectively [18] (see supplementary Figure S1). As to the type of interactions between CDR-loops and pMHC, Macdonald et al. [18] reported predominantly van der Waals interactions and H-bonds, including one salt bridge between Glu52β of CDRβ2 and Arg79 of the α1 helix. These crystallographic studies have also shown that CDR1α and CDR2α loops make contacts with the α2 helix of pMHC, whereas contacts to the α1 helix are established by the CDR2β and CDR3α loops (CDR3α also makes interactions with the α2 helix and the peptide). While the CDR1β loop participates only minimally in the TCR-pMHC interactions, the contacts with the peptide are dominated by the CDR3β loop, which sits centrally above the peptide-binding cleft, adjoining the CDR3α loop (see also Figure 12C).
(A)   Only minor changes were observed in the inter-domain distance between TCR V α relative to TCR V β and MHC α1 relative to MHC α2 as a consequence of binding. Fluctuations within domains (TCR V α , TCR V β , MHC α1, and MHC α2) are virtually unaffected by the binding state, see Figures 4 and 7. In addition, angles between corresponding eigenvectors of MHC α1 and MHC α2 show larger variation in the unbound states, indicating a wider range of the groove between the two α-helices, see Figures 5 and 8.
A study based on static crystallographic data [42] demonstrated that the angle between TCR V α and TCR V β domains varies depending on the presence of peptides, whereby domains show larger variations in the unbound state and smaller variations in the bound state of the TCR. These findings are consistent with our MD-based results regarding the angles between vectors v 1 and w 1 .
As to the design of our MD experiments, we performed 10 independent 40 ns runs for each of the 3 configurations studied (pMHC free, TCR free, TCR-pMHC complexed). This relatively short simulation length was chosen because it is more efficient (in terms of sampling the conformational space) to run various short independent MD simulations than a single long one [43], see, e.g., the variations of RMSD-values between the runs in the boxplots of Figure 7. On the other hand, RMSD mean values as calculated from ten pooled 40 ns runs were consistent with the respective values obtained from one 144 ns run (TCR-pMHC complex) und two 400 ns runs (TCR and pMHC free), see run 11 in Tables 1  and 2.
There are several limitations in the present study, in particular the influence of the lipid bilayer on TCR-pMHC interactions. Since the MHC and TCR are membrane-embedded proteins, the influence of the membrane environment itself is an essential ingredient for TCR-pMHC interaction and T cell response [44]. Various advanced MD studies have been performed to evaluate the role of the lipid bilayer for TCR-pMHCII systems: In a ground-breaking study, Wan et al. [45] performed MD simulations of a membrane-embedded TCR-pMHC-CD4 complex. Although the duration of the simulation was limited to 10 ns due to the size of the system (more than 300.000 atoms), the computed structural and thermodynamic properties, such as binding free energies of TCR-pMHC and pMHC-CD4, were in fair agreement with experimental data. Bello et al. [46] studied the influence of the membrane on the dynamics and energetics of peptide-bound and peptide-free MHC class II molecules in aqueous and membrane-bound environments by 100 and 150 ns MD simulations and found that the presence of the membrane might restrict the conformational flexibility of the peptide-binding cleft. More recently, Bello et al. [47] explored the energetic and dynamic behavior of a pMHCII-TCR complex embedded in two opposing lipid membranes using three independent 300 ns-long unbiased MD simulations. This study provided evidence of the main contributors to the pMHCII-TCR complex formation and identified key residues involved in this molecular recognition process.
The formation of the immunological synapse incorporates co-receptors, integrins, and various signaling proteins that trigger signaling cascades resulting in T-cell activation. Therefore, future MD studies of the 3KPS system should not only include the lipid membrane, but also co-receptors (CD8, CD3 complex) and other co-stimulatory molecules to assess their influence on the system [6].
Another limitation of the present study relates to the simulation time of 40 ns, which is way too short to reliably cover all processes involved, e.g., in TCR-signaling or even dissociation of the TCR-pMHC complex. Several authors have used coarse-grained simulations and elastic network models for modeling protein flexibility and interactions, see the recent reviews by Kmiecik et al. [48,49]. Moreover, steered MD [50] has been applied to extend simulation times to more realistic regimes.

Conclusions
We have performed MD simulations with different starting configurations to monitor differences in intramolecular domain movements between bound and unbound TCR-pMHC molecules of the LC13-HLA-B*44:05-pEEYLQAFTY complex. The results provide some insights into alterations of the dynamics as a consequence of complex formation: For the domains MHC α1 and MHC α2 the angles between corresponding eigenvectors show larger variation in the unbound state, indicating a wider range of the groove between the two α-helices. For the TCR V α and TCR V β domains, angles between corresponding eigenvectors indicate little movements in their main orientation but significant rotations around these major axes. The MHC α2-helix showed larger RMSD flexibility for unbound states. The observed structural changes in the CDR3α loop upon binding together with the large RMSF-values of the same region in the unbound state (see Figure 12A) are consistent with the view that structural changes involved in protein binding correlate with intrinsic motions of proteins in the unbound state [51]. Moreover, the structural transition from a 3 10 -helix to a coil in the CDR3α loop upon complex formation supports the hypothesis that (at least for the CDR3α loop) TCR-pMHC binding can be characterized by the induced fit model, whereby structural rearrangements of the TCR CDRs ensure to achieve the most stable complex [3]. Differences in conformational flexibility (as characterized by RMSD-and RMSF-values) between free and bound forms both for the MHC α-helices and the peptide are consistent with the view that the TCR perceives the peptide and the MHC as a single, composite ligand, making only little distinction between them [7]. Although even our 30 relatively short MD simulations of 40 ns each revealed interesting dynamical features not available from crystallographic data alone, the length of simulations should be extended in subsequent studies to cover processes with longer characteristic time scales.
Author Contributions: R.K. performed MD simulations, implemented key parts of the algorithms, performed evaluations and wrote parts of the manuscript. C.S. performed MD simulations, evaluations and wrote parts of the manuscript. N.I. performed initial MD simulations and contributed concepts in preparation of the project. W.S. contributed concepts for evaluation and coordinated the study.
Data Availability: For this study data have been downloaded from the publically available protein data bank (https://www.rcsb.org/). Funding: Work was carried out by staff of the Medical University of Vienna, no external funding was used.