HLA-DM Stabilizes the Empty MHCII Binding Groove: A Model Using Customized Natural Move Monte Carlo

MHC class II molecules bind peptides derived from extracellular proteins that have been ingested by antigen-presenting cells and display them to the immune system. Peptide loading occurs within the antigen-presenting cell and is facilitated by HLA-DM. HLA-DM stabilizes the open conformation of the MHCII binding groove when no peptide is bound. While a structure of the MHCII/HLA-DM complex exists, the mechanism of stabilization is still largely unknown. Here, we applied customized Natural Move Monte Carlo to investigate this interaction. We found a possible long-range mechanism that implicates the configuration of the membrane-proximal globular domains in stabilizing the open state of the empty MHCII binding groove.


■ INTRODUCTION
The MHC class II complex (MHCII) presents potentially harmful peptides to the immune system. 1 The structure of the peptide-loaded MHCII binding groove is well documented. 2 However, due to its dynamic and unstable nature no structure has been solved for the peptide-free MHCII complex to date; 3 when peptide is absent, MHCII complexes rapidly take on an inactive peptide-averse state. 4 The peptide-exchange factor human leukocyte antigen DM (HLA-DM) is known to stabilize the empty MHCII complex 5 and promote a peptidereceptive state. 6 A structure of HLA-DM in complex with MHCII has been published; 7 however, the mechanism by which HLA-DM stabilizes the peptide-receptive state of the MHCII binding groove is still a subject of debate. 8 Given that one of the key binding sites for HLA-DM on MHCII includes the β2 globular subunit 7,9 and that antibodies detect peptide-induced structural changes in the β2 globular subunit, 10 it was suggested that the HLA-DM chaperoning mechanism involves the long-range transmission of structural changes from the globular subunits to the binding groove ( Figure 1). 11 The hypothesis states that HLA-DM binding causes structural changes in the membrane-proximal β2 subunit of MHCII that can be propagated to the binding groove to modulate peptide binding. 10 We evaluated this experimentally suggested hypothesis and investigated the underpinning mechanistic details of the HLA-DM/MHCII (without peptide) interaction with customized Natural Move Monte Carlo (cNMMC). 12 cNMMC is a protocol for hypothesis based modeling using Natural Move Monte Carlo (NMMC), 13 an established molecular simulation method that has been applied to investigate several biophysical research questions including receptor signaling, 14 prediction of nucleosome occupancy, 15 energy contributions of epigenetic DNA marks, 16 RNA mutations, 17 and refinement of projection images. 18 Importantly, NMMC has previously been validated on the MHC class I complex. Coarse-grained NMMC simulations, using a knowledge potential, showed high agreement with experimental results (AROC = 0.85) in the prediction of peptide binders and nonbinders. 19 Our results suggest that conformational shifts in the MHCII globular domains enable a larger range of motion in the two αhelices lining the binding groove, thereby increasing the likelihood of binding-groove collapse and occlusion of peptide binding. In follow-up simulations we show that the peptidechaperoning protein HLA-DM stabilizes the interaction between the two MHCII globular domains, thereby promoting the open state of the empty binding groove via an indirect long-range mechanism.

■ RESULTS
The Effect of HLA-DM on the MHCII Binding Groove. HLA-DM is known to catalyze the peptide exchange in MHCII molecules. It binds laterally to the α1 region in the binding site and the α2 and β2 globular domains (+DM in Figure 2A) and stabilizes the peptide-receptive form of the MHCII. The structural mechanism by which this stabilization occurs is largely unknown.
Based on the experimentally suggested hypothesis ( Figure  1), we designed a set of customized Natural Moves to investigate a potential involvement of the MHCII globular domains in the stabilization of the binding groove. In a previous study we identified a set of NMMC parameters that led to the collapse of the binding groove largely driven by the plasticity of the β1-1 kink in the β1 helix. 12 Here, we used the same structure (PDB accession code: 4GBX -MHCII in complex with HLA-DM) and NMMC parameters as a starting point for investigating the mechanism underlying MHCII binding-groove stabilization. First, we simulated the collapse of the binding groove by performing simulations in the absence of HLA-DM. In a second set of simulations, we used the same parameters and added the peptide-loading chaperone HLA-DM to test its effect on the stability of the binding groove. While the simulation without HLA-DM (-DM) showed a welldefined bimodal distribution of binding-groove widths indicating the potential for a collapse (leading to a peptideaverse state), the simulation with HLA-DM (+DM) resulted in a single mode around the open state of the binding groove ( Figure 2C). Furthermore, the movement of the MHCII globular domains was restricted by HLA-DM ( Figure S1).

Article
To test whether the configuration of the globular domains was linked to the plasticity of the binding groove, we generated a subset of the original simulation trajectories without HLA-DM (-DM) that excluded states with globular domains that had separated from each other and the binding groove (distance between centers of mass of globular domains α 2 and β 2 ≥ 26 Å and nearest distance between the binding groove and the globular domain ≥6 Å) (-DM*). These cutoffs were chosen based on the distance distributions observed during simulation ( Figure S2). Figure S3 shows the distance between the binding groove and the globular domain that was measured. Interestingly, the resulting conditional distribution matched the binding-groove width distribution seen in the simulation where HLA-DM was present (+DM in Figure 2C).
To confirm this result, we designed a customized Natural Move simulation in which the two globular domains α2 and β2 were grouped (using original positions observed in the crystal structure) -DM**, thereby effectively preventing the separation of the globular domains throughout the simulation (unlike in the -DM simulation, where the globular domains were allowed to move independently). Grouping the globular domains was sufficient to prevent the collapse of the binding groove in the absence of HLA-DM as shown by the -DM** distribution in Figure 2C. The results closely matched the distributions of the simulation where HLA-DM was bound (+DM) and the subset of the simulation with freely moving globular domains that only included interacting globular domains (-DM*). Figure 3 shows that the effect we discovered cannot be traced back to a simple interaction between the top or bottom part of HLA-DM and MHCII. The complete structure of HLA-DM is required for the stabilizing effect to be seen. Note that the globular domains of MHCII are less mobile in the presence of the bottom part of HLA-DM than in the presence of the top part of HLA-DM. However, neither fully restricts the movement of the MHCII globular domains in our simulations; this only happens in the presence of the full HLA-DM structure. Note, we ran initial simulations of the MHCII/HLA-DM complex with a mutated residue in the βchain globular domain (B:V186K), and while binding of the globular domain was affected, the effect was not sufficient to reproduce the binding-groove collapse (data not shown). A more detailed visualization of the HLA-DM/MHCII complex is shown in Figure S4.
This suggests a long-range mechanism that implicates the role of the MHCII globular domains in the modulation of the MHCII binding groove (Figure 4). When the peptide-loading chaperone HLA-DM is bound, the MHCII globular domains are kept in a compact configuration, and the binding groove remains in an open state. In the absence of HLA-DM, the MHCII globular domains are free to move, which allows increased movement of the β1 helix leading to an increased collapse of the binding groove.
We tested the hypothesis that binding-groove plasticity may be linked to the configuration of the globular domains. We performed the following simulations: 1. Simulation of MHCII in the absence of HLA-DM with independently moving globular domains; 2. Simulation of MHCII in the presence of HLA-DM with independently moving globular domains; and 3. Customized Natural Move simulations of MHCII in the absence of HLA-DM, in which the globular domains were propagated as a unified segment.
In our initial simulations without HLA-DM and independently moving globular domains, the globular domains regularly visited structural states, which showed increased occurrence of a narrowed binding groove and coincidentally separated globular domains. We used customized Natural Moves to investigate the importance of the globular domains in the collapse of the binding groove. Using cNMMC, we propagated the globular domains collectively rather than independently, which prevented globular domain separation and bindinggroove collapse. When the simulation with independent globular domains (and no HLA-DM) was filtered for states with compact globular domains, the proportion of open binding-groove states was also significantly increased, suggesting a dependency between the configuration of the globular domains and the conformational state of the binding groove.

■ DISCUSSION
The experimental observations from the literature together with our results suggest a mechanism where, in the absence of HLA-DM, the MHCII globular domains change conformation, resulting in a long-range effect which causes the binding groove to collapse. Specifically, it seems that the β1-1 kink on the β1 helix plays a major role in this plasticity. HLA-DM prevents this collapse indirectly by restraining the globular domains in a compact conformation and thereby stiffening the β1 helix and stabilizing the open binding groove.
In our simulations we found that HLA-DM modulates the receptiveness of the MHC binding groove indirectly via a necessary contact to the globular α2/β2 regions of the MHC. This finding is supported by several lines of experimental evidence. Carven et al. have used chemical labeling and mass spectrometry to characterize residues that are involved in peptide-induced configurational changes in the globular domains. 20 Furthermore, studies comparing MHCII crystal structures have shown conformational diversity in the globular domains, specifically rigid body motions of the β2 domain of up to 15 degrees. 21,22 HLA-DM binds MHCII laterally, with the main interaction sites being at the N-terminal side of the α1 helix and on the two globular domains, 7 distal from the structurally labile β1 helix. 20 Several groups have observed increased structural movement in the β1 helix in simulations of empty MHCII molecules. [10][11][12]22,23 The question arises how HLA-DM may counteract this structural lability if its binding site is on the opposite of the MHCII molecule.
In 2004, Carven et al. generated four antibodies that detected peptide-induced structural changes. Three of these conformation-sensitive antibodies bound to an epitope on the β1 helix (β53-67), and the fourth antibody, MEM-266, bound to the last five residues on the β2 globular subunit. 10 As the most important epitope residues are positioned in a linear contiguous fashion and similar MEM-266 affinities were observed for the empty protein and the corresponding epitope peptides, it was suggested that the epitope was nonconformational. Therefore, selective antibody binding was most likely regulated by configurational changes of the globular domains rather than local residue rearrangements.
In 2000, before the first MHCII/HLA-DM structure was published, Doebele et al. found a set of residues on MHC that, if mutated, led to the disruption of DM activity. Among this set was a mutation in the β2 globular domain (V186 K), which is distant from the binding groove and is one of the sites that make up the epitope region for MEM-266.
Given this evidence Carven et al. suggested a HLA-DM chaperoning mechanism involving the long-range transmission of structural changes from the globular subunits to the binding groove. Specifically, they proposed that HLA-DM binding causes structural changes in the membrane-proximal β2 subunit of MHCII that is propagated to the binding groove to modulate peptide binding. 10 The results presented in this study, that is the stabilization of the MHCII binding groove via a long-range interaction involving the globular domains of both MHCII and HLA-DM, corroborate the observations described above. Furthermore, our results point toward a mechanism that could explain how the stabilization takes place despite the main site of interaction between HLA-DM and MHCII being at the opposite end of the labile MHCII β-chain helix.

■ CONCLUSION
The biophysical mechanism, by which HLA-DM stabilizes the peptide-receptive state of the MHCII binding groove, is a subject of debate. Experimental results in the literature suggest that the mechanism involves long-range structural changes that propagate from the membrane-proximal globular subunits to the membrane-distal binding groove. In our Natural Move Monte Carlo simulations we were able to reproduce the MHCII binding-groove stabilization by HLA-DM. By running simulations with various structural parts of HLA-DM, we have found that neither the globular domains nor the membranedistal part that interacts with the MHCII binding groove stabilizes the MHCII binding groove in isolation. The full HLA-DM structure was required to stabilize the MHCII globular domains and by extension the MHCII binding groove. This result corroborates existing hypotheses about the nonintuitive stabilization of the peptide-receptive state of the MHCII binding groove.

■ METHODS
Natural Move Monte Carlo. In Natural Move Monte Carlo (NMMC), the system is decomposed into a priori defined structural segments or groups of segments, which are called regions. The set of residues that are not part of segments constitutes the molten zones. Within the current model (see below) a molten zone is a continuous set of coarse-grained residues connecting adjacent segments of the same molecular chain.
Formally, all degrees of freedom, X, are partitioned into independent (X i ) and dependent (X d ) degrees of freedom (DoF). Here, X i represents independent translational and orientational arrangement of structural segments or regions (groups of segments). In addition, X i also contains the internal DoFs of segments, such as torsional and bond angles. The independent motion of segments (e.g., within one chain) may lead to chain breaks, which are repaired using a chain closure algorithm 24 that adjusts the position of coarse-grained centers of molten zone residues by changing corresponding torsional and bond angles, which are the dependent DoFs, X d . In this way X d are instantaneously adjusted (following each proposed change in X i ) to facilitate exploration along X i and simultaneously preserve the integrity of the molecular chain(s) through chain closure(s).
The basic principle is that each new configuration during a proposal step is obtained via a combined chain breakage closure algorithm. This composite proposal kernel includes a stochastic proposal to update X i followed by finding an optimal (with respect to the new X i ) arrangement along X d so that all chain breaks are repaired. More details about the method can be found in the following publication. 24 For the current model used, the exact definition of independent X i and dependent X d degrees of freedom is described in detail. 24 In addition, we provide an online tutorial with movies (www.cs.ox.ac.uk/mosaics) that illustrates the independent degrees of freedom in each simulation.
Note that the algorithm aims to sample the conformational space along user defined independent degrees of freedom as described in ref 12. Given this initial choice, the method generates canonical distributions over an effective energy

Journal of Chemical Information and Modeling
Article surface. Since the proposal kernel is symmetrical, we use classical Metropolis Monte Carlo, which satisfies detailed balance, to sample the different states.
Customized Natural Move Monte Carlo. Traditionally NMMC is used to explore the conformational landscape along a particular set of degrees of freedom chosen by the researcher. However, the initial choice of degrees of freedom might not always be optimal. Additionally, if the objective is to investigate the causality of functional motions, it may be informative to perform NMMC simulations for a variety of sets of degrees of freedom. NMMC allows for the modulation of translations and rotations of segments as well as torsion and bend angles of bonds. Thus, users can compare different sets of customized Natural Move Monte Carlo (cNMMC) simulations to infer causal relationships in functional motions. cNMMC was introduced in ref 12 as a practical framework for systematically designing a set of said simulations that together are targeted at probing a specific biological question. In this way the cNMMC approach entails an ensemble of NMMC simulations each with a unique set of degrees of freedom aimed at testing a specific hypothesis.
Practical Implementation of the Decomposition. Natural Move Monte Carlo requires the decomposition of the structure into segments and flanking molten zones. The segments are propagated at each MC iteration, and the molten zones serve to close the chain break. Thus, the molten zones represent labile areas in the structure and should be carefully chosen. In the MHCII structure, residues β53-68 on helix β1 are part of epitopes for conformation-sensitive antibodies that are selective for the empty binding groove. 10,25 This region has been shown to undergo local structural changes by CD spectroscopy. 10 MD simulations and comparison of experimental MHCII structures revealed structural variability around a sharp kink in this region. 11,22,23 Furthermore, we have shown using customized Nature Moves that this region contributes largely to the collapse of the binding groove. Given these observations, we introduced a molten zone at the Nterminal kink of the β1 helix (β1-1 in Figure 2B), as well as at the regions flanking the β1 helix.
The Model. The coarse-grained model has three centers per residue. For each residue, the three atoms chosen are the C, the carbonyl oxygen, and a single side-chain atom that is the closest to the center of mass of the side chain. The three center per residue model inherits an interaction potential from previous work, a carefully parametrized knowledge-based interaction potential 26 trained on 4500 known protein structures. Our model (the current three center per residue model combined with the above interaction potential) has been benchmarked to reliably describe the nativelike properties of most known protein folds. 27 Our knowledge-based potential (similar to most statistical potentials trained on experimental structures) implicitly incorporates the effects of solvent. Furthermore, the current model (3-point per residue coarse-grained description combined with the above-mentioned knowledge-based interaction potential) has been shown to reproduce functionally relevant and often experimentally seen molecular motions in several applications on systems including molecular chaperonines, 18 diabodies, 14 and even MHC complexes. 12,19 Simulation Details. All simulations were carried out with the MOSAICS software package. 28 All distributions were plotted with matplotlib 29 and pandas 30 using a bandwidth of 0.1. NMMC simulations were initiated from an X-ray structure of HLA-DR (MHCII) in complex with HLA-DM at a resolution of 2.6 Å (PDB: 4GBX). 7 The peptide was removed, and the MHCII/HLA-DM complex was equilibrated for 10,000 steps at 300 K. The structures were coarse-grained using a 3-point per residue protein model. 27 The MHCII model was generated by removing the HLA-DM part of the structure file and equilibrating for a further 10,000 iterations at 300 K. In order to ensure extensive conformational sampling, we performed Parallel Tempering using six replicas at temperatures 300, 336, 376, 421, 472, and 529 K. We ran 15 independent repeats for each test case. Each repeat was run for 6 × 1,000,000 Monte Carlo iterations each, and the replica exchange rate was 0.1. The average acceptance rate within replicas was 0.5 and 0.3 for MHCII and MHCII/HLA-DM simulations, respectively. The average acceptance rate for jumps between replicas was 0.16. All data were collected at a canonical temperature of 300 K. The average run time per simulation was 16 h (1 h equilibration, 15 h production) using individual Intel E5-2640v3 Haswell cores on the ARC Oxford Supercomputer. Distances were calculated with MDAnalysis. 31 Structural Distances. In the analysis of our simulations we define the binding-groove width as the distance between the centers of mass of residues chain A:60-65 and chain B:65-70 (as defined by Kumar et al. 32 ). The distances between the two globular domains α2 (A:84-182) and β2 (B:  were calculated between the residues closest to their respective centers of mass. The distance between globular domain β2 and the binding groove was calculated between the two residues on chain A and chain B that were closest in the starting structure (A:29 and B:151, see Figure S3).

Journal of Chemical Information and Modeling
Article