Selector function of MHC I molecules is determined by protein plasticity

The selection of peptides for presentation at the surface of most nucleated cells by major histocompatibility complex class I molecules (MHC I) is crucial to the immune response in vertebrates. However, the mechanisms of the rapid selection of high affinity peptides by MHC I from amongst thousands of mostly low affinity peptides are not well understood. We developed computational systems models encoding distinct mechanistic hypotheses for two molecules, HLA-B*44:02 (B*4402) and HLA-B*44:05 (B*4405), which differ by a single residue yet lie at opposite ends of the spectrum in their intrinsic ability to select high affinity peptides. We used in vivo biochemical data to infer that a conformational intermediate of MHC I is significant for peptide selection. We used molecular dynamics simulations to show that peptide selector function correlates with protein plasticity, and confirmed this experimentally by altering the plasticity of MHC I with a single point mutation, which altered in vivo selector function in a predictable way. Finally, we investigated the mechanisms by which the co-factor tapasin influences MHC I plasticity. We propose that tapasin modulates MHC I plasticity by dynamically coupling the peptide binding region and α3 domain of MHC I allosterically, resulting in enhanced peptide selector function.


S1.1 Atomistic molecular dynamics simulation protocol
All the molecular dynamic simulations were performed using the IRIDIS High Performance Computing Facility, and we acknowledge the associated support services at the University of Southampton in the completion of this work.
The starting conformations are experimentally determined structures from the RSCB Protein Databank. These were the X-ray crystal structures of HLA-B*44:02, PDB id: 1M6O and HLA-B*44:05, PDB id: 1SYV. The in silico W147A point mutation to the structure of HLA-B*44:05 was performed using MODELLER (1).
The GROMACS version 4.5.3 (2, 3) molecular dynamics package was used for the all atom simulations. The simulations used the Amber99SB-ILDN (4) force field and TIP3P (5) explicit water molecules using the Simple Point Charge water system (6), and Sodium counter ions were added to neutralise the charge of the system. The protein structures were placed in rhombic dodecahedron shaped box centred at 1.5 nm from the edge with periodic boundary conditions. Covalent bond lengths were constrained using the P-LINCS algorithm (7) and the water angles were constrained using the SETTLE algorithm (8) allowing an integration time step of 2 fs to be used.
Nosé-Hoover temperature coupling (9,10) and Parinello-Rhaman pressure coupling (11,12) used a time constant of 0.5 ps with reference baths of 300 Kelvin and 1 bar respectively to maintain the average thermodynamic properties of the protein and solvent comprising the system. Electrostatic interactions use a cut-off of 1 nm with the interactions beyond this cut-off treated using the particle mesh Ewald method (13). Van der Waals forces used a cut-off of 1 nm. The neighbour list is updated every five steps.
Each system initially underwent an energy minimization over 1000 steps of 2 fs to relax the structure and remove the forces from the systems that were introduced by the protonation of the molecule and addition of solvent. This was followed by a 5 ns equilibration of the water surround the protein with the protein atoms restrained using a randomly generated initial starting velocity.
Full production runs were performed with the position restraints released. To analyse conformational dynamics, concatenated trajectories of 420 ns were created from three independent repeats of 150 ns, with the first 10 ns of each simulation discarded. Two additional control simulations were performed over 45 ns each.
For the simulations using distance restraints, a simple force constant of 7.437 kJ mol −1 nm −2 was applied to the selected C α atoms for the restraint. This is equivalent to 3 k b T, where k b is Boltzmann's constant and T is temperature in Kelvin. For the control simulations we used a strong restraint of 100 kJ mol −1 nm −2 equivalent to 40 k b T.
The system components are summarised in Table S1. Parameter files for equilibration and the initial production runs were kindly provided by Tom Piggot (University of Southampton, UK) along with his assistance in the installation and invaluable advice in using GROMACS. Quality assurance and post processing was performed using a combination of the suite of utilities provided with GROMACS. In some cases these utilities have been adapted as described in the text where relevant. Additional post-processing tasks were performed using MATLAB™ and bespoke UNIX awk scripts. Visualisation of the protein structures and molecular dynamics trajectories was performed using the VMD (14) and USCF Chimera (15)

S1.2 Simulation stability
The stability and convergence of the molecular dynamics simulations towards an equlibrium state was assessed by calculating the Root Mean Square Fluctuation for blocks of each trajectory as plotted in Figure  S1. These indicate that the RMSF between blocks is less than 1 Å and that all the simulations were stable.

S1.3 Global motions
Global motions of the molecular dynamics simulations were analysed in three ways: Normalised Covariance Analysis, Functional Mode Analysis and Probability Density.

S1.3.1 Normalised covariance analysis
We analysed the degree of correlation between the motions of pairs of atoms of the MHC I complex over the course of the combined simulation trajectory. Highly correlated atoms are moving together and this therefore this indicates they are forming more a rigid body-like structure. A mass weighted variance-covariance matrix was built using the C α atoms. This is a symmetric 3N × 3N matrix comprising of the fluctuation of the atom positions with coordinates x as a function of the trajectory such that: where <> indicates the conformational ensemble average. This matrix C therefore contains as elements, for each atom pair, the difference between the mean product of their atomic positions and the product of their mean atom positions i.e. the difference between their average position as a pair and the product of their individual average positions. Atom pairs moving together in the same direction give rise to positive covariances and pairs moving in the opposite direction give rise to negative covariances. Non-correlated atoms give near zero covariances. The variance for each atom is contained on the main diagonal. We calculate the normalized covariance of the atoms by summing the x, y and z components of the matrix and normalizing with the self-covariance of the atoms. Re-expressing equation 1 for atoms i and j as: the normalized covariance is: This yields a matrix containing the correlations between the atoms in terms of the Pearson Correlation Coefficient for each pair of atoms. Extracting atom pairs above a given magnitude of correlation e.g. greater than 0.7, it is possible to create an image of the correlated atoms on the three dimensional structure as a web of connections (16). This gives an indication of which atoms are moving together during the simulation which in turn indicates parts of the structure that are acting as rigid bodies or where there is potential communication between parts of the structure. It does not however indicate the magnitude of the motions, only that there are correlated above a certain threshold i.e. there is a linear relationship between the motion of correlated atoms as defined by the strength of correlation coefficient chosen. Therefore this analysis cannot tell us anything about any non-linear relationships between the atoms of the protein. For this analysis we used an amended version of the GROMACS g_covar utility and the g_anaeig utility. We used a common peptide free reference structure with the combined trajectories for the analysis. This therefore excludes correlations with the peptide in the peptide bound state. Fig. S2 shows the web of correlations between all pairs of C α atoms with a normalised covariance of greater than 0.7 for the peptide bound and peptide free HLA-B44 simulations. It's important to reiterate this is simply a dimensionless number indicating the degree above which we observe a correlation between a pair of atoms and that the choice of 0.7 is some what arbitrary. The choice of C α atoms is simply to provide clarity in visualization of the correlations.  Normalised covariance webs C α atom pairs with a covariance yielding a correlation coefficient greater than 0.7 are shown connected by red lines for the peptide bound and peptide free simulations. The red lines indicate that connected atoms move together during the simulation, but not the direction or magnitude of the motion. Correlations with the peptide have been excluded for clarity. The calculations were done using an amended version of the GROMACS g_covar utility and the g_anaeig utility.

S1.3.2 Functional mode analysis
Having performed Principal Component Analysis using the GROMACS g_covar and the g_anaeig utilities, as previously described (17), to extract top 50 modes accounting for ∼90% of the total atomic motion, we performed Functional Mode Analysis (FMA) of peptide free MHC, as proposed and implemented by Jochen Hub and Bert de Groot, . It aims to detect the collective atomic motions directly connected to protein function. It therefore requires the selection of a "functional quantity" that describes the functional state of the protein to be correlated with the protein motion. A full description of the theory and the method is detailed in the publication by Hub (18).
In the absence of experimental data to guide the choice of a functional quality, we make the observation that in the X-ray crystal structure the peptide is buried inside the peptide binding groove with no easy means of entry or exit. Yet we know that peptides load and can be exchanged. We identified a potential hinge in the α 2−1 helix and measured the distance fluctuations in this region to define a functional quantity for peptide binding and unbinding. It only makes sense to consider this for peptide free MHC as in the peptide bound state there are no significant distance fluctuations.
Correlating the changes in distance with the overall motion of the molecule as described by the principal components yields a single collective motion for MHC I associated with the conformational changes in the peptide binding groove. The method is implemented using the common reference structure created for principal component analysis and the FMA package developed by Jochen Hub (18): 1. From identification of the two hinge points common in the α 2−1 helix, two regions were defined on either side of the F-pocket of the MHC molecule peptide binding groove. Residues 135 to 156 on the α 2−1 helix and residues 69 to 85 on α 1 helix as coloured red in Fig. S3 2. For each combined trajectory in the peptide bound and peptide free states the distance was measured between the centres of mass of these two regions over the duration of the concatenated simulation. This distance is a function of time d α 1 α 2 and is the functional quantity.
3. Assuming d α 1 α 2 is approximately a linear function of the principal components, a collective vector a was constructed from the principal components derived by diagonalization of the covariance matrix.
As the first 50 principal components account for ∼90% of the atomic mean square atomic fluctuations, this subset of the principal components was considered a reasonable set for the construction of a.

4.
Quantifying the correlation using the Pearson correlation coefficient R, the motion along a is maximally correlated to the change in d α 1 α 2 to yield the Maximally Correlated Motion (MCM) as function of time given by projection p a , such that: 5. The process of maximizing the R generates a model for d α 1 α 2 as a function of the principal components. The model is cross validated (Fig. S3) by dividing the trajectory into model building frames and cross validation frames. We used 350 ns for model building and 70 ns for cross validation. We test the ability of the model to predict the value of R in the cross validation set as compared with that calculated from the data (Fig. S3). This was analysed to determine contributions from the different principal components to the variance of the model and therefore the influence of individual principal components on d α 1 α 2 (Fig. S3) 6. The MCM was then further optimized to find the most probable motion given the input ensemble of structures. This yields the ensemble-weighted Maximally Correlated Motion (ewMCM) in accordance with the free energy landscape described by the input trajectory (18). This estimates the most probable collective motion for the MHC Class I molecule that achieves a substantial change in d α 1 α 2 which is visualised as a projection of the trajectory onto this eigenvector (Fig. S3).

S1.3.3 Probability density analysis
Probability densities were calculated to characterise the conformations explored by the MHC structures in two ways, one using intra-molecular distances and the other combining intra-molecular distance and angle between the heavy chain domains. Intra-molecular distances were chosen from the identification of the two common hinge points in the α 2 helix from the conformational angles analysis, two regions are defined on either side of the F-pocket of the MHC molecule peptide binding groove. Residues 135 to 156 on the α 2 helix and residues 69 to 85 on α 1 helix as coloured red in Fig. S4A. The domain-domain distance was defined as between residues 96-100 in the peptide binding groove platform and the flexible α 3 domain region 220-227, also coloured red in Fig. S4A. These distances were measured for the combined simulations using the GROMACS utility g_dist and joint probability densities for these two distances were then calculated using MATLAB ™ as plotted in Figures 4 and 5 in the main manuscript.
Restrained control region 1 in blue Restrained control region 2 in blue

S1.3.4 Interpretation of twist angle
For the angle of twist between peptide binding domain and the α 3 two planes represented as discs in Fig. S5A were defined through the centre of mass of the heavy chain by the C α atoms of residues 118, 174 and 252 and the centre of mass of the α 3 domain by residues 199, 209 and 260. The angle between the normal to these planes θ during the 420 ns combined simulations measures the twisting angle in degrees. This angle was calculated using GROMACS utility g_sgangle. The range of this twisting angle for each HLA-B*44 molecule is shown in Figure 1F of the main text. Joint probability densities for the F-pocket distance defined in Fig. S4A and the domain-domain twisting were then calculated using MATLAB ™ as plotted in Fig. S5C-E.

One conformation model
Previously, we constructed a model of MHC class I that could describe the time-dependent optimisation of peptide cargo (19). Emanating from this work was the prediction that HLA-B molecules vary in their intrinsic ability to load high affinity peptides, which influences the extent to which tapasin confers an additional optimisation benefit via skewing the competition between tapasin binding and peptide binding.
In this study, we sought to determine whether this model could be extended to incorporate a conformational intermediate, and understand more specifically how peptide loading might be altered in different HLA-B molecules. Before defining the two conformations models, we introduce the original oneconformation model, detailing some alterations that were made and propagated to the two-conformations models.
1. As peptide supply and turnover is much faster than the other kinetics in the system, we assume that the concentration of peptide is in equilibrium. i.e. [P i ] ≡ P i .
2. After collecting temporally resolved data for the three HLA-B molecules investigated herein, we noticed that the abundance of recoverable radioactive molecules sometimes increased between the data-points at 15 minutes and 30 minutes after the termination of radio-labelling. Therefore, on this timescale, some process must occur that changes MHC I molecules from an immature state to a mature state (in terms of its recognition by antibody). This could be the effect of β 2 m binding, calnexin unbinding, both, or neither. To account for the effect (which could be explored in detail in another study), we modelled this as a first order reaction in which unrecoverable molecules (M u ) transition into mature molecules (M) with rate m, following endoplasmic reticulum (ER)entry of immature MHC I molecules at rate g M .
3. Changing the v T reaction to reversible, with rate a T , Applying these two alterations to the original model gives rise to Assuming mass action, the corresponding ODEs are given by Solving the system at steady state, with a T = 0, leads to the following expression for the populations of egressed peptide-MHC I complexes, as previously derived (19) [

S2.1.2 Two-conformations (unbinding) model
To incorporate the possibility of a conformational change that alters peptide binding and unbinding, we first proposed a model that extends the one-conformation model simply by using two conformational states for all MHC I molecules ( Fig. 1B and S6). These are distinguished as molecules that are able to load peptide ("open") and those that are not able to load peptide (i.e. "closed"). Transitions between open and closed states proceed by opening rates o and o e for peptide-bound and empty molecules respectively, and closing rates c and c T for tapasin-free and tapasin-bound molecules respectively. Additionally, we assumed that MHC I molecules transition to the closed state (M c ) immediately following maturation from the state M u . As we explicitly model a maturation process during which time no peptide binding can occur, we felt that this assumption would lead to very minor differences to the alternative assumption of M u → M o . As for the one-conformation model, we allowed tapasin to enhance peptide dissociation by a factor q. The underlying reactions for the two-conformations (unbinding) model are given by (see Fig. S6a for a graphical depiction)

Relationship to the one-conformation model
In order to determine the consequences of the additional conformational state, we derive the one-conformation model from these equations. Intuition suggests that the one-conformation model is approximately equivalent to the closing and opening rates being infinitely fast. Therefore, we apply the following substitutions: which is precisely equivalent to the one-conformation model with By observing these relationships, we can interpret the tapasin enhancement of the peptide on-rate as and the Tapasin enhancement of the peptide off-rate aŝ This enables us to re-interpret the filter principle in the two-conformations model, assuming that closing and opening occur on a faster timescale than peptide binding/unbinding (note that for comparison with the one-conformation model, we assume a T = 0 in the following). This gives where e = e · c+o o and x = v T q · c T +o o . Note that this equation has the same form as for the one-conformation model (6), so the filtering relation is unchanged for the two-conformations (unbinding) model.

S2.1.3 Two-conformations (opening) model
Based on an argument related to thermodynamics, it is more likely that the peptide-dependent step of peptide-MHC I unbinding is the rate at which a "closed" molecule becomes "open". i.e. rate o in the twoconformations model. This argument rests on the observation that the extent of the interactions a peptide makes with a MHC I molecule would be greater in a closed state than an open MHC I state and therefore the number and quality of these interactions would determine the opening rate of MHC I in a peptidedependent manner. Therefore, we propose a variation to the two-conformations model that is as close as possible to the other models, except that the rate of opening of a peptide-bound MHC I molecule is peptidedependent, while the rate of peptide-MHC I disassociation is homogeneously assigned. We propose that the effects of tapasin are to increase peptide unloading by increasing the peptide-dependent opening rate (o i → q · o i ), and to modify the rate of closing (c → c T ), though this time having no effect on peptide binding (rate b) and unbinding (rate u).
The underlying reactions for the two-conformations (opening) model are given by (see Fig. S6b for a graphical representation) Using exactly the same procedure as above, we obtain the one-conformation model with Therefore, the tapasin enhancement of the peptide on-rate remains the same as for the two-conformations (unbinding) model as but the tapasin enhancement of the peptide off-rate becomes peptide-specific as We now derive an expression equivalent to (6) and (9) but cast in terms of the peptide-specific opening rate o i to quantify the filtering achieved by this mechanism. We obtain where x 1 = e u+e , ρ = v T v T +u (< 1) and x 2 = c T q .
Using the relationships between the unbinding, closing and opening rates derived in this document, we can plot representative curves that compare the effective dissociation in each model (Fig. S6C). The twoconformations (unbinding) model exhibits the behaviour assumed in the one-conformation model, with tapasin increasing peptide unbinding by a constant factor. A more complex relationship was observed with the two-conformations (opening) model, with tapasin having no effect for higher values of o i (and therefore less stable peptides). Given this fundamental difference between the two variants of the two-conformations model, it would not be unexpected that different dynamical behaviours could be observed. Therefore, we reasoned that the models would differ in their ability to reproduce experimental data when attempting to fit the underlying model parameters.

S2.2 Simulation of thermostability and endoglycosidase H resistance
All simulations were performed as done previously (19). The one-conformation and two-conformations (unbinding) models used representative high, medium and low affinity peptides, instantiated as three different peptide off-rates u i . For the two-conformations (opening) model, representative peptides were instantiated with peptides for which peptide-MHC I opening o i differed. To relate the simulations to the experimental measurements, 50 • C data were compared with the high affinity peptide alone, 37 • C data were compared with the sum of the high and medium affinity peptide-MHC I complexes, while 4 • C data were compared with the total quantity of MHC I molecules in the system.
For endoH resistance, we considered the total quantity of egressed MHC I molecules to be resistant, and the intra-ER MHC I molecules to be susceptible. Therefore, to calculate % endoH resistance, we used the formula endoH = 100 · resistant susceptible + resistant

S2.3 Bayesian parameter inference for the kinetic models
To assess the plausibility of a model of a specific circuit, we first determine optimal parameter values using probabilistic inference techniques, similar to the technique used previously (19). In particular, we seek to approximate the likelihood of parameters taking on specific values, given a model hypothesis and some observation data. i.e. we attempt to approximate the posterior density Pr(θ|H, D), where θ is the vector of parameters to be inferred, H is the model hypothesis, and D is the set of experimental data used for inference. The posterior distribution is related to an evidence or likelihood function Pr(D|θ, H) according to Bayes' rule: Pr(θ|H, D) = Pr(D|θ, H)Pr(θ|H) Pr(D) We obtained approximations of the posterior distributions using the Filzbach software, available from the author's website (http://research.microsoft.com/science/tools), which uses a Metropolis-Hastings (MH) Markov Chain Monte Carlo (MCMC) sampling routine. MCMC is a stochastic search strategy that forms a Markov chain of proposal parameter vectors, moving to new proposal vectors based on the ratio of the likelihoods of the proposal and previous points. By biasing the stochastic search in parameter regions of high probability mass, we converge on the true joint posterior distribution of the parameters more efficiently than would be possible with a purely random search.
To define the likelihood of a parameter vector, we assumed that the experimental measurements were noisy samples drawn independently from Gaussian distributions centered on the model predictions of the fluorescence. Therefore, in the optimal case that the model precisely describes the underlying biological behaviour, deviations of the measurements result purely from experimental error. Consequently, a data point y i is distributed as y k ∼ N(x k , σ 2 ), where x k is the model prediction at t = t k and σ 2 is the variance of experimental error. We assume that the variance of experimental error is proportional to the measured fluorescence signal, as quantified in Figs. 2B, 5B (i.e. σ = α √ y k for some α). Therefore, the likelihood function is formed as the product of the probabilities of each data-point with α left as a parameter to be inferred. As the model simulates the concentration of peptide-MHC I complexes and not fluorescence directly, we made the simplifying assumption that there is a linear relationship between the two. Therefore, we compared the simulated output to the data with a scale factor, which was calculated using linear regression.
We determined the extent to which each model hypothesis (one-conformation, two-conformations (unbinding) and two-conformations (opening)) could reproduce experiments measuring time-dependent optimisation via thermostability (Fig. 6A) and cell surface transit via endoglycosidase H resistance (Fig. 6B). As before (19), we examined a variety of hypotheses for which parameters could be nominated as allele-specific, then used our parameter inference procedure to find optimal values. We then calculated the Bayesian Information Criterion (BIC), which attempts to find a compromise between having as few variable parameters as possible while minimising the deviation between model simulation and experimental observation. The BIC is defined as BIC := −2 log L (θ * ) + k ln (n) where θ * is the maximum likelihood estimate of the parameters (i.e the vector that maximises L(θ)), k is the number of variable parameters and n is the number of experimental observations.