Understanding Allosteric Stabilization of pMHC by Peptide Ligands Through Molecular Dynamics Simulations

Introduction

The adaptive immune response is a critical component of the immune system, providing highly specific and long-lasting protection against pathogens. Cytotoxic T lymphocytes play a crucial role in this response by eliminating infected or cancerous cells (1,2). Cytotoxic T lymphocytes recognize and kill target cells through the interaction of their T cell receptors (TCRs) with peptide-loaded major histocompatibility complex (pMHC) molecules presented on the surface of antigen-presenting cells (APCs) (3). This interaction is highly specific, with the TCR recognizing both the peptide antigen and the MHC molecule, a phenomenon known as MHC restriction (4,5).

There are two classes of MHC molecules: MHC class I (MHC-I), which presents endogenously derived peptides to CD8+ T cells, and MHC class II (MHC-II), which presents exogenously derived peptides to CD4+ T cells (6). MHC-I molecules consist of a heavy chain, a light chain (β2-microglobulin, β2m), and a bound peptide. The heavy chain comprises three domains (α1, α2, and α3), with the α1 and α2 domains forming the peptide-binding groove. The peptide is anchored within this groove, while β2m interacts with the binding groove and α3 domain of the heavy chain, stabilizing the overall structure (7,8). 

Protein dynamics play a critical role in the interaction between TCRs and pMHC complexes (9). Both MHC and TCR molecules exhibit inherent flexibility, which is essential for their function (10,11). The dynamics of MHC-I molecules thus influence peptide binding, presentation, and T cell recognition. Once bound, the peptide ligand itself can also modulate the dynamics of pMHC, affecting its flexibility (12), with significant implications for TCR recognition and subsequent T cell activation (13).

The stability of pMHC complexes is also essential for efficient TCR recognition. Peptide ligands were shown to stabilize pMHC by interacting with specific pockets (A, B, C, D, E, and F) within the peptide-binding groove (14). Different peptides modulate the dynamics and flexibility of pMHC complexes in distinct ways, thereby altering the protein’s energy landscape and influencing overall stability. Notably, pMHC stability has been positively correlated with its flexibility (15). 

Furthermore, β2m binding and dissociation have been used as accurate measures of pMHC stability, highlighting the importance of the human leukocyte antigen (HLA)-β2m interface (16,17). These findings suggest that peptides can allosterically affect the HLA-β2m interface, thereby influencing the stability of the pMHC complex. However, the exact mechanisms underlying this allosteric effect remain largely unknown.

Molecular dynamics (MD) simulations provide detailed insights into the conformational changes and interactions that govern pMHC stability and flexibility (18-22). Given the potential of MD simulations to unravel the allosteric mechanisms by which peptides stabilize pMHC complexes, this study aimed to investigate these mechanisms using MD simulations and shed light on the complex interplay between peptide binding, pMHC dynamics, and TCR recognition.

Materials and Methods

Molecular Dynamics Simulations of Peptide-HLA*02:01 Complexes

Explicit-solvent MD simulations were performed using GROMACS version 2023.4 (GROMACS Development Team, Stockholm, Sweden) (23). Hydrogen mass repartitioning was employed to enable a 4 fs integration time step (24). Each system was simulated three times, with each simulation covering a total time of 750 ns. The AMBER99SB-ILDN force field (25) was used in combination with the TIP3P water model (Jorgensen Laboratory, Yale University, New Haven, CT, USA). 

Systems were neutralized with sodium/chloride ions to a final ionic strength of 0.15 M. Simulations were run in the isothermal-isobaric (NPT) ensemble with the pressure maintained at 1 bar and the temperature at 310 K. A non-bonded interaction cutoff of 12 Å was applied for van der Waals interactions. A conformation was saved every 5000 steps. 

Initial protein structures were obtained from Protein Data Bank (PDB; Research Collaboratory for Structural Bioinformatics, Piscataway, NJ, USA) (26) with accession codes of 3UTQ, 5N1Y, and 5C0F (24). The peptide-free system was obtained by simply removing the peptide ligand from the MHC structure deposited in 3UTQ.

Root Mean Square Deviation and Root Mean Square Fluctuation Analysis 

Root mean square deviation (RMSD) was calculated for each simulation to assess the deviation of the molecular structure from its initial conformation over time. The calculation was based on alpha carbon coordinates, with the first frame of the trajectory as the reference. The equation used for RMSD calculation was as follows:

where N is the number of atoms, xi is the position of atom i in the current frame, and xiref is the position of atom i in the reference frame.

Root mean square fluctuation (RMSF) was computed to quantify the atomic flexibility around their average positions during the simulation. Root mean square fluctuationfor each atom was calculated using the following equation:

where T is the total number of frames, xit is the position of atom i at time t, and xi is the average position of atom i over the entire simulation. Root mean square fluctuation was computed using only the final 500 ns of each trajectory.

Essential Dynamics Analysis

Essential dynamics analysis (EDA), also known as principal component analysis (PCA), was performed to identify cooperative global motions of amino acids. First, a covariance matrix (C) of alpha carbon coordinates was computed:

where X is the coordinates of all alpha carbon atoms, with a dimension of 3N by 3N. <X> denotes the average position matrix over the equilibrium portion included in the RMSF analysis as well. This matrix is then diagonalized to identify its eigenvalues and eigenvectors:

Here, λ is a diagonal matrix consisting of eigenvalues, and T denotes eigenvectors. The eigenvectors are principal components (or individual modes of motion) that describe the essential dynamics of the protein complex.

From each mode of motion, individual fluctuations of alpha carbons can thus be obtained. Using these fluctuations, elements of dynamic cross-correlation matrices (DCCMs) can be constructed using the following equation:

Here, Cij is the cross-correlation value between residues i and j, ∆ri and ∆rj represent the positional fluctuations of alpha carbons in the respective modes of motion. The values of Cij range from -1 to 1, denoting the maximum levels of anti-correlation and correlation, respectively. To simplify the analysis and focus on the top global modes of motion that describe the majority of variance in motion, C matrices were computed for only the top 20 global modes of motion. A weighted average of these matrices was then computed using eigenvalues as weights to yield a single DCCM for each simulation trajectory. 

Root mean square inner products (RMSIPs), also termed subspace overlaps, of two simulation trajectories were also computed to quantify the levels of (dis)similarities between them, using the following equation:

Here, ηi and νj are eigenvectors associated with modes i and j, respectively. 

Calculation of Pairwise Residue Interaction Energies

Non-bonded interaction energies between contacting amino acid pairs throughout the simulation trajectories were computed using a modified version of gRINN (get Residue Interaction eNergies and Networks) (28). First, a cutoff distance of 10 Å was used to identify pairs of amino acids to be included in the calculation, based on positions of their alpha carbons. Then, the “rerun” feature of gmx mdrun was used to compute non-bonded van der Waals and electrostatic interaction energies between these identified amino acid pairs. Energies were initially reported in kJ/mol by GROMACS and subsequently converted to kcal/mol before further analysis. 

Construction of Protein Energy Networks

Once interaction energies were computed, a protein energy network (PEN) was constructed by taking each amino acid as a node, each interaction as an edge, and van der Waals interaction energies as an edge weight. Edges were added between two residues if they were covalently-bound or interacted with each other in the respective conformation with a van der Waals interaction energy of at least 1 kcal/mol in magnitude (either ≤ –1 kcal/mol or ≥ 1 kcal/mol). Edge weights were given a value of 1 if the two interacting pairs are covalently bound, and a value of Xij if otherwise. Xij was computed as follows:

Here, Xij is the edge weight, εij is the non-bonded van der Waals interaction energy, and εatt is the array of negative interaction energies. Using this formula enables the assignment of higher weights to more negative (or more attractive) interactions. A PEN was constructed for each conformation in the simulation trajectory. Once a PEN was constructed, betweenness-centralities (BCs) of each node were computed using Brandes’ algorithm (29), based on shortest paths identified using Dijkstra’s algorithm (30).

Statistical Analysis

Statistical significance was defined as p<0.05. One-way analysis of variance (ANOVA) and Student’s t-test were used as appropriate.

Tools and Software Packages

All analyses were performed within a Python-based data analysis environment using the following packages: NumPy (31), Matplotlib (32), Seaborn, Pandas, and ProDy (33).

Results

Weaker Interactions with the F Pocket Lead to Lower pMHC Conformational Stability

Three peptides (ALW, MVW, and RQW) bound to HLA-A*02:01 were studied using molecular dynamics simulations. Each system underwent three 750-ns simulations, including a peptide-free control as reference. The selection of these particular peptides was based on comprehensive experimental data reported by Hopkins et al. (31). This study is notable as it uniquely provides not only half-life measurements, indicative of thermal stability, but also detailed insights into the conformational flexibility of these pMHCs through high-pressure perturbation experiments. Such extensive experimental data, covering both stability and dynamic conformational responses, is uncommon in pMHC literature, where studies often report only affinity or thermal/cell surface stability, leaving the role of conformational flexibility less defined. 

Given that all three pMHC systems (ALW, MVW, and RQW with HLA-A*02:01) were found to behave differently in high-pressure or temperature perturbation experiments and previous molecular dynamics simulations (31), they serve as ideal benchmarks for a more granular investigation of the effect of peptide binding on the overall pMHC structure, as intended here. Notably, the experimental data from Hopkins et al. (31) highlighted distinct biophysical characteristics for these peptides; for instance, the ALW-HLA-A*02:01 complex showed physicochemical indicators (a significantly higher ΔCp) suggestive of greater conformational flexibility compared to the others. Furthermore, these peptides are valuable benchmarks due to observed discrepancies between their experimental stability and affinity data and those predicted by standard computational tools, underscoring the need for detailed mechanistic studies. A detailed description of the characteristics of these peptides can be found in Supplementary Text.

Figure 1A shows stable fluctuations after equilibration, with the ALW peptide and peptide-free system displaying higher fluctuation levels. This suggests ALW's effect on pMHC dynamics is less significant than that of other peptides. The effect on stability of the whole complex is often directly linked to the binding stability of the peptide throughout the simulations. To probe the stability of each peptide residue, we next computed the RMSF values of each peptide residue (position) in each simulation. 

Figure 1:

Conformational stability and dynamic profiles of simulated peptide-loaded Major Histocompatibility Complex (pMHC) systems

A) Root mean square deviation (RMSD) profiles over time for peptide-free and peptide-bound HLA-A02:01 systems. Darker lines represent the mean RMSD from three replicate simulations; shaded areas indicate standard deviation ranges.
B) Root mean square fluctuation (RMSF) profiles of individual peptides in each simulation. Each line corresponds to one replicate simulation.
C) Root mean square fluctuation profiles of HLA-A02:01 and β2-microglobulin (β2m) residues. Darker lines indicate the average per-residue RMSF across three replicates; shaded areas denote the standard deviation.
D) Subspace overlap (Root Mean Square Inner Product, RMSIP) values between individual replicate simulations. An RMSIP of 1 denotes identical essential dynamics, while 0 indicates no similarity. Axis labels show PDB four-letter codes for each peptide–HLA-A*02:01 complex, followed by the replicate number. 3UTQ corresponds to the ALW peptide, 5N1Y to the MVW peptide, and 5C0F to the RQW peptide.

Figure 1:

Conformational stability and dynamic profiles of simulated peptide-loaded Major Histocompatibility Complex (pMHC) systems

This analysis revealed partial peptide detachment in one simulation of each system: at the C-terminus (F pocket) for ALW and MVW, and at the N-terminus (B pocket) for RQW. In peptide-free and ALW simulations, HLA and β2m residues showed higher variation in positions 70-85 and 140-150, corresponding to α-1 and α-2 regions of the F pocket. Results indicate that strong peptide C-terminus-F-pocket interactions are crucial for conformational stability.

Peptide Ligand Leads to Reproducible Global Dynamics Profiles 

Our findings showed that peptide sequence and interactions affect fluctuations of peptide and HLA residues, and influence regions beyond the peptide binding groove. This is not clearly visible in the RMSF plot (Figure 1C), except in parts of the α3 domain, possibly due to thermal fluctuations masking potential dynamic differences. We performed EDA of simulation trajectories to better understand similarities/differences between global dynamics profiles of different systems.

Essential dynamics analysis provides motional modes and their contributions to overall variance. We analyzed the top 20 motions describing global coordinated movements and calculated similarities between these modes using the RMSIP metric. Figure 1D shows that the replica simulations of each system exhibited the highest similarity, indicating strong reproducibility. This enables a reliable comparison between different systems. Subspace overlaps between different systems exhibit moderate similarity levels (0.5–0.6), indicating that while peptides do not create unique cooperative domain motions, they produce distinct motional profiles. RQW and MVW peptide simulations showed relatively higher similarities compared to others.

Stronger Peptide Binding in F Pocket Leads to More Extensive Changes in Dynamic Residue Couplings

Having confirmed that replica simulations explored similar essential subspaces, we next investigated residue couplings, comparing peptide-loaded and peptide-free systems. For this, we generated average DCCMs for each system using the top 20 global motional modes. Subtracting the DCCM of the peptide-free system from each peptide-loaded system yielded difference DCCMs (Δ­DCMMs), visualizing changes in residue couplings upon peptide binding (Figure 2). 

Figure 2:

Residue coupling patterns and peptide-induced differences in residue coupling in HLA-A*02:01 systems.

Dynamic cross-correlation matrices (DCCMs) and difference DCCMs (ΔDCCMs) are shown for peptide-loaded and peptide-free systems.
Top left: DCCM (upper triangle) and ΔDCCM (lower triangle) for the ALW peptide system.
Top right: DCCM (upper triangle) and ΔDCCM (lower triangle) for the MVW peptide system.
Bottom left: DCCM (upper triangle) and ΔDCCM (lower triangle) for the RQW peptide system.
Bottom right: Schematic legend defining intra- and inter-domain interaction regions on the DCCM heatmap. Domain labels indicate the peptide-binding groove (BG; α1 and α2 domains), α3 domain, α2-microglobulin (α2m), and inter-domain regions.

Figure 2:

Residue coupling patterns and peptide-induced differences in residue coupling in HLA-A*02:01 systems.

Consistent with similar RMSF profiles and RMSIP values above 0.5 (Figure 1D), the systems exhibited similar overall residue coupling profiles. However, ΔDCMM plots also revealed subtle peptide-specific differences. ALW binding induced relatively smaller coupling changes, increasing both correlations and anti-correlations between the HLA binding groove and α3/β2m. Notably, ALW binding showed minimal changes in intra-α3 couplings, with increased intra-β2m correlations but decreased β2m-α3 correlations. 

Conversely, MVW and RQW exerted similar, but stronger, effects than ALW. They increased positive correlations between the HLA binding groove and α3/β2m, shifted intra-α3 couplings toward anti-correlation, and increased β2m-α3 coupling. In summary, stronger F pocket peptide binding led to more extensive changes in MHC dynamic coupling, particularly an increase in the positive correlation between the peptide-binding groove and distal ­domains.

Peptide Binding Induces Perturbations in 

Intra- and Inter-Chain Residue Interaction Energies

We identified "consistent pairs" that exhibited significantly similar interaction energies across replicas using one-way ANOVA. In other words, pairs whose interaction energies did not differ significantly (p>0.05) across replica simulations (defined as those with a p-value > 0.05 from the ANOVA test) were identified (the full list of pairs and their corresponding p values can be found in Supplementary Data file under Tab “Consistent Pairs”). Finally, we identified pairs that showed statistically significant (p<0.05, Student’s t-test) energy shifts upon peptide binding, defined as those with a p-value < 0.05 from a Student’s t-test, and which also had  with an absolute interaction energy difference of at least 2 kcal/mol upon peptide binding toward attractive or repulsive ranges (Figure 3) (the full list of pairs and their corresponding p values can be found in Supplementary Data file under Tab “Affected upon Peptide Binding”). 

Figure 3:

Residue pairs exhibiting significant shifts in interaction energies upon peptide binding.

Residue pairs that showed statistically significant shifts in nonbonded interaction energies upon binding of ALW, MVW, or RQW peptides are shown.
(A, B) ALW-bound system; (C, D) MVW-bound system; (E, F) RQW-bound system.
Red spheres represent residue pairs with interactions shifted toward stronger attraction; purple spheres represent residue pairs with interactions shifted toward increased repulsion.

Figure 3:

Residue pairs exhibiting significant shifts in interaction energies upon peptide binding.

Strikingly, ALW exhibited a smaller effect compared to the other peptides. ALW induced shifts toward attraction only within the binding groove, not extending to distal domains or the β2m-α3 interface. It also caused shifts toward repulsion in this interface, decreasing attraction between Val12(BG)-Ser34(β2m) and Ile23(BG)-Phe63(β2m). Conversely, MVW and RQW showed much more extensive effects, both within the binding groove and on interacting pairs between β2m and α3, as well as at the bottom of the binding groove. Their effects were also highly similar, affecting the same pairs. Interactions becoming more attractive included Phe9(BG)-Phe57(β2m), Ser11(BG)-Glu55(β2m), and Tyr11(β2m)-Pro235(α3), while those becoming more repulsive included Val12(BG)-Ser34(β2m), Thr10(BG)-Gly56(β2m), and Val25(BG)-Gly56(β2m).

Stabilization of the F Pocket by Peptides Leads to Stabilized Interactions Near the Core of the MHC 

Motivated by the differential effects of peptides on the pMHC protein-protein interface, we investigated how these changes translate to individual residue behavior. Since edge weights reflect interaction attractiveness or repulsiveness, PENs are useful for probing stability. We computed betweenness centrality (BC) values for each residue, representing the average number of shortest paths through that residue. Here, residues with optimized interactions (and thus predominantly negative/attractive non-bonded interaction energies with their surroundings) have higher BC values.

These scores were subsequently compared to overall residue conservation scores obtained from the CONSURF database (35). A scatter plot of conservation scores versus mean BC levels is shown in Figure 4. A key takeaway from this plot is the high conservation levels of HLA residues with the highest BC levels. These include Pro235(α3), Tyr209(α3), Pro210(α3), Trp244(α3), Ala245(α3), and Val28(BG).

Figure 4:

Relationship between residue betweenness centrality and evolutionary conservation.

Scatter plot of mean median betweenness centrality (BC) values for residues from peptide-bound simulations plotted against their conservation scores obtained from the CONSURF database (PDB ID: 5N1Y). Lower conservation scores indicate higher evolutionary conservation. Chain identifiers following underscores denote the respective chains in the simulated complex; chain A corresponds to HLA-A*02:01.

Figure 4:

Relationship between residue betweenness centrality and evolutionary conservation.

Finally, we identified residues with statistically significant BC differences between systems, which were defined as those having a p-value < 0.05 from a using one-way ANOVA test. From this group of significant residues, we (p<0.05) and selected the top 10 residues with the largest BC increases or decreases upon peptide binding (Figure 5), based on median BC levels (the full list of residues with their corresponding p values can be found in Supplementary Data file under Tab “Betweenness Centralities”).

Figure 5:

Residues with the highest peptide-induced changes in betweenness centrality.

Box plots show residues (Residue_id) with the greatest increases (top) and decreases (bottom) in median betweenness centralities (BC) upon peptide binding. Data are presented for ALW, MVW, and RQW peptide-bound systems and the peptide-free control. Chain identifiers denote the respective chains: chain A for HLA-A*02:01 and chain B for β2-microglobulin (β2m). The p values denoting statistical significance of the BC differences per residue are denoted above each box plot

Figure 5:

Residues with the highest peptide-induced changes in betweenness centrality.

Two key findings emerged from this analysis. First, residues 64 B (β2m threonine) and 235 A (α3 proline), both located in distal domains (Figure 6), showed the largest BC increases upon peptide binding. All peptides increased 64 B's betweenness, suggesting an allosteric effect near the β2m/pMHC core, reinforcing surrounding interactions. However, 235 A's BC increase was observed only with RQW and MVW, not ALW. 

Figure 6:

Spatial distribution of residues most affected by peptide binding based on betweenness centrality.

Residues corresponding to the top (red) and bottom (blue) ranked residues from Figure 5 are mapped onto the HLA-A*02:01 structure. Red spheres denote residues with the largest BC increases; blue spheres indicate residues with the largest BC decreases upon peptide binding.

Figure 6:

Spatial distribution of residues most affected by peptide binding based on betweenness centrality.

The first finding leads to the second key finding: for most residues in Figure 4, BC patterns were shared between the ALW and peptide-free systems, as well as between the RQW and MVW peptides, consistent with the lack of structural effects upon ALW binding. Consequently, the behavior of Pro235(α3) is particularly noteworthy; its high mean BC as shown in Figure 4 underscores its general importance to the structural integrity of the complex, while its significant and selective increase in BC upon binding of the strongly stabilizing peptides RQW and MVW further pinpoints it as a key residue involved in mediating peptide-specific allosteric effects on pMHC stability.

Discussion

This study investigated the influence of different peptide ligands on the structural dynamics and conformational stability of HLA-A*02:01 using molecular dynamics simulations. Essential dynamics and protein energy network analyses provided insights into peptide-induced changes in MHC-I and the modulation of complex dynamics.

The stability of pMHC has long been recognized as a critical factor for immunogenicity (39,40). Consequently, characterization of pMHC stability is of great interest for immunotherapeutic applications (41). It is understood that peptide binding is not solely dictated by interactions within the peptide-binding groove. Instead, the entire MHC complex, including more distal domains, contributes to the binding event and the overall stability of the pMHC assembly (16,42). Indeed, previous simulation studies have shown the importance of regions beyond the canonical binding site in modulating peptide-MHC interactions (19,42).

Despite this broader understanding, a significant portion of computational studies aimed at predicting or characterizing pMHC stability has predominantly focused on the interactions within the peptide-binding groove. This focus is often a necessary simplification to enable predictions at larger scales (44-47). However, such a narrow perspective may overlook crucial contributions from other regions of the complex to its overall stability. To our knowledge, the present study is the first to employ molecular dynamics simulations to systematically evaluate pMHC stability with a specific focus on protein-protein interaction interfaces throughout the entire complex (e.g., the heavy chain-β2m interface and intramolecular domain contacts). This approach enabled the identification and characterization of allosteric mechanisms through which peptides differentially modulate stabilization across these distal interfaces, extending our understanding beyond direct peptide-groove interactions.

Our results are consistent with prior research that highlights the tight connection between peptide ligands and MHC-I. Peptides typically anchor deeply into the MHC-I binding groove, resulting in MHC-I stability and flexibility depending on the peptide cargo (49). Peptide-free MHC is known to be unstable (49), and peptide binding stabilizes and dampens complex dynamics, affecting distal β2m and α3 domains. MHC conformational stability depends on peptide C-terminus interactions with the F-pocket (50). Our findings offer new insights into this stabilization. We have shown that stable F-pocket binders (MVW and RQW) induce larger increases in positive correlations between the binding groove and α3, as well as in both positive and negative correlations between the binding groove and β2m. This strengthens crosstalk between distal domains and the binding groove, potentially through increased attraction between HLA and β2m, as revealed by our non-bonded interaction energy analysis. Weaker ALW binding even increased repulsiveness at the β2m-HLA interface.

The highly polymorphic nature of MHC-I molecules has significant implications for their structural dynamics and stability (20,51-53). Certain MHC alleles are more dependent on the chaperone protein Tapasin for proper peptide loading compared to others (54-57). Micropolymorphisms within the F-pocket region of the peptide-binding groove may lead to differential dependencies on Tapasin for binding of high-affinity peptides (57-60). These subtle sequence variations in the F-pocket can lead to distinct global conformational dynamics of the MHC complex, even in the peptide-free state. Specifically, some F-pocket polymorphisms may favor a more stable and rigid conformation of the MHC, while others may result in a more flexible and dynamic structure (61). This allele-specific behavior of the peptide-free MHC likely has consequences for the conformational stabilization induced by the binding of different peptide ligands, as observed in our findings. Future studies that explore the interplay among F-pocket polymorphisms, Tapasin dependency, global dynamics, and complex stability will provide valuable context for interpreting the conformational stabilization mechanisms revealed by our molecular dynamics’ simulations.

Furthermore, major stability differences among HLA genes themselves may also imply differential behaviour to peptide binding. For instance, HLA-C alleles are known to exhibit lower cell surface expression and stability levels compared to HLA-A and HLA-B alleles (62,63). This reduced stability of HLA-C has also been associated with less efficient responses against viral infections (64). Our previous local frustration analysis of human MHC-I allele structures highlighted a less stable F pocket in HLA-C (65). Weaker interactions between the F pocket and peptide cargo may thus represent a general phenomenon in HLA-C, resulting in a failure to strengthen interactions with distal domains of the complex, as observed in the case of the ALW peptide in this study.

Finally, it is important to note that we did not observe a strong correlation between experimental or theoretical measures of thermal stability (i.e., melting temperature) or kinetic stability (i.e., half-time) and our findings related to conformational stability (see Supplementary Text). This may be explained by discrepancies between experimental measurements on the same peptide-HLA pair, or, more specifically, a lack of sufficient data points, which restricts the reliability of such comparisons. Alternatively, the thermal and conformational stability of MHC molecules may not be directly related, as is also the case for many other proteins (66,67). Previous work in this field has explored both kinetic (i.e., over time) and thermal stability of pMHC, and a positive correlation between these two types of stability metrics has been identified (68,69). It is also worth emphasizing that classical MD simulations in the isothermal-isobaric ensemble (i.e., at constant temperature and pressure) are inherently limited to mimic experimental conditions used in both types of stability measurements since simulations cover neither the full-time scale (often hours) nor the temperature ranges (increased over time) used in these experiments. To this end, predicting and evaluating MHC thermal stability alongside molecular simulation findings remains a challenge. Application of a similar computational strategy to a higher number of pMHC complexes with reliable experimental structures and both thermal and kinetic stability data may lead to more comprehensive insights into peptide-induced stabilization in this regard.

Conclusion

In summary, this study highlights the intricate relationship between peptide ligands and the conformational dynamics of HLA-A*02:01, emphasizing the crucial role of stable interactions within the F-pocket for achieving structural integrity and stability. The observed enhancements in dynamic coupling and crosstalk between the peptide-binding groove and distal domains underscore the complexity of MHC-I-peptide interactions and their influence on overall protein behavior. These enhancements may result in a shift of attractive interactions towards the HLA-β2m interface, ultimately contributing to increased conformational stability of the whole complex. To this end, our findings also challenge the assumption that thermal stability directly correlates with conformational stability, highlighting the need for further research to elucidate these relationships. The insights gained from our molecular dynamics simulations not only reinforce existing literature but also pave the way for future studies aimed at understanding the nuances of MHC-peptide interactions, which could ultimately inform therapeutic strategies in immunology and vaccine development.