Coordinates of B*44:02 loaded with the HLA DPA*0201 9-mer peptide were taken from a 1.6-Å resolution X-ray crystal structure (PDB ID 1M6O) (40). The α subunit contains 276 residues and excludes the membrane-spanning helix and cytosolic tail, while the invariant β₂m subunit contains 99 residues. Water molecules resolved in the crystal structure were retained. PROPKA 3.1 (41) was used to determine protonation states at pH 7. D156 in the α₂₋₂ helix segment was the only residue with a non-standard protonation state, sharing a proton with D114 in β₆ at the bottom of the binding groove.

The structure was first subjected to 60 steps of steepest-descent (SD) energy minimization. Peptide removal, truncation, and side chain replacement were then used to generate from MHCI^PL the starting structures for MHCI^PD, MHCI loaded with six truncated DPA peptides (Δ1N, Δ2N, Δ3N, Δ1C, Δ2C, Δ3C) and MHCI loaded with three peptides in which one or two residues were replaced by glycines (δ1N, δ2N, δ1C); all peptides were capped with standard (charged) amine and carboxyl termini. All systems were solvated in periodic rhombic dodecahedron cells, with a 10-Å minimal distance between the solute and the box edge. Random water molecules were replaced by Na⁺ and Cl⁻ ions at a 0.15 M concentration, yielding a neutral net charge. The systems contained ≈69,000 atoms in cells of ≈280,000 ų. Finally, all systems were again subjected to 500 steps of SD energy minimization.

The tapasin–MHCI^Δ1C system was built from the results of a previous simulation of tapasin–MHCI (14, 15) in a solvated periodic rhombic dodecahedron cell with a 10-Å minimal distance between the solute and the box edge, and a 0.15 M concentration of Na⁺ and Cl⁻ ions. The system contained ≈140,000 atoms in a ≈1,480,000 ų cell. The DPA peptide in this system was truncated and capped with a carboxyl terminus prior to 500 steps of SD energy minimization.

Simulations were carried out with GROMACS 2016.1 (42, 43). The Amber99SB*-ILDNP protein forcefield (44–49) and TIP3P water model (50) were used. The SETTLE (51) and LINCS (52) algorithms were applied to constrain the internal degrees of freedom of water molecules and the bonds in other molecules, respectively. In combination with virtual site hydrogens (53), this allowed for a 4-fs integration time step. Short-range non-bonded Coulomb and Lennard-Jones 6–12 interactions were treated with a Verlet-buffered pair list (54) with potentials smoothly shifted to zero at a 10-Å cutoff. Long-range Coulomb interactions were treated with the PME method (55) with a grid spacing of 1.2 Å and cubic spline interpolation. Analytical dispersion corrections were applied for energy and pressure to compensate for the truncation of the Lennard-Jones interactions. The thermodynamic ensemble was nPT. Temperature was kept constant at 300 K by a velocity-rescaling thermostat with a stochastic term (56), with coupling time constant 0.1 ps. For constant 1.0-bar pressure, an isotropic Berendsen barostat (57) was used with coupling time constant 0.5 ps and 4.5 × 10⁻⁵/bar compressibility.

For each system, five independent 1.0-µs trajectories were acquired (by generating random initial velocities from a Boltzmann distribution), with coordinates recorded every 10 ps. All of these simulations were preceded by 10 ns of equilibration, during which all heavy (non-hydrogen) protein atoms were position-restrained by harmonic potential energy functions with force constants of 1,000 kJ/(mol nm²).

The first 10 ns of simulation without position restraints were considered further equilibration time and were not used for analysis of averaged RMSD, binding groove width, RMSF, configurational entropy, and pairwise forces between residues. In addition, one trajectory involving the Δ1C peptide was excluded from all analyses since the peptide dissociated entirely from the binding groove (see Results for further details).

Configurational entropies were calculated using the quasi-harmonic approximation (QHA) as formulated by Schlitter (58), (1)S=0.5 kB ln  det(1+kBTe2ℏ−2 M1/2 C M1/2)>Strue, which provides an upper bound to the configurational entropy. In this equation, k_B is the Boltzmann constant, T the temperature, e Euler’s number, ℏ the reduced Planck constant, and M the 3N-dimensional diagonal mass matrix of the N particles. The matrix C is the covariance matrix of particle positions, (2)C=〈(x−〈x〉)(x−〈x〉)T〉, where the 3N-dimensional vector x represents the Cartesian coordinates of the N particles for which the entropy is calculated after removing overall translation and rotation by fitting to a reference structure. The starting structure of our simulations was used as the reference structure for this fit. The coordinates of the C_α atoms were used to construct C. By comparing the entropy accumulated for the full 5.0 µs of sampling for each system to the value obtained when discarding a random 10% of the data, we estimated the error and verified that our sampling enables us to evaluate configurational entropies with sufficient statistical reliability.

Force distribution analysis (FDA) was performed using the GROMACS-PF2 implementation (59). Pairwise atom–atom forces were recomputed from the trajectories and were time- and residue-averaged. For the latter, the norm of the sum of the interatomic force vectors between all atoms i and j of residues v and w were calculated, Fvwres=|∑ij Fij|, with i ∈ v and j ∈ w. Attractive and repulsive forces were distinguished by sign.

Figure 2 shows the C_α RMSD to the X-ray crystal structure for all studied peptides. Deletion of a single C-terminal residue (Δ1C) is sufficient to trigger partial peptide dissociation from the binding groove on the simulation timescale. The C-terminal F9 residue, whose side chain is buried in the F-pocket, is thus crucial for binding to B*44:02. In its absence, central residues cannot stabilize the peptide, and a large segment (G4–S8) dissociates from the MHCI binding groove (Figure 2, Δ1C). Further deletions (Δ2C, Δ3C) have no additional effect; the same G4–S8 region unbinds from the groove. These results are consistent with the observation (36) that removing or modifying the side chain of the peptide C-terminal residue accelerates dipeptide-catalyzed peptide exchange.

In contrast to the C-terminus, a single N-terminal deletion (Δ1N) has a less pronounced effect and leads only to the partial dissociation of residues E2–F3 (Figure 2). Removal of E2, which is buried in the A-pocket, is necessary for dissociation of the F3–R5 region, comparable to the loss of interactions observed for C-terminal deletions. Thus, the peptide is held in place by its termini at E2 and F9, while the residues in between contribute less to MHCI binding. The importance of peptide termini was shown previously from the thermal stability of MHCI molecules loaded with modified peptides (60). Our results are also consistent with the conclusions of a previous MD study of MHCI allele H2-K^b loaded with truncated peptides (37), although the shorter simulation times of that study precluded the observation of partial peptide dissociation.

Due to the thermal fluctuations, even the full-length DPA peptide deviates from its position in the X-ray crystal structure. While the N- and C-terminal residues remain firmly anchored to the A- and F-pocket, respectively, central residues G4, R5, and A6 exhibit an average C_α RMSD of about 1.5 Å, consistent with the idea that residues E2 and F9 are the main determinants of peptide–MHCI stability. In line with this observation, even though the dehydrated environment of the crystal might naturally favor a tight binding groove, crystallographic Debye–Waller factors (B-factors) are slightly higher for residues 3–8 than they are for E1, E2, and F9. This peptide flexibility, observed in our simulations of MHCI in solution, could have implications for antigen recognition. MHC recognition by TCR, reviewed recently (61), involves three TCR loops (CDR1, 2, and 3); while CDR1 and 2 contact MHC, CDR3 is located directly atop the antigen peptide and thus interacts mostly with residues at the center of the peptide. For instance, interaction hotspots at positions 4, 5, and 8 in a 9-mer peptide have been identified in the structure of a TCR complexed with an MHCI molecule loaded with a cytomegalovirus antigen (62). Conformational changes observed in CDR3 upon binding to MHCI^PL (61) were interpreted in terms of an induced fit mechanism that could be facilitated by small-scale rearrangements and fluctuations of the peptide central residues.

In all but one of the simulations performed in the present study, the peptide remained partially bound to B*44:02 and did not fully dissociate from the binding groove. However, in one of the Δ1C peptide simulations, the peptide C-terminus dissociated from the groove (after 5 ns), before the N-terminus also dissociated (after 255 ns). This trajectory was excluded from all further analyses. Investigating the possibility of complete peptide dissociation would require much more extensive simulations. To determine the thermodynamic equilibrium, one would need to fully sample both unbinding and rebinding events, which is clearly beyond the scope of the present work.

Partial dissociation is reversible, with truncated peptides mostly sampling conformations away from the binding groove, but also inside the groove to a lesser extent. Figure 3 shows RMSD timeseries from simulations of Δ1C and Δ2N, the two smallest deletions leading to partial dissociation of the central residues. For the C-terminally truncated peptide (Figure 3A), we observe three dissociations and two rebinding events; the peptide spends around 5% of the trajectory time with both termini bound. In the case of Δ2N (Figure 3B), the peptide spends roughly 35% of the trajectory time in the vicinity of the binding groove, although it never goes back to its initial position. Even though a similar behavior was observed in each of the five independent simulations, significantly longer trajectories would be required for statistically more accurate results. Larger peptide truncations (Δ2C, Δ3C, Δ3N) show a behavior similar to that of Δ2N, with the peptide sampling two sets of conformations, either away from or close to the groove.

Taken together, our results suggest that MHCI loaded with low-affinity peptides can exist in two states: one where both peptide termini are near or inside the binding groove, and another where one terminus is displaced from the groove (see structures in Figure 3). These results are consistent with NMR experiments of MHCI molecules loaded with variant peptides (63), which showed that MHCI at the cell surface exist as a mixture of a major state and a minor state, the former being more hydrated. Since the same study showed that the melting temperature of MHCI increases with the prevalence of the less hydrated, minor state, we speculate that partial peptide dissociation might be the first step in the internalization and recycling of surface-exposed MHCI molecules.

Simulation studies of MHCI often compare dynamics in the PL and PD forms (14, 15, 25, 33). However, in vivo, MHCI molecules are bound to (usually low-affinity) peptides after their synthesis and assembly with β₂m and calnexin (39), before peptide editing takes place. We compared the dynamics of the PD and PL forms to those of MHCI loaded with various truncated peptides to assess if MHCI^PD is an appropriate proxy for low-affinity peptides.

Figure 4 shows the width of the B*44:02 binding groove, measured at three positions, as a function of peptide cargo. Peptide release widens the binding groove by 2.5 Å at d₁, measured near the F-pocket, and by 0.5 Å at d₂, measured at the center of the groove. However, no significant change is observed at d₃, near the A-pocket. In addition to this opening, the amplitude of the fluctuations of the groove dimensions is larger in the absence of peptide. Interestingly, the effects of N-terminal deletions are very different from those of full peptide removal. The groove remains mostly unchanged at d₁ (F-pocket) and d₂ (center), but narrows by 1 Å at d₃ (A-pocket). This partial collapse suggests that MHCI^PD is not an appropriate surrogate for studying the effects of N-terminal peptide editing. The effects of C-terminal deletions are more similar to those of peptide removal, albeit less pronounced. As more residues are removed, the binding groove widens at d₁ (F-pocket), as expected, and fluctuations around the average increase at d₂ (center); d₃ (A-pocket) remains mostly unchanged. In conclusion, these results post a cautionary note when using MHCI^PD to study peptide editing; while it might be warranted for C-terminal editing by tapasin (since the binding groove widens in the PD form as it does for the ΔC peptides), groove dynamics near the N-terminus depend not only on the loading state of the A-pocket, but also on that of the opposite F-pocket.

To understand which domains of B*44:02 are influenced by peptide truncation and to what extent, we analyzed protein C_α RMSF for each system under study. These are shown along B*44:02 sequence in Figure 5 for the α₁α₂ binding groove domain, and in Figure 6 for the α₃ domain. As expected, the binding groove is the region whose fluctuations are most impacted by peptide cargo modification. The α₂₋₁ helix segment known to interact with tapasin (7, 14, 15) and the N-terminal portion of α₁ located in the vicinity of the A-pocket are its most affected regions. Stretch 51–54, in α₁ N-terminus, was proposed as a hinge that pivots in the transition from peptide-receptive to peptide-loaded MHCI (30, 31) in murine allele H2-L^d. In previous work (14, 15), we proposed a peptide-editing mechanism according to which peptides and tapasin compete for the α₂₋₁ region. The results presented here yield a refined view, in the sense that they suggest there could be two sensing mechanisms at work to distinguish properly loaded MHCI at the A- and F-pockets. The fluctuations in the two affected regions are independent: N-terminal peptide truncation recovers the increased RMSF in α₁ (Figure 5B), while C-terminal truncation yields increased fluctuations in α₂₋₁ (Figure 5C).

Possible long-range effects of cargo modification on α₃ are important to consider given that this domain contacts tapasin in the PLC and that allosteric effects involving α₃ have been proposed (24) as a mechanism for peptide sensing by tapasin. However, apart from the Δ1N peptide (see below), we could not observe any significant effect of peptide truncation on α₃ fluctuations (Figure 6). It is possible, however, that signal transduction events are slow compared to the length of our trajectories (1.0 µs) and that increased sampling could yield a different picture. Another possibility is that subtle conformational changes in α₃ are not reflected in the magnitude of its fluctuations.

To complement our analysis of fluctuations, we computed configurational entropies for the α₁α₂ and α₃ domains of B*44:02. Results are shown in Figure 7, where ΔS is the difference between the entropy of each system and the entropy of the reference system MHCI^PL. The important increase of entropy (100–300 J/[K mol]) observed for α₁α₂ (Figure 7A) is consistent with the increased fluctuations of the individual residues, as discussed above. The small statistical uncertainties (less than ±10 J/[K mol] for most systems) indicate convergence, i.e., that 5.0 µs of MD sampling is sufficient to assess B*44:02 α₁α₂ dynamics.

Configurational entropy differences in the α₃ domain (Figure 7B) are also consistent with residue RMSF: ΔS is very small for most systems. A modest but significant increase (63 J/[K mol]) is observed, however, for Δ1N. While this could be an indication of a signal being transmitted between the binding groove and α₃, it also raises the question of why this difference is not observed for any other N- or C-terminal deletion. One possibility is that, since allostery-related protein dynamics often take place on the microsecond–millisecond timescale, our simulations lack the sampling necessary to describe such slow events. Another possibility is that a certain stiffness of the binding groove is required for signal transduction, and that the increased α₁α₂ plasticity caused by deletions larger than a single residue may disrupt this sensing mechanism; however, this would not explain the apparent absence of allosteric signal in the Δ1C system. Finally, a third possibility is that the presence of tapasin, the tapasin–ERp57 conjugate, or other accessory proteins forming the PLC, is necessary to induce conformational changes in MHCI that mediate the allosteric mechanism.

Antigenic peptide specificity is often driven by their N- and C-terminal residues (65), which also determine their stability in the MHCI binding groove. In B*44:02^DPA, residues E2 and F9 extend their side chains toward the floor of the binding groove, occupying the A- and F-pockets, respectively. Since their truncation leads to partial peptide dissociation in our simulations (see Figure 2), we expect their side chains to be key determinants of peptide stability. To test this hypothesis, we have carried out additional simulations of B*44:02-bound modified DPA peptides. These have had one or two residues replaced by glycines at either the N- or the C-terminus. We denote these systems MHCI^δ1N, MHCI^δ2N, and MHCI^δ1C. As for the other modified peptides, trajectories totaling 5 µs of MD simulations (5 × 1 µs) were acquired for each system.

Partial peptide dissociation was observed for two of the three systems under study (Figure 8). Interestingly, side chain removal at the N-terminus leads to a partial dissociation pattern that is significantly different from the one observed for the corresponding truncated peptides. δ1N does not dissociate at all, showing a stable binding almost identical to that of unmodified DPA. δ2N exhibits partial dissociation, but to a much lesser extent than Δ2N: only E1 and E2 vacate the binding groove. Therefore, at the N-terminus, both the peptide main chain and the E1 and E2 side chains contribute to MHCI binding at the A-pocket. Conversely, at the C-terminus, the δ1C dissociation pattern is similar to that of Δ1C, with residues 4–9 dissociating from the groove. The F9 side chain is thus a critical determinant of peptide binding at the F-pocket.

We performed 5 additional, independent, 1-µs MD simulations of MHCI^Δ1C in complex with tapasin, using our previously determined structure of tapasin–MHCI (14, 15) and truncating the full-length DPA peptide present in that structure. This allows us to compare the dynamics of B*44:02 loaded with a low-affinity peptide in the free and tapasin-complexed forms. Results (Figure 9) show that the binding groove width near the F-pocket (d₁, as shown previously in Figure 4), is wider by about 1.0 Å in the presence of tapasin. This is consistent with the “tug-of-war” mechanism we previously described (15), according to which tapasin selects peptides depending on their C-terminus by pulling MHCI β₇ underneath α₂₋₁ to open the groove and thus accelerate peptide release.

Next, to obtain a more detailed and spatially resolved picture, we used force distribution analysis (FDA) to compute the residue-residue pairwise forces between B*44:02 and the Δ1C peptide in tapasin-complexed and -free MHCI (Figure 10). Overall, forces between MHCI and the peptide are on average less attractive (by 150 pN) in the tapasin–MHCI complex, showing that the opening of the groove induced by tapasin decreases the affinity of MHCI for its low-affinity cargo. Interestingly, considering individual pairwise forces at the residue level reveals that tapasin affects the peptide not just in the vicinity of the α₂₋₁ region (C-terminus), but along its whole sequence. Indeed, the largest force difference is observed between the peptide at position R5 and B*44:02 R97 at the bottom of the binding groove; in the tapasin–MHCI complex, this force is more repulsive (by 148 pN) than observed for B*44:02 in the free form. A less attractive binding force is also observed at the peptide N-terminus (E1). Thus, tapasin-induced widening of the binding groove near the F-pocket (see Figure 9) has long-ranging effects that destabilize peptide binding, increasing the latter’s off-rate. Full peptide release, however, was not observed in the present simulations. Much longer timescales would likely be required to investigate the tapasin-catalyzed unbinding event itself.

Binding groove widening has been suggested before (9) in a biochemical study of tapasin–MHCI. In their mechanism of complex assembly/disassembly, the authors proposed two possible conformational transitions taking place in the MHCI binding groove. The first is an equilibrium between “closed” and “open” states, which exist both in presence and absence of tapasin and could correspond to a tightly packed and loose binding groove, respectively (9, 24, 66). The second is a switch between a canonical and a “tapasin-disrupted” binding groove conformation, which could correspond to the widening we observe in our simulations. Partial peptide dissociation of the Δ1C peptide, in turn, would prevent the transition from the “open” to the “closed” state. Interestingly, the same study (9) also found that mutating certain high-affinity peptides at their N-terminus caused them to exhibit higher tapasin sensitivity. This indicates that tapasin impacts MHCI across the full length of the peptide, which is consistent with our above FDA results that reveal destabilizing effects at E1 and R5.