Markov state modeling reveals alternative unbinding pathways for peptide-MHC complexes

Abstract

Peptide binding to major histocompatibility complexes (MHCs) is a central component of the immune system, and understanding the mechanism behind stable peptide-MHC binding will aid the development of immunotherapies. While MHC binding is mostly influenced by the identity of the so-called anchor positions of the peptide, secondary interactions from nonanchor positions are known to play a role in complex stability. However, current MHC-binding prediction methods lack an analysis of the major conformational states and might underestimate the impact of secondary interactions. In this work, we present an atomically detailed analysis of peptide-MHC binding that can reveal the contributions of any interaction toward stability. We propose a simulation framework that uses both umbrella sampling and adaptive sampling to generate a Markov state model (MSM) for a coronavirus-derived peptide (QFKDNVILL), bound to one of the most prevalent MHC receptors in humans (HLA-A24:02). While our model reaffirms the importance of the anchor positions of the peptide in establishing stable interactions, our model also reveals the underestimated importance of position 4 (p4), a nonanchor position. We confirmed our results by simulating the impact of specific peptide mutations and validated these predictions through competitive binding assays. By comparing the MSM of the wild-type system with those of the D4A and D4P mutations, our modeling reveals stark differences in unbinding pathways. The analysis presented here can be applied to any peptide-MHC complex of interest with a structural model as input, representing an important step toward comprehensive modeling of the MHC class I pathway.

Keywords: Markov state modeling; adaptive sampling; competitive binding assay; peptide–MHC binding stability.

Fig. 1.

Overview of the simulation framework. (A) The exploration stage involves running umbrella-sampling simulations along the z-dist reaction coordinate, which approximates the unbinding direction. $B_i$ is the energy bias, while $k$ is the force constant. The $\beta$-sheet floor of the MHC (light blue) is aligned to the XY plane, and then the Z coordinate is used to define z-dist. The truncated portion of the MHC (light gray) is not included in any of the simulations. (B) The connection stage involves running unbiased simulations in an adaptive sampling fashion until most of the states are connected. Restarting conformations are chosen by analyzing the trajectories in a dimensionality-reduced space using TICA that adequately captures the binding/unbinding pathway. Then the selection of conformations is biased toward the less densely sampled regions of the TICA space.

Fig. 2.

$\mathrm{\Delta }\mathrm{\Delta }G$ predictions from the mutational analysis. The black dashed line represents the predicted $\mathrm{\Delta }{G}_{WT}$ of $-7.19$ kJ/mol. The gray dashed line represents the separation between predicted binders and nonbinders. Alanine mutations in positions 2, 4, and 9 are all predicted to significantly impair binding, while alanine mutations in positions 1, 5, and 7 are predicted to reduce the binding affinity.

Fig. 3.

Competitive binding assays to determine the ranking of WT, D4A, and D4P. Based on the relative position of the WT curve (green plus) versus the positive control (blue circle), we see that QFKDNVILL is indeed a weak binder to HLA-A*24:02 (${\text{IC}50}_{\text{WT}}\hspace{0.17em}=\hspace{0.17em}$1,600 nM). Upon mutation of D4 to an alanine, inhibition is significantly reduced (${\text{IC}50}_{\text{D}4\text{A}}\hspace{0.17em}>\hspace{0.17em}$6,000 nM) as the D4A curve (red cross) is most similar to the negative control (purple triangle). Upon mutation of D4 to a proline, inhibition is increased (${\text{IC}50}_{\text{D}4\text{P}}\hspace{0.17em}=\hspace{0.17em}$600 nM) as the D4P curve (orange square) is most similar to the positive control.

Fig. 4.

Flux network of unbinding trajectories for the WT system. States 0, 1, 2, and 3 denote the set of associated states that have the peptide in contact to the MHC. State 4 represents the dissociated or unbound state. Size of the nodes (depicted in red) indicates the equilibrium probabilities of each state (${\pi }_i$). (A) The WT system prefers to unbind through detaching first on the C-term end (state 0 to state 1 transition) due to the stronger interactions on the N-term end, which include the aspartic acid in position 4. (B) With a single mutation, the D4A system prefers to unbind through detaching first on the N-term end (state 0 to state 2 transition), and the accessibility of both detachment pathways favors the instability of the D4A system. Note that the MSM model includes all transitions between nodes, in all directions. However, this flux network depicts only trajectories starting from state 0 and reaching state 4 (i.e., unbinding pathways).

Fig. 5.

Representative conformations in the WT system from state 0 (native state) and state 1 (N-term bound state). A and B depict the side views of states 0 and 1, respectively. These states can be distinguished by the location of the C-term of the peptide relative to the MHC binding cleft (i.e., proximity to the F pocket). C and D depict the top views of states 0 and 1, respectively. Peptide’s p4 residue (aspartic acid, D) is depicted in magenta (carbon atoms in magenta; oxygen atoms depicted in red). Other peptide positions are depicted in green. Key MHC residues predicted to interact with p4 are depicted in yellow (carbon atoms in yellow; oxygen atoms depicted in red; nitrogen atoms in blue; hydrogen atoms in white), including lysine 66 (K66), threonine 163 (T163), tyrosine 159 (Y159), and glutamine 155 (Q155). Hydrogen bonds involving any of these residues are depicted in yellow dashed lines.

Fig. 6.

(A) Flux network of unbinding trajectories for the D4P system. The introduction of a proline forces the unbinding starting from the N-term side (state 2). (B) (Blue contour) Phi/psi angles (in radians) of position 4 from WT/D4A unbinding trajectories where the C-term side unbinds first. The bottom region covers states 0 and 1, while the top region covers state 3. (Orange border) Ramachandran plot of accessible phi/psi angles of proline. Unbinding trajectories during the transition from state 1 to state 3 lie in regions that do not overlap with the accessible phi/psi angle of proline. Thus, the unbinding trajectories adopt backbone conformations of p4 that are incompatible with the rigidity of proline. Note that the MSM of D4P (A) includes transitions from state 0 to state 1 and from state 1 back to state 0. However, these transitions are not depicted in the flux network, since none of the paths passing by state 1 were able to progress to state 4.