A Peptide Filtering Relation Quantifies MHC Class I Peptide Optimization

Major Histocompatibility Complex (MHC) class I molecules enable cytotoxic T lymphocytes to destroy virus-infected or cancerous cells, thereby preventing disease progression. MHC class I molecules provide a snapshot of the contents of a cell by binding to protein fragments arising from intracellular protein turnover and presenting these fragments at the cell surface. Competing fragments (peptides) are selected for cell-surface presentation on the basis of their ability to form a stable complex with MHC class I, by a process known as peptide optimization. A better understanding of the optimization process is important for our understanding of immunodominance, the predominance of some T lymphocyte specificities over others, which can determine the efficacy of an immune response, the danger of immune evasion, and the success of vaccination strategies. In this paper we present a dynamical systems model of peptide optimization by MHC class I. We incorporate the chaperone molecule tapasin, which has been shown to enhance peptide optimization to different extents for different MHC class I alleles. Using a combination of published and novel experimental data to parameterize the model, we arrive at a relation of peptide filtering, which quantifies peptide optimization as a function of peptide supply and peptide unbinding rates. From this relation, we find that tapasin enhances peptide unbinding to improve peptide optimization without significantly delaying the transit of MHC to the cell surface, and differences in peptide optimization across MHC class I alleles can be explained by allele-specific differences in peptide binding. Importantly, our filtering relation may be used to dynamically predict the cell surface abundance of any number of competing peptides by MHC class I alleles, providing a quantitative basis to investigate viral infection or disease at the cellular level. We exemplify this by simulating optimization of the distribution of peptides derived from Human Immunodeficiency Virus Gag-Pol polyprotein.


Introduction
MHC class I molecules are encoded within the genetic region known as the Major Histocompatibility Complex and are present in all nucleated human cells. MHC class I molecules direct cytotoxic T lymphocytes (CTL) to destroy virus-infected or cancerous cells, thereby preventing disease progression [1]. They provide a snapshot of the internal contents of a cell by binding to peptides arising from intracellular protein turnover and presenting these peptides at the cell surface, where the peptide-MHC complex can be recognized by CTL (Fig. 1). Most cells will present an array of tens of thousands of different peptides at their surface, some of which will be unique to virus-infected or cancerous cells. The efficacy of a CTL response to these peptides depends to a large extent on the ability of MHC class I molecules to select only a limited number of the potentially billions of different peptides that are generated by the hydrolysis of all intracellular proteins [2]. Peptides are selected for presentation on the basis of their ability to form a stable complex with MHC class I, by a process known as peptide optimization. A better understanding of the optimization of peptides is important for our understanding of immunodominance [1], the predominance of some CTL specificities over others, which can determine the efficacy of an immune response, the danger of immune evasion, and the success of vaccination and immunotherapeutic strategies.
Peptide binding to MHC class I occurs in the endoplasmic reticulum (ER) and is assisted by multiple cofactors that are thought to enable the optimization process [3]. The transporter associated with antigen processing (TAP) supplies the lumen of the ER with peptides generated in the cytosol and forms the backbone of the peptide loading complex (PLC). A number of chaperone molecules also comprise the PLC, namely calreticulin, calnexin, ERp57 and importantly tapasin, which bridges the gap between MHC class I and TAP [4] (Fig. 1). Of these chaperones, only tapasin is known to influence the extent of peptide optimization, in such a way as to skew the cell surface cargo towards peptides with low off-rates [5], and this influence is known to vary across different MHC class I alleles [6].
While a range of interactions between tapasin and MHC class I have previously been identified [7,8], the effects of these interactions on peptide optimization are still not well-understood. A recent study used computational modeling to distinguish between different hypotheses of tapasin function within the ER [9], but the model assumed that only peptides with low off-rates could egress to the cell-surface, and was therefore unable to predict the optimization of peptides with different off-rates. As a result, the model was unable to account for observed effects of tapasin on peptide optimization, both over time [6] and at steady state [5].
In this paper we present a dynamical systems model for predicting MHC class I peptide optimization. We include interactions with the chaperone molecule tapasin, and propose a relation of peptide filtering to quantify peptide optimization as a

Author Summary
Major Histocompatibility Complex (MHC) class I molecules bind to protein fragments (peptides) within the cell and present these fragments at the cell surface, thus providing a snapshot of the cell contents that can subsequently be used to trigger an immune response. Only a fraction of the potentially billions of peptides inside a cell are selected for presentation, and the process is optimized to select for peptides that form a stable complex with MHC class I. The mechanisms of the optimization process are important for predicting the efficacy of an immune response and for designing effective vaccines, yet are still not wellunderstood. In this article we present a dynamical systems model of peptide optimization by MHC class I. We show that peptide optimization can be quantified and mechanistically explained by a peptide filtering relation, which relates cell surface abundance to peptide supply, peptide unbinding and interactions with the chaperone molecule tapasin. The filtering relation also accounts for differences in optimization across MHC class I alleles. Finally, we show how the filtering relation can be used to quantify the cellsurface presentation of virus-derived peptides for immune system surveillance.
function of peptide supply and peptide off-rates. Using a combination of published and novel experimental data, together with a combination of Bayesian inference and kinetic analysis, we show that tapasin can improve both the rate and extent of peptide optimization by accelerating peptide off-rate, and that differences in optimization across MHC class I alleles can be explained by an allele-specific peptide on-rate. Our filtering relation provides a mechanistic interpretation for recent experimental observations of peptide optimization both over time [6] and at steady state [5]. Finally, we demonstrate how the filtering relation can be used to quantify optimization of a large set of competing peptides in the context of an immune response, by simulating the cell surface abundance of Human Immunodeficiency Virus (HIV) peptides in complex with different MHC class I alleles.

A model of MHC class I peptide optimization
We formulated a dynamical systems model of MHC class I peptide optimization using a biological modeling language (SPiM [10], Fig. S1 in Text S1) and exported the model to an equivalent set of biochemical reactions for further analysis (Fig. 2). The model characterizes the interactions between MHC, peptides and tapasin within the endoplasmic reticulum, together with the dynamics of egressed peptide-MHC complexes at the cell surface.
A number of simplifying assumptions were made when constructing the model: (i) Peptides P i are supplied to the ER at rate g i and then degraded or removed from the ER at rate d P .
Since different peptides P i can have different levels of abundance within the cytoplasm and different rates of TAP transport, each peptide is associated with its own generation rate g i . (ii) MHC class I heavy chain and b 2 m are assumed to represent a single unit, where b 2 m dissociation from empty class I heavy chain is interpreted as a form of MHC degradation. (iii) All peptides are assumed to have a similar rate of binding b to MHC, such that peptide affinity is determined by a peptide-specific rate of unbinding u i , and is defined as 1=u i . This is motivated by the measurements in [11]. (iv) Since MHC, tapasin and peptide continually cycle between the ER and Golgi apparatus [12,13], we do not explicitly represent the Golgi as a separate compartment. Instead, we consider our ER compartment to include both the ER and Golgi, where the rate of egress e represents the rate of transit from the Golgi to the cell surface. By representing this process as a firstorder reaction, we are making the simplifying assumption that the quantity of peptide-MHC complexes which egress is related to the quantity of complexes in the ER. (v) MHC can load peptides in the presence of tapasin at a higher rate c, which implicitly models the stabilization of TAP molecules by tapasin, but we neglect egress of tapasin-bound MHC, since tapasin retains MHC by bridging it to the TAP transporter [14]. (vi) Tapasin can increase the rate of peptide unbinding from MHC by a factor q [8], while peptide can increase the rate of tapasin unbinding from MHC by a factor v [15]. Tapasin has been shown to increase the peptide off-rate to a similar extent for peptides with a range of off-rates, though some variation has been shown for certain classes of peptide [8]. (vii) We neglect egress of empty MHC, which is retained and recycled in the ER by the chaperone calreticulin [16,17]. (viii) Furthermore, we assume that b 2 m dissociation from peptide-loaded or tapasin-bound class I heavy chain is negligible compared to b 2 m dissociation from empty class I heavy chain. (ix) Once at the cell-surface, peptide unbinds from MHC irreversibly at rate u i , and empty MHC is degraded at rate d Me . These assumptions can be refined in future iterations of the model.
Predicting peptide optimization over time MHC class I HLA-B alleles were previously shown to differ in their ability to optimize their peptide cargo over time, both in the presence and absence of tapasin [6]. Specifically, the HLA-B4402 (B4402) allele was shown to be highly dependent on tapasin for peptide optimization, while the HLA-B2705 (B2705) and HLA-B4405 (B4405) alleles were shown to be less tapasin-dependent. B2705 is of particular interest because it is a susceptibility factor for certain autoimmune diseases and is associated with long-term non-progression of HIV [18]. Therefore, we sought to use our peptide optimization model to explain the variation in tapasindependence between HLA-B alleles, through a combination of model simulation and Information Theory.
We simulated pulse-chase experiments [6] using the peptide optimization model of Fig. 2 , with representative peptides of low, medium and high affinity (Text S1). The experiments followed the thermostability of fixed cohorts of MHC class I complexes over time, making use of the known correlation between the thermostability of complexes and the affinity of their peptide cargo. Specifically, complexes stable at 50 0 C were shown to contain only high affinity peptides, complexes stable at 37 0 C were shown to contain a combination of medium and high affinity peptides, while all complexes were shown to be stable at 4 0 C, including empty MHC. Since the measurements correspond to both ER-localized and egressed peptide-MHC complexes, our assessment of the model was performed by comparing total peptide-MHC complexes with the 4 0 C measurement, total medium and high affinity complexes with the 37 0 C measurements, and total high affinity complexes with the 50 0 C measurements. Since many of the kinetic parameters of the model have not previously been measured directly, due to the technical difficulties involved in obtaining such measurements, we used heuristic search methods to infer the parameter values from the experimental data [6] (see Methods). Essentially, this involved finding values for the parameters which minimized the deviation between the experimental data and the corresponding model simulation. Using this approach, we investigated how allelic variation in HLA-B might affect peptide optimization, by distinguishing between allele parameters, which were allowed to vary between alleles, and fixed parameters, which were assumed to be invariant between alleles. Each hypothesized set of allele parameters defined a variant of the model, which possessed a different intrinsic ability to reproduce the observed dynamics.
The Bayesian Information Criterion (BIC) [19] was used to quantify the performance of each set of allele parameters (equation (18) in Methods). BIC incorporates a term which penalizes the deviation of the simulation from the data, and a second term which penalizes increasing numbers of allele parameters. Therefore, BIC can be used to assess a range of models by taking into account the added cost of additional unconstrained variables. Since the dynamics of peptide optimization varied considerably between HLA-B alleles in the absence of tapasin [6], we reasoned that at least one allele parameter must be tapasin-independent. To incorporate this insight whilst focusing on the principal contributors to allelic variation, we examined combinations of up to two allele parameters, with at least one tapasin-independent parameter selected from b, d M , d Me , e, g i and u i .The best BIC scores were obtained when the peptide on-rate b was the only allele parameter (470.38), and when both b and the rate of egress e were the allele parameters (469.39; Fig. 3), suggesting that at least peptide on-rate is allelespecific. However, having both b and e as allele parameters required unrealistically fast egress of B2705 and B4405 complexes to obtain a closer fit to the data ( Fig. S2 in Text S1). Therefore a single allele parameter ( Fig. 4) b was used, which was able to effectively account for the experimental data [6]. Specifically, in the absence of tapasin B4402 exhibited worse time-dependent optimization than both B2705 and B4405 ( Fig. 4 A), while in the presence of tapasin B4402 exhibited better time-dependent optimization than both of these alleles (Fig. 4 B). To ensure that the MCMC search algorithm was robust to random variations, and could reproducibly generate consistent parameter estimates, we produced 10 different chains for each model hypothesis. For the allele-specific b model, 8 out of 10 chains converged to BIC values between 470.38 and 470.69, while the other two chains performed poorly. We next plotted the mean and standard deviation of the posterior distributions of the model parameters for each of the 10 chains, which revealed that the 8 high performing chains had overlapping posterior distributions ( Fig.  S3 in Text S1), and were therefore producing consistent parameter estimates.
To understand the effects of an allele-specific peptide on-rate on peptide optimization, we plotted MHC complexes with high, medium, and low affinity peptide separately, and distinguished free and tapasin-bound MHC complexes within the ER and at the cell surface ( Fig. S4 in Text S1). For B4402 without tapasin, an intrinsically low peptide on-rate meant that the majority of B4402 complexes remained in the ER without peptide, resulting in very low optimization. For B2705 and B4405 without tapasin, an intrinsically high peptide on-rate meant that these alleles rapidly bound their peptide cargo and exhibited good time-dependent optimization. In contrast, for B4402 with tapasin, most complexes first bound to tapasin and were subsequently able to rapidly bind peptides and optimize their peptide cargo, presenting almost exclusively high affinity peptides at the cell surface. For B2705 and B4405 with tapasin, the intrinsically high peptide on-rate meant that peptide out-competed tapasin for binding to MHC, such that a higher proportion of peptides followed the non-tapasin pathway, Figure 3. Selection of HLA-B allele parameters. The horizontal axis indicates the set of parameters that were allowed to vary between alleles. The vertical axis quantifies the Bayesian information criterion (BIC) of the best parameter values for a given set of allele parameters. BIC penalizes deviations of the model simulation from the experimental data, whilst also penalizing models with more variable parameters, implying that low BIC values correspond to more representative models. The best parameter values for a given set of allele parameters were inferred using the Filzbach MCMC software (see Methods). doi:10.1371/journal.pcbi.1002144.g003 resulting in reduced optimization. Thus, variation in the intrinsic ability of free HLA-B alleles to bind peptide in the absence of tapasin was shown to be the most likely explanation for allelic variation in peptide optimization, both in the presence and absence of tapasin. For all three alleles, cell surface optimization could not be improved by modifying most other parameters in the model (c, u T , v, q and e) (Fig. S5 in Text S1). This indicates that the balance between peptide binding b and tapasin binding b T is a major determinant of peptide optimization, achieved by controlling the effectiveness of the tapasin-mediated pathway. The prediction that allele-specific tapasin dependency results from variations in peptide binding to MHC class I molecules is consistent with analysis from molecular dynamics simulations, which suggest that tapasin stabilizes peptide-receptive conformations [20]. This stabilization in the presence of tapasin is represented in our model by setting the binding rate c to be allele-independent and greater than or equal to the binding rate b. In the absence of tapasin, MHC class I molecules of different alleles may have varying levels of peptide receptiveness, which is represented in our model by allowing b to vary between alleles. Corresponding experimental results [6] are also reported (circles). Simulations were conducted for representative low, medium and high affinity peptides with a separate dissociation rate u i and generation rate g i for each peptide P i , and a separate peptide binding rate b j for each HLA-B allele M j (Table S1 in Text S1; Protocol S1). doi:10.1371/journal.pcbi.1002144.g004 Kinetic control of peptide optimization Having established a hypothesis which explains how MHC alleles experience differential tapasin-dependence, we sought to identify the mechanisms that determine the extent and rate of peptide optimization, both in the presence and absence of tapasin. Peptide optimization is the process by which high affinity peptides are selected for presentation at the cell surface [6]. Peptide-MHC complexes generally need to be stable for many hours or days at the cell surface in order to effectively elicit an immune response [21], yet peptide optimization in the ER is typically limited to tens of minutes [3,6,22]. This requires optimization beyond the limit that would be obtained in equilibrium. How such high optimization is achieved in so little time is still not well-understood [3].
One way to increase the extent of peptide optimization is for peptide-MHC complexes to be retained in the ER for an extended period prior to egress, so that unstable peptides have an opportunity to unbind [22]. However, delaying egress also increases the time for complexes to reach the cell surface. Therefore, a trade-off exists between the extent of optimization and the rate at which this optimization can be achieved. We quantify this trade-off by calculating the relative probabilities of MHC egress and peptide unbinding.
Consider an MHC complex containing a peptide with off-rate u i (Fig. 5 A). The complex can either egress to the cell surface at rate e, or the peptide can unbind at rate u i . The probability of each event is proportional to its rate, such that the probability of egress is given by e= u i ze ð Þ. The competition between unbinding u i and egress e defines a peptide filtering step, where the basic filtering mechanism is comparable to principles of kinetic proofreading [23]. Let MeP i ½ Ã denote the expected number of peptide-MHC complexes that egress to the cell surface before the peptide can escape (Fig. 5 A). If there are N i MHC complexes containing peptides with off-rate u i in the ER, we expect N i e= u i ze ð Þ to egress and the remainder to release their peptide cargo. For very high e, all N i complexes will egress irrespective of their peptide cargo. For very low e, the number of egressed complexes will tend to N i e=u i . Let O MeP i ð Þdenote the proportion of egressed MHC complexes containing peptides with off-rate u i (Fig. 5 A). This defines a measure of peptide optimization. For very high egress we observe no optimization, where the proportion of peptides at the cell surface is equal to the proportion of peptides in the ER. For very low egress we observe maximum optimization, where the proportion of peptides at the cell surface varies inversely with the peptide off-rate.
The introduction of tapasin provides an additional filtering step (Fig. 5 B), involving a competition between peptide unbinding u i q and tapasin unbinding u T v. Let ½MeP i Ã denote the expected number of MHC complexes that unbind from tapasin and egress to the cell surface before the peptide can escape (Fig. 5 B). For very high e and x~u T v=q, all N i complexes will egress irrespective of their peptide cargo. For very low e and x, the number of egressed peptide-MHC complexes will tend to N i ex=u 2 i . The number of egressed complexes therefore varies with 1=u 2 i in the presence of tapasin (Fig. 5 B), compared with 1=u i in the absence of tapasin (Fig. 5 A). This implies that tapasin enhances presentation according to peptide affinity 1=u i , in agreement with experimental results [5]. Since the proportion of egressed complexes now varies inversely with the square of peptide off-rate, low affinity peptides are much more likely to escape than high affinity peptides, resulting in improved peptide optimization.
The peptide filtering relation presented above also holds for the full dynamical systems model of Fig. 2 , in which peptides can bind and unbind multiple times to MHC. By translating the reactions of Fig. 2 to a set of ordinary differential equations, we obtained the following expression for the steady-state concentration of peptide-MHC complexes at the cell surface (see Methods): where x~u T v=q. The equation includes the ER peptide filtering steps described in Fig. 5, together with peptide optimization at the cell surface given by 1=u i , where peptides with a lower off-rate u i are more likely to remain bound to MHC. The equation also quantifies the ratio of egressed peptide-MHC complexes that are loaded in the presence and absence of tapasin, given by the ratio of ½TM Ã cx=(u i zx) to b½M Ã . Assuming peptide loading takes place via the tapasin pathway (½TM Ã &½M Ã ) and that peptides have a high turnover in the ER [24], characterized by high generation and degradation rates (½P i Ã &g i =d P ), we can simplify equation (1) as where C~c½TM Ã =d P . This corresponds to an upper bound on peptide optimization in the presence of tapasin. In the absence of tapasin, the equation for ½MeP i Ã is the same as (2) but without the tapasin optimization step x=(u i zx). This implies that tapasin enhances peptide presentation according to peptide affinity 1=u i , in agreement with the analysis of Fig. 5 and experimental findings [5].
To further place our insights in a biological context, we used the dynamical systems model to identify the mechanisms that determine the rate of peptide optimization. Consider the filtering step between peptide unbinding u i q and tapasin unbinding u T v. Tapasin can enhance peptide optimization to the same extent either by increasing the peptide off-rate by a given factor q, or by decreasing the tapasin unbinding rate by the same factor. However, decreasing the tapasin unbinding rate essentially delays the transit of MHC to the cell surface, resulting in slower optimization. In contrast, increasing the peptide off-rate allows tapasin to increase the extent of peptide optimization while still maintaining a rapid flux of peptide-MHC complexes to the cell surface.

Predicting peptide optimization at steady-state
To further probe the applicability of our model, we investigated whether it could be used to predict peptide optimization at steady state. Previously, the effects of tapasin on steady-state peptide optimization were measured for peptides in the MHC class I allele H22K b (K b ) [5]. The experiments were conducted using four target peptides, obtained by performing substitutions at positions 5 and 8 of the amino acid sequence SIINFEKL. Peptide off-rates were measured in RMA-S cells and each of the target peptides were introduced as minigenes into a tapasin-deficient cell line (.220) and into the same cell line transfected with tapasin (.220.Tpn). Steady-state levels of cell surface peptide-MHC complexes were measured by flow cytometry using mAb 25.D1 (Fig. 6 C, Table S3 in Text S1), which specifically recognizes the SIINFEKL peptide variants bound to K b [25]. Total cell-surface MHC was also measured with mAb Y3, which recognizes empty and peptide-occupied K b (Fig. 6 D, Table S3 in Text S1).
We complemented previous experimentation [5] by measuring the off-rates of the target peptides in .220 cells directly.  Fig. 6 A). We also used Y3 to measure total MHC during BFA incubation of cells with no target peptide, in the presence and absence of tapasin, to provide an indication of the off-rates of endogenous peptides presented by .220 and .220.Tpn cells (Fig. 6 B).
We simulated the above experiments using the peptide optimization model of Fig. 2 with a parameter set specific to H22K b (Table S2 and Fig. S6 in Text S1). The full range of endogenous peptides was characterized by two representative peptides with off-rates u 1 and u 2 , and supply rates g 1 and g 2 (Text S1). Each experiment was simulated using one of the target peptides with off-rate u i , together with the representative endogenous peptides. Since target peptides were expressed at approximately equal levels inside cells [5], we assumed that they were generated at the same rate g i~g0 . The model simulations agreed with the trends observed experimentally, accurately recapitulating the enhancement of steady state optimization conferred by tapasin (Fig. 6 C, D, Table S3 in Text S1). However, we observed that the model did not fit the experimental data for SIINFEKM as well as for the other target peptides. We hypothesized that the poor fit could be caused by increased TAP transport of the SIINFEKM peptide, due to a change in the terminal residue at position 8 [26]. To explore this idea, we increased the generation rate of SIINFEKM by a factor of 2.5, which gave a better fit to the experimental results (Fig. 6 D). This  (Table S2 in Text S1; Protocol S2). Each simulation computes the steady state of the model with three types of peptide: two background peptides P 1 and P 2 and one of the four SIINFEKL peptide variants P i (u i estimated from data in panel A).  [5]. (D) Total steady-state peptide-MHC complexes (cell surface), comparing simulation with measurements of Y3 from [5]. Simulated values were scaled by a proportionality factor for optimally overlapping the 25.D1 data (with SIINFEKM removed) and the Y3 data (all points) (Text S1). (B-D) The x-axis shows the relative affinity of peptides given by the inverse of the off-rate. Steady state concentrations were obtained by equating the right hand sides of the ODEs to zero. Steady state concentrations in tapasin-deficient cells were simulated by setting g T~0 . (E-G) For quantifying egression of peptide-MHC complexes, .220.K b (E) and .220.K b .Tpn (F) were pulsed for 10 min with 35 S-Met/Cys and chased for the indicated times (min). Y3 immunoprecipitates were digested with endoglycosidase-H (EndoH) and SDS-PAGE hypothesis further highlights the potential importance of peptide supply in predicting relative presentation levels [27], as can be seen from the peptide filtering relation (2). Although experimental measurements were only obtained for four distinct peptides, the model predicts the presentation levels for a continuum of peptide off-rates over a broad range, which can be checked in future experiments.
To distinguish between optimization resulting from delayed tapasin unbinding versus enhanced peptide off-rate, we measured the time taken for a fixed cohort of pulse-labeled MHC complexes to reach the cell surface by measuring endoglycosidase-H (EndoH) resistance (Fig. 6 E-G). By taking into account the temporal constraints of the EndoH data, we found that enhanced peptide off-rates were required to allow increased peptide optimization in the presence of tapasin without significantly delaying the transit of peptide-MHC complexes to the cell surface (see Fig. S7 in Text S1). Further parameter variation analysis indicated that cell surface optimization is nearly maximal in the H22K b model with respect to q, u T , v and c, but could be improved by reducing e (Fig.  S8 in Text S1). However, reducing e decreases the export of peptide-MHC complexes, suggesting a possible trade-off between optimization and the efflux of new information concerning cellular protein content.

Predicting optimization of viral peptides
To illustrate how the peptide optimization model may be used in more realistic scenarios, we simulated the presentation of HIVderived peptides using our models for the HIV-associated allele HLA-B2705 (B2705), and HLA-B4402 (B4402) for comparison. Peptides between 8 and 10 amino acids in length were identified from the Gag-Pol polyprotein (UniProt; accession P03367) and assessed for their off-rates using the BIMAS prediction algorithm [28]. For B2705, the slowest off-rate identified was for the immunodominant KRWIILGLNK (positions 262-272) epitope [29] (off-rate: 3:8508|10 {7 s {1 ). For B4402, the allele parameters of the BIMAS algorithm were not available, so we quantified off-rates based on the BIMAS algorithm parameters for the closely related allele B4403. The off-rates identified were generally higher than for B2705 (Fig. 7 A). When comparing simulations of B4402 in the presence and absence of tapasin, the highest affinity peptide AETGQETAY (positions 1250-1258; off-rate: 2:1393|10 {5 s {1 ) was enhanced by a factor of 445 by tapasin (Fig. 7 B), though cell surface presentation was over 25 times less than the presentation of KRWIILGLNK by B2705 (Fig. 7 B). Despite B2705 being a relatively tapasin-independent allele [6] (Fig. 4), tapasin significantly enhanced presentation of peptide KRWIILGLNK by a factor of 120 (Fig. 7 C). and autoradiography were performed. Arrows indicate K b heavy chain resistant (R) and sensitive (S) to EndoH digestion. EndoH analysis of H22K b was performed as described previously [6]. (B-D, G) The solid lines indicate model simulations and triangles indicate measured data-points. The experimental data for (A,B,E-F) is novel, while the experimental data for (C,D) is from [5]. doi:10.1371/journal.pcbi.1002144.g006 Figure 7. Simulation of cell surface presentation of HIV virus peptides by HLA-B2705. The sequence of the HIV-1 polyprotein Gag-Pol was obtained from the UniProt online resource (accession P03367). All peptides between 8 and 10 amino acids in length were then derived from the sequence and assessed for their off-rates using the BIMAS prediction algorithm (http://www-bimas.cit.nih.gov/molbio/hla_bind [28]). The peptides were then simulated by assuming that they are all supplied into the ER via TAP at an equal rate, such that the total supply rate is equal to the total supply rate estimated for Fig. 4. As the algorithm predicted many peptides to have the same off-rate, peptides were clustered for ease of computation. (A) The number of peptides with a given peptide off-rate, as calculated by BIMAS. (B) Steady-state cell surface presentation of peptide-MHC complexes as a function of peptide off-rate. Peptide supply was assumed to be constant for each individual peptide. Therefore, the supply rate associated with a particular off-rate is simply scaled by the number of peptides with that off-rate, as quantified in A. The lowest off-rate (highest affinity) peptides for B2705 (KRWIILGLNK) and B4403 (AETQCETAY) are indicated. Simulations were performed for the presence and absence of tapasin, as indicated in the key. (C) Enhancement of cell surface presentation by tapasin was computed by dividing simulated tapasin-sufficient presentation by simulated tapasin-deficient presentation for each peptide. The results of the HIV simulations illustrate the extent to which tapasin can affect a downstream immune response. Theoretically, tapasin can enhance presentation by up to a factor 1=u, where u is the off-rate of the peptide from MHC (Fig. 5). However, the characteristics of the MHC allele, such as the allele-specific peptide on-rate, can significantly alter the effect of tapasin on the presentation of a given peptide. Our model allows differences in presentation levels to be quantified by taking into account peptide supply and peptide off-rate, together with the effects of tapasin and the binding properties of the MHC class I allele under consideration. In particular, our analysis of the HIV-1 Gag-Pol polyprotein provides a specific quantitative prediction for the cell surface presentation of the immunodominant KRWIILGLNK by HLA-B2705. By simulating the range of peptides derived from Gag-Pol, representing a range of off-rates, we observe that the enhancement by tapasin is independent of peptide supply, instead being wholly determined by the peptide off-and on-rates (Fig. 7 C). doi:10.1371/journal.pcbi.1002144.g007

Discussion
The optimization of peptide-MHC class I complexes at the surface of antigen presenting cells is one of the key factors that determines the hierarchy of the T-cell response to a complex antigen [1]. Peptide optimization is also important for vaccine design, where vaccine peptides compete with endogenous peptides for presentation [1,30]. In this paper we propose a dynamical systems model of MHC class I peptide optimization, which takes into account the supply of peptides in the cytosol, the affinity of peptides to MHC and the interactions between peptide and MHC at the different stages of the optimization process, both within the ER and at the cell surface. The model also incorporates the effects of tapasin, which is known to increase peptide optimization [5] and to affect different MHC class I alleles to different extents [6]. This variation in tapasin dependence may protect from viral immune evasion strategies such as tapasin inhibition by an adenovirus [31].
The dynamical systems model is firmly grounded in experimental data, and techniques already exist to measure many of the model parameters [5][6][7][8]. The model therefore allows a multitude of experimental results to be unified within a common framework, so that a range of mechanistic hypotheses can be formulated and tested. We derive a peptide filtering relation which, for the first time, provides a mechanistic explanation for experimental data on MHC class I peptide optimization, both over time [6] and at steady state [5]. Specifically, it suggests that tapasin enhances peptide off-rate in order to improve peptide optimization without significantly delaying the transit of MHC to the cell surface.
We have also shown that an allele-specific peptide on-rate is the most likely mechanistic explanation for differences in peptide optimization across HLA-B alleles. A possible interpretation is that differences in peptide on-rate are due to allelic differences in molecular conformation. For example, alleles such as B4402 could adopt a closed conformation, reducing the ability of peptides to bind MHC, while alleles such as B2705 could adopt a more open conformation, allowing peptides to readily bind MHC, as suggested in [32]. When tapasin binds to MHC the peptide binding groove may then adopt a peptide receptive conformation, allowing MHC to bind peptides more readily, as suggested in [20]. Although allelic differences in the conformation of MHC class I are largely peptide-independent, variations in the on-rates of different peptides have nevertheless been observed. These variations can be incorporated in future versions of the model by allowing a separate on-rate for each peptide. However, published estimates indicate that variations in the affinity of peptide-MHC interactions are mostly governed by variations in peptide off-rate [11], supporting our assumption that the on-rate is allele-specific and largely peptide-independent. Although the current model makes a number of simplifying assumptions on the antigen presentation process, the model can be readily extended to incorporate additional details as more data are acquired. These details could include the explicit contribution of TAP transport, proteasomal cleavage and cytosolic protein abundance to ER peptide supply [26,33,34]. At present these mechanisms are only implicitly represented in the model via peptide-specific supply rates g i . Further extensions could also include conformational changes in MHC [7], and chaperones such as ERp57 and calreticulin which are known to influence total cell-surface presentation [15,17]. Since the mechanisms by which additional chaperones interact with MHC class I are only partially known, we can investigate a variety of hypotheses by using our Information Theoretic framework to assess allele-specific chaperone-dependency. In the future, coupling model analysis with additional experimental measurements will enable quantitative predictions of peptide optimization for a wide range of MHC class I genotypes. Having a robust model, known to make accurate predictions, will improve our ability to assess the efficacy of vaccines involving multiple peptides, and will provide a quantitative means to prioritize different vaccination strategies.
The current work is part of a broader research programme to use experimental data to build credible mathematical models of immunological processes, ranging from relatively simple examples to complex systems such as organ-specific autoimmunity. The resulting models can then be used to make specific and testable predictions that relate directly to immunological function. Subsequent iterations offer an opportunity to refine or develop the models from the simple to the complex, or from the static to the time-resolved, at the molecular, cellular or organ level.

Ordinary differential equation representation of the dynamical systems model
The chemical reaction representation of the dynamical systems model of Fig. 2 is as follows:

Me
By assuming mass-action kinetics, we converted the system of reactions to a set of ordinary differential equations (ODEs), where ½X denotes the concentration of a given species X and ½X 0 denotes rate of change in concentration.

Steady-state analysis of the ODE representation
Equation (1) was derived by considering the steady state (equilibrium) solutions of the ODE representation. By equating ½MP i 0 , ½TMP i 0 and ½MeP i 0 with zero, we obtained the following expressions for ½TMP i Ã and ½MP i Ã : By substituting (12) in (11) we obtained the expression for the steady-state concentration of egressed MHC complexes with a given peptide P i in equation (1) ½MeP where x~u T v=q. The equation incorporates the peptide filtering steps described in the main text, where e=(u i ze) denotes the optimization of free MHC in the ER, while x=(u i zx) denotes the optimization of MHC bound to tapasin in the ER. For a given peptide with off-rate u i , the optimization of free and tapasin-bound MHC is therefore fully determined by e and x, respectively. Moreover, MHC complexes optimized in the presence of tapasin will always be subject to an additional optimization step after tapasin unbinding. The equation also incorporates optimization of MHC at the cell surface given by 1=u i , where peptides with a lower off-rate u i are more likely to remain bound to MHC.

Optimizing model parameters with respect to experimental data
Heuristic search methods were used to fine-tune model parameters, based on the available experimental data. Our approach was to minimize a cost function, defined as the sum of the squared differences between experimental data y (i) d (i~1,n) and corresponding model simulated output y (i) m , subject to an arbitrary proportionality constant a. i.e.
where P is the space of search parameters, which may be the full parameter set or a subset thereof. For the inner minimization problem (14), it is possible to assign an optimal a~a Ã (=0), by equating the partial derivative of D (with respect to a) to 0.
Approximate minimizers of the multi-dimensional objective function (13) were found using a Metropolis-Hastings (MH) algorithm.
During execution of the MH algorithm, a Markov chain of proposal parameter sets is formed. Starting from an initial parameter set p~p 0 with associated objective function value D(p 0 )~D 0 , the algorithm iteratively searches neighboring parameter sets by accepting or rejecting new proposal parameter sets at each step. Neighboring pointsp p are proposed with probability Q(p,p p), according to a jump rulê The chain moves to the new pointp p according to an acceptance criterion, which makes a probabilistic choice about whether to acceptp p. Given an observation u drawn from U(0,1), the proposal point is accepted providing uv min Q(p,p p)P(p p) Q(p p,p)P(p) ,1 where P(p) is the probability that the parameter set p matches the data and Q(p,p p) is the proposal density (which we fix to be symmetric). Assuming the deviations from the experimental data are Gaussian distributed and that D(p) makes only small jumps, the acceptance ratio is approximately given by P(p p) P(p) &e (D(p){D(p p)) Note that when D(p)wD(p p) we always move top p, because the exponential function of a positive number is greater than one. The algorithm iterates until some termination condition is reached, such as a maximum number of iterations or a convergence condition. All parameter searches were performed using 'Filzbach', a software library for carrying out Metropolis-Hastings Markov chain Monte Carlo parameter estimation in C++ or C#.
Filzbach is under development in the computational science lab at Microsoft Research Cambridge and is available for download, complete with a suite of example uses, via http://research. microsoft.com/science.

Model selection using the Bayesian Information Criterion
To assess which model parameter(s) should be allele-specific, different hypotheses were compared using the Bayesian Information