BACKGROUND
class ii mhc is an essential protein complex in the adaptive immune system that is involved in activating helper t cells.
it is found on the surfaces of specialized immune system cells, such as dendritic cells and macrophages, where it binds fragments of extracellular peptides and presents them to cd4+ helper t cells.
class ii mhc is highly polymorphic and each allotype generally has different peptide binding preferences.
the computational prediction of such binding affinities can be applied to the important problem of finding promiscuous epitopes that bind to multiple allotypes for use in rational vaccine design.
also these predictions may aid in understanding autoimmune diseases, many of which have been linked to particular class ii mhc alleles  <cit> .
characteristic properties of peptide binding to class ii mhc molecules can be inferred from available high-resolution x-ray structures of peptide-mhc complexes.
the binding cleft in class ii mhc is open at both ends so that it generally binds longer peptides, 15- <dig> residues in length, with the peptide n- and c-terminii extending out of the cleft.
this is in contrast to class i mhc, which has a peptide binding cleft that is closed at both ends so that it only accommodates shorter peptides of 8- <dig> residues.
an examination of x-ray crystal structures of peptide bound class ii mhc molecules reveals that the peptides bind in an extended polyproline ii helix with the conserved backbone structure maintained by conserved hydrogen bonds with the mhc  <cit> .
peptide side chains within a contacting nine residue core segment of the peptide bind into pockets on the mhc surface  <cit> .
many of the polymorphic residues occur in these pockets and so determine the distinctive peptide binding preferences of each mhc allotype.
although the core residues of the bound peptide assume similar backbone conformations, each peptide can generally bind in multiple registers.
this makes the prediction of peptide binding affinities for class ii mhc particularly challenging since, unlike class i mhc, all possible binding registers must be considered.
prediction methods for class ii epitopes remains an active area of research, motivated both by their biomedical importance and the difficulty of achieving high accuracy.
early prediction methods fit the total peptide binding energy  <cit> , binding motif  <cit> , geometric average binding affinity  <cit> , or sequence alignment profile  <cit>  in a particular register to a linear combination of contributions from individual residues, and represented them as binding profile matrices.
the scores for all possible peptide binding registers were calculated and either the maximum value or sum were used as a total peptide binding score.
later methods employed various machine learning and data fitting approaches to prediction including partial least squares   <cit> , gibbs sampling  <cit> , linear programming  <cit> , support vector machines   <cit> , and kernel methods  <cit> .
one prediction method, called smm-align, successfully combined two methods, binding profile matrices and gibbs sampling  <cit> .
we introduce a new model for predicting binding affinities of peptides to class ii mhc called the regularized thermodynamic average, or rta, model.
three features are incorporated into the model:  independent binding of individual peptide residues so that the total binding energy is a sum of the binding energies of each core peptide residue,  a thermodynamic average over all possible peptide binding registers, and  regularization of the model to prevent overfitting.
feature # <dig>  the approximate independence of peptide binding on peptide residue interactions, is supported experimentally  <cit>  and is a basic assumption in most prediction models to date.
indeed, it is the primary basis for binding profile methods.
the importance of contributions of suboptimal peptide binding registers to the total binding affinity motivates feature # <dig>  experimental studies have found such cases in which a peptide binds with significant affinity to class ii mhc in multiple registers  <cit> .
in fact, such peptide binding in multiple registers may be involved in the origin of autoimmune diseases  <cit> .
the thermodynamic average in the rta model insures that contributions from suboptimal binding registers, which may be significant, also contribute to the total binding affinity.
regularization  is especially important in prediction models that contain a relatively large number of parameters, such as ours.
regularization is implemented in rta using an l <dig> penalty function, like in lasso regression  <cit> .
this has the advantage over alternative penalty functions, such as the l <dig> penalty employed in smm-align, that many coefficients vanish at lower cutoff values so that the number of model parameters is reduced.
in addition, unlike in many previous prediction methods, the p <dig> residue for peptides binding to hla-dr is not restricted to be hydrophobic in the rta model.
finally, we use deterministic local optimization methods to fit the model parameters, which we found to converge to better solutions than commonly employed stochastic optimization methods.
the rta prediction method was tested using cross-validation sets from an smm-align study  <cit> , for comparison with results from  <dig> different prediction methods examined in that study, and also using larger cross-validation sets created from binding data from the iedb database  <cit> .
the rta method was also evaluated using additional peptide binding data sets from lin et al.
<dig>  <cit>  and el-manzalawy et al.
<dig>  <cit> .
a novel method to efficiently create data sets with minimal sequence overlap was developed to create the new cross-validation sets.
the performance of rta in predicting peptide binding affinities and in predicting the correct peptide core binding residues were tested.
finally, we examined the relative contribution of suboptimal peptide binding registers to the overall binding affinity.
methods
peptide-mhc binding data set
peptide binding data for the  <dig> mhc allotypes in table  <dig> were downloaded from the iedb database http://www.immuneepitope.org <cit> .
all ic <dig> values from either radioactively or fluorescently labeled peptide competition binding assays were collected.
the median value was used for the few cases in which multiple experimental values were measured.
the total numbers of data for each allotype are given in table  <dig>  these datasets were then divided into minimally overlapping 5-fold cross-validation sets using the method described in "creating independent cross-validation sets".
prediction model
as mentioned above, the prediction of peptide - class ii mhc binding affinity involves not only estimating the affinity of the bound core portion of the peptide but also accounting for multiple binding registers.
the model begins with the approximation that the contributions to the total binding free energy of binding individual peptide residues to each mhc pocket are approximately independent.
this assumption has been tested experimentally  <cit>  and is the basis of matrix-based prediction methods in which the total peptide-mhc binding energy is a sum of pocket specific residue binding energies  <cit> .
furthermore, x-ray structures of peptide - class ii mhc complexes show that the 9-residue core bound portions of the peptides adopt similar backbone conformations so that the total score over an ungapped 9-residue segment of the peptide sequence is appropriate.
the peptide-mhc binding data from competition experiments are given as ic <dig> values.
ic <dig> is not a direct measure of binding affinity but it can be related to the inhibition constant, ki, through the cheng-prusoff equation
is the concentration of the free labeled peptide and kd is its binding affinity to the mhc molecule  <cit> .
it is assumed that the labeled peptide concentration is sufficiently low, i.e.
= kd, so that ki ≈ ic <dig>  this is an implicit assumption in other prediction methods as well.
according to this approximation, experimental binding affinities were calculated as Δgexp = - kt log, in which kt is the boltzmann constant times the absolute temperature, which is approximately  <dig>  kcal/mol at room temperature.
each peptide can bind the mhc molecule in different registers defined by which 9-residue core peptide segment contacts the mhc.
the length of peptide k is denoted by l. the register will be denoted by the index m, which varies from  <dig> for the first  <dig> residues as the core segment to l -  <dig> for the last  <dig> residues as the core segment.
the sequence of peptide k is represented by matrix elements , which are equal to  <dig> if residue type j appears at position i in the peptide sequence and are equal to  <dig> otherwise.
the affinity of peptide k binding to the mhc molecule in register m is then the linear function
in which the coefficients βij are to be fit using the experimental data.
the apparent total binding affinity for peptide k can be calculated as
in which  is the fraction of bound peptides binding in register m and  is the affinity for binding in register m, calculated in the model using eq.
<dig>  an entropic contribution to the free energy in eq.
<dig> was omitted because it did not improve prediction performance .
the fraction of peptide bound in a particular register is related to the binding affinities in eq.
<dig> by
, , and  are the concentrations of unbound peptide, unbound mhc and peptide bound to the mhc in register m, respectively.
is the dissociation constant for peptide k binding to the mhc in register m. by combining eqs.
<dig> and  <dig> one arrives at an expression for the total binding affinity as an exponentially-weighted average of the affinities for the peptide binding in different registers
many prediction methods for peptide-mhc binding make the further approximation that the largest binding affinity dominates the sum in eq.
<dig> so that only the largest binding affinity over all possible registers is included in the prediction model.
because suboptimal terms are comparable in magnitude to the optimal term in many cases, we will retain all terms in eq.
<dig>
the final component in the prediction model is regularization.
this is needed to avoid overfitting, since the model contains  <dig> parameters, which is relatively large compared to the number of binding data for each mhc type.
regularization is accomplished by imposing an upper bound on the l <dig> norm of the coefficient amplitudes
with the cutoff parameter t. this type of constraint is also used in lasso regression  <cit>  and has the desirable property that as the cutoff t is lowered, an increasing number of coefficients become zero so that the number of nonzero parameters in the model is reduced.
by combining these components, the binding affinity for peptide k, Δg, is calculated as
the parameters in this model, {βij, i =  <dig>  k,  <dig>  j =  <dig>  k, 20}, are fit to the experimental binding affinities, , by minimizing the mean square error
subject to the constraint in eq.
<dig>
optimization method
to solve the optimization problem minimizing eq.
<dig> under the constraint in eq.
<dig> one first splits the variable βij into the difference of two nonnegative variables  and , so that βij =  - .
the constraint eq.
<dig> then has as the left side the sum over both new variables.
overall one has a linearly constrained nonlinear and nonconvex optimization problem.
in principle, global optimization methods would have to be applied to find the global minimizer.
due to the dimensions of the problems these would be non-deterministic methods such as some of the many metaheuristics, which include simulated annealing and genetic algorithms.
these methods require a large number of evaluations and due to their stochastic character would have to be run several times in order to increase the likelihood of finding the global optimum, although a guarantee for that is impossible.
after initial tests we decided to instead use local solvers.
these are very efficient and through various measures the chances of getting very good local minima can be increased substantially.
using local solvers with "multistart" or several often randomly generated starting guesses is another way of solving global optimization problems.
it is also implemented in several software packages.
we instead used simple constant starting guesses.
all of our results were obtained using the constant starting values of  <dig> or  <dig> for all components of β+ and β-.
we did, however, use a sequence of different t values in the regularization constraint eq.
<dig> and then chose the best  minimum found for the solvers used.
it is significant that the best solutions were obtained with values of the regularization parameter t that were small enough to restrict the fit βij values.
this shows that the additional l <dig> constraint in eq.
<dig> helps alleviate overfitting.
in order to be able to easily call a variety of solvers we phrased the problem in the modeling language ampl  <cit> .
in order not to have to list many separate citations we state that we used the applicable  solvers installed at neos  but run locally, not through this free service in which we  are also heavily involved.
in the way described above we generated the values in tables  <dig> and  <dig>  for table  <dig> it was possible to just run the solver that had best solved the corresponding case in table  <dig> on the data set and evaluate for the test sets.
we did not use the optimal coefficients found but applied our method to the full training set in each case.
the rta results were calculated using the cross-validation sets from ref.
<cit> .
the smm-align results were taken from the same reference.
"✔" indicates a correct prediction and "✘" indicates an incorrect one.
"na" denotes that a prediction is not available either because the method does not have a prediction model for that mhc type or the method cannot be applied to a peptide shorter than  <dig> residues.
creating independent cross-validation sets
the most relevant measure of a prediction model's accuracy is how well it performs on novel peptide sequences that are dissimilar to those used to fit the model.
a model with a large number of parameters, such as most peptide-mhc binding prediction models, can potentially overfit the data so that its accuracy is good on the training set but poor on unrelated peptide sequences.
cross-validation is one procedure that can be used to estimate the accuracy of the prediction model for novel data.
we have used 5-fold cross-validation in which the data is divided into  <dig> approximately equal size subsets and predictions are made in turn on each subset using a model fit to the remaining data in the other  <dig> subsets.
because many of the peptide sequences in the experimental binding data are similar it is important to create the cross-validation sets with the fewest possible sequence similarities between subsets.
the binding data from the iedb database was divided into 5-fold cross-validation data sets using a method that minimizes the number of common 9-mer subsequences shared between different cross-validation sets while keeping the size of the  <dig> sets approximately the same size.
the method was able to achieve perfect separation, i.e.
no common 9-mers between sets, for all mhc allotypes except for hla-drb1* <dig>  which had only  <dig> common subsequences .
we briefly outline the method here; details are given in additional file  <dig>  the peptide sequences are first divided into sets of similar sequences such that every sequence within a given set shares a common 9-mer with at least one other sequence in the same set but does not share any common 9-mers with sequences in other sets.
if these sets are sufficiently small then perfectly separated 5-fold cross-validation sets can be created simply by partitioning these sequence sets into  <dig> parts, each part containing sets with approximately 1/ <dig> of the total number of sequences.
because there are no common 9-mers between sequence sets and because each sequence set is entirely contained within a cross-validation set there are no common 9-mers between cross-validation sets.
the only impediment to generating perfectly separated cross-validation sets, without any inter-set sequence similarities, is if some of the sets contain more than 1/ <dig> of the sequences and so must be split.
the peptide sequences are divided into sets sharing common 9-mers using a graph representation of their similarities.
the nodes in the graph represent sequences and edges connect nodes for sequences sharing a common 9-mer.
the sequences corresponding to the nodes in each connected component of this graph are then the initial sets of similar sequences referred to above.
if one of these connected subgraphs is too large then it needs to be split in two in a way that minimizes the number of similar sequence similarities between the two new parts.
since edges represent sequence similarities, this is accomplished by finding a graph cut that has few edges between the two new subgraphs.
this graph cut problem is formulated as the sparsest cut problem in which the goal is to find a set of nodes s that minimizes the sparsity sp defined by
in which the numerator is the number of edges between nodes in s and nodes outside of s and the denominator is the minimum number of nodes in either subset.
this measure insures that not only does the cut cross few edges, it also divides the graph approximately in half.
it would be easy, but not very productive, to find a cut that separates only one node, for example.
furthermore we chose to solve this np-hard problem using an approximate algorithm based on spectral graph theory.
although there are alternative algorithms that achieve tighter worst-case approximation bounds, the spectral method has worked well in applications.
also it is much easier to implement, since it only requires calculating matrix eigenvalues and eigenvectors, for which we used svd.
the approximate solution of the sparsest cut problem is based on the bound on the sparsity
in which dmax is the maximum node degree and λ <dig> is the second smallest eigenvalue of the graph laplacian matrix l defined by
with d the degree of node i.
this bound implies that an efficient cut, with low sparsity, is guaranteed if λ <dig> is small  <cit> .
furthermore, the cut can be found using the associated eigenvector v <dig>  or so-called fiedler vector.
the cut is defined by dividing the graph into subgraphs s and  according to the eigenvector components with
for some cutoff c  <cit> .
the approximate solution to the sparsest cut problem is found by calculating the sparsity for all |v|- <dig> such cuts, defined by sc with cutoff c equal to each vector component, and choosing the one with the minimimum sp.
the overall algorithm to find optimal cross-validation sets begins by grouping the peptide sequences by similarity according to connected components or, in cases in which the groups are too large, connected components split using the spectral partitioning algorithm just described.
the groups only needed to be split for hla-drb1* <dig> and only a single iteration of splitting was necessary.
each group is then iteratively added to the cross-validation set with the fewest sequences starting with the largest groups.
the total number of pairs of similar sequences between the training and test sets, which ideally should be small, is equal to the total number of edges removed in the sparsest cut portion of the algorithm.
this means that perfect separation is achieved for all mhc allotypes except for hla-drb1* <dig>  the generated cross-validation sets are provided in additional file  <dig>
RESULTS
binding affinity prediction
the rta prediction method was evaluated both with the 5-fold cross-validation sets in the smm-align reference  <cit> , so that its performance could be compared with the results for the  <dig> different prediction methods examined in that paper, and also with the generally larger 5-fold cross-validation sets created using the method described above, which minimizes the overlap between training and test set peptide sequences.
the prediction performance was evaluated using the area under the receiver operating characteristic  curve  <cit> , or auc, the root mean squared error , and the pearson correlation coefficient.
the latter two statistics were not calculated for the smm-align data sets because the data for low affinity binders in those data sets were assigned a fixed cutoff value, which renders those statistics less meaningful.
the prediction results for the new data sets are given in table  <dig> and the results for the smm-align data sets are in table  <dig>
it may be seen from table  <dig> that rta had better performance than smm-align on the same data sets for all mhc allotypes examined.
the improvement is especially large for hla-drb3* <dig>  with an auc for rta of  <dig>  compared with only  <dig>  for smm-align.
however this smm-align data set is exceptional since it contains only  <dig> strong binders, with ic <dig> <  <dig> nm, and has the highest percentage  of thresholded low-affinity values .
ref.
<cit>  also provides auc values obtained by  <dig> other prediction methods applied to the same smm-align data sets.
the methods examined were a gibbs sampler  <cit> , tepitope  <cit> , svrmhc  <cit> , mhcpred  <cit> , arb  <cit> , netmhcii  <cit> , predbalb/c  <cit> , and a variant of smm-align that includes information on the peptide length and flanking residues called smm-prf  <cit> .
a total of  <dig> of the methods could make predictions for all  <dig> hla-dr allotypes, while  <dig> other methods could only be used for a subset of these allotypes.
only netmhcii could be used for all  <dig> murine allotypes while two other methods, arb and predbalb/c, could be applied to some of these allotypes.
overall, rta was more accurate than the other  <dig> prediction methods.
rta yielded the same auc value as tepitope for hla-drb1* <dig> and the same value as arb for hla-drb1* <dig>  it gave a lower auc value than one of other the other methods only for hla-drb1* <dig>  hla-drb1* <dig>  and h2-ias.
significantly, rta yielded the highest auc values for all  <dig> methods examined in that study for  <dig> out of the  <dig> allotypes, even though  <dig> of the methods were evaluated without cross-validation so that some of the training data may have been included in the test sets.
in general, the auc values for the iedb data sets are slightly lower than for the corresponding smm-align data sets for the same mhc allotypes.
one factor that likely contributes to this difference is that the composition of the iedb and smm-align data sets are quite different.
whereas the median percentage of low affinity peptides, with ic <dig> >  <dig> μm, in the iedb data set is only  <dig> %, it is 52% for the smm-align sets.
one common trend between the results for the iedb and smm-align data sets is that the auc for the murine mhc allotypes is generally higher than for the human ones.
the rta prediction results for the iedb data sets were not compared with results using other prediction methods because these other methods were fit to peptide-mhc binding data that likely significantly overlap the iedb data, especially considering that iedb is one of the most comprehensive databases.
the occurrence of similar test set data in the training data can significantly increase the apparent prediction accuracy, leading to overly optimistic accuracy estimates.
in the extreme case in which the rta model was trained on both the training and test set data, close to perfect auc values of >  <dig>  were obtained .
this also emphasizes the importance of using minimally similar training and test data sets in order to reliably estimate the prediction accuracy for novel peptide sequences.
the iedb cross-validation sets, created using the algorithm described above, had no overlap at all, defined by identical 9-residue segments, for all but one mhc allotype.
likewise, the smm-align cross-validation data sets used for comparison with that method had low overlap of about  <dig> -2%.
finally, we point out that using more cross-validation sets, will allow a larger fraction of the limited data to be used for training but it will also generally increase the overlap between corresponding training and test sets.
in the extreme case of leave-one-out cross-validation this overlap will be large, leading to overly optimistic performance estimates, unless all similar peptide sequences are removed.
the rta method was also evaluated by making predictions for the peptide binding data in lin et al.
<dig>  <cit> .
that study evaluated  <dig> class ii mhc binding prediction servers using binding affinity data for  <dig> peptides binding to seven different allotypes.
the peptides comprise overlapping segments of four allergens and were chosen as an unbiased test set for comparing different binding prediction methods.
rta predictions were made for all allotypes except for hla-drb1* <dig>  for which insufficient binding data was available for training.
the prediction performance, as evaluated by auc values, is shown in table  <dig>  a comparison with the results for the prediction servers in ref.
<cit>  shows that rta yields among the most accurate results for hla-drb1* <dig>  drb1* <dig>  drb1* <dig>  and drb1* <dig> and about average results for hla-drb1* <dig> and hla-drb1* <dig>  one potential uncertainty in this evaluation is how many of the peptides in this test set were used to train the individual prediction methods.
as mentioned above, significant overlap between the test and training sets is expected to give an overly optimistic assessment of accuracy.
this problem can be reduced through cross-validation, however this procedure is infeasible for evaluating prediction servers on the internet.
finally, the rta method was evaluated using the iedb similarity reduced data sets from el-manzalawy et al.
<dig>  <cit> .
these data sets were compiled from binding data downloaded from the iedb database by removing similar peptides at different similarity thresholds: unique sequences , no common 9-mer subsequences , and pairwise sequence identity < 80% .
the prediction results for rta as well as the three prediction methods evaluated in ref.
<cit>  are given in tables s1-s <dig> in additional file  <dig>  these results show that the rta method performed significantly better than the other three prediction methods.
rta had the highest auc out of all four methods for 9/ <dig> mhc allotypes for the upds data set and the highest auc for all  <dig> mhc allotypes for both the srds <dig> and srds <dig> sets.
peptide binding register prediction
the peptide-mhc binding prediction model is based on the observation in x-ray structures of complexes that only a portion of the peptide binds to the mhc.
the model predicts peptide binding at ambient temperature by summing the contributions of all possible 9-mer segments binding to the mhc.
the difference between the binding affinity of the strongest binding 9-mer segment and that of other possible binding segments will be larger at the low temperatures used to determine x-ray structures and presumably is represented in the structure of the complex.
the available x-ray structures of peptide-mhc complexes can be used to check this assumption by comparing the prediction of the binding register, defined by the largest term in the boltzmann-weighted average in eq.
<dig>  with the actual binding register observed in the x-ray crystal structure.
the peptide binding registers inferred from the x-ray structures provide additional experimental data that is unrelated to the binding affinity data and so allow an independent estimate of the model's accuracy.
finally, we note that there has been some confusion in the literature since the sequence of the peptide appearing in the x-ray crystal structure may in fact be different from the sequence of the actual co-crystallized peptide.
this is because the peptide n- and c-termini may be disordered and so missing from the structure.
in addition, some peptide residues with disordered side chains are modeled as alanines in the structure.
because of this potential discrepancy the peptide sequences in table  <dig> were obtained from papers reporting the pdb structures rather than the structure itself.
the 9-residue core segment of the peptide that binds to the mhc molecule, which begins with the p <dig> residue, is highlighted in boldface in the peptide sequence.
these structures show which peptide segment binds to the mhc molecule and so are used to evaluate how well the rta method predicts the correct peptide binding register.
peptide binding register predictions were made for all unique peptide-mhc complexes with available x-ray structures.
these are shown in table  <dig>  the 9-mer segment in the structure begins with the residue binding in the p <dig> pocket and the predicted 9-mer segment extends from residue mmax + <dig> to mmax +  <dig> in which mmax is the m index of the largest term in the boltzmann sum appearing in eq.
<dig>  all peptides similar to those for which the predictions were made were removed from the training set.
the results showed that the rta method predicted the peptide binding registers accurately.
the predictions for all peptides in table  <dig> were correct except for one, the peptide rgisqavhaahaeineagr binding to h2-iad from pdb entry 1iao.
interestingly, the core segment for the clip peptide is the same for binding to both hla-drb1* <dig>  and h2-iab , as has been observed for most other human and murine mhc allotypes  <cit> .
clip binds to all class ii mhc molecules as an intermediate step in mhc processing and peptide loading.
for comparison, peptide binding register predictions were made for the peptide-mhc complexes listed in table  <dig> using seven other prediction methods implemented as web servers: arb  <cit> , mhc2pred  <cit> , nn-align  <cit> , propred  <cit> , rankpep  <cit> , smm-align  <cit> , and svrmhc  <cit> .
each of these methods uses different computational approaches to predicting peptide-mhc binding affinities.
arb, propred, and rankpep use position specific scoring matrices , svrmhc and mhc2pred use support vector machines , nn-align uses artificial neural networks , and smm-align uses a combination of scoring matrices and a novel gibbs sampler method  <cit> .
the overall prediction results for rta and these other methods are given in table  <dig>  the performance of the different prediction methods could not be directly compared since most methods could not be applied to all peptide-mhc complexes.
however, rta is clearly among the most accurate methods for this data set, with  <dig> out of the  <dig> peptide binding cores correctly predicted.
a recent study by wang et al.
<cit>  performed a similar analysis using an overlapping but slightly different set of peptide-mhc complexes and obtained qualitatively similar results for the same prediction methods examined in that study.
that study concluded that the tepitope method  <cit> , as implemented by the propred web server  <cit> , achieved the best overall performance.
we found similar results, with only rta and propred correctly predicting all  <dig> hla-dr peptide binding registers.
nn-align achieved similar performance with  <dig> correctly predicted hla-dr peptide binding cores.
significantly, some prediction methods, which generally attain accuracies that are comparable to or better than the earlier tepitope method in predicting peptide binding affinities, had significantly worse accuracy than tepitope in predicting the peptide binding cores.
this demonstrates that a prediction model that achieves high accuracy for peptide binding affinity does not necessarily achieve corresponding high accuracy for predicting the peptide binding core.
the fact that our rta method does perform well for both prediction tasks suggests that the underlying approximate physical model for the contribution of each binding register to the total binding affinity is correct and accurate.
importance of suboptimal terms in the boltzmann sum
a common approximation in other prediction models, such as tepitope, is to use only the peptide binding register with the highest affinity and so neglect the contributions of suboptimal binding registers.
under this approximation eq.
<dig> would become
this approximation is accurate if only one term dominates the boltzmann average in eq.
<dig>  in order to estimate the importance of the suboptimal terms included in our model, we calculated the ratios  and |t2/t1|, in which ti are the terms in the boltzmann sum in eq.
<dig> ordered so that |t1| ≥ |t2| ≥ l ≥ |tl-8|.
t <dig>  are then the terms with the largest  absolute values.
both ratios vary from  <dig>  to  <dig> .
these ratios were calculated for the largest iedb data set, containing  <dig> peptides binding to hla-drb1* <dig>  a low value of  indicates that suboptimal terms make a significant contribution to the sum so that the approximation using only the optimal term is inaccurate.
the cumulative distribution of this quantity is shown in figure  <dig>  for example, suboptimal terms contribute at least 50% to the boltzmann sum for about 20% of the peptide sequences.
likewise figure  <dig> shows the cumulative distribution for the relative contribution of the leading suboptimal term relative to the optimal term, |t2/t1|.
for example, it can be seen that this ratio is greater than  <dig>  for about 40% of the sequences.
the ratios for other data sets are expected to have a similar distribution.
to the extent that these predicted values are accurate, the distributions of these ratios indicate that suboptimal terms make significant contributions to the total binding affinity for many peptides so that the assumption that only one binding register dominates is incorrect for these peptides.
the error in the predicted binding affinity for these peptides in a prediction model including only the optimal term may not be as large as expected since the model parameters are fit directly to experimental binding affinities and so the model may compensate to some degree for the contribution of missing suboptimal terms.
in any case, to the extent that these estimates of the contributions of suboptimal terms are correct, prediction models, such as rta, that incorporate these terms are expected to be more accurate.
CONCLUSIONS
the prediction results for the rta method show  that suboptimal peptide binding registers make a significant contribution to the total peptide binding affinity for many peptides,  regularization with the l <dig> parameter constraint improves prediction accuracy for novel data,  the method achieves good accuracy in predicting peptide-mhc binding affinities with generally higher auc values than other prediction methods evaluated on the same data sets, and  the method accurately predicts the 9-residue binding core for all but one out of  <dig> peptide-mhc complexes examined.
while peptide-mhc binding affinities can be quickly calculated using eq.
<dig>  the efficient solution of the difficult global optimization problem for fitting parameters using local optimization methods contributed to the success of the method.
the accuracy of the method is expected to continually improve as more experimental binding data become available, particularly for certain mhc allotypes with limited current data.
one direction for future work is to extend the rta model to a multi-allotype model by including another parameter dimension reflecting the mhc peptide binding determinants, either through unique binding pockets, as in tepitope, or through polymorphic pocket residues, as in netmhciipan  <cit> .
another possible extension is to include interaction terms for nearby peptide residues.
the approximation that contributions of individual peptide residues to the total binding free energy are largely independent is supported by side chain scanning experiments  <cit>  and by the accuracy of prediction models, such as rta, that are based on this approximation.
however, recent experimental studies  <cit>  have found evidence for cooperative effects between peptide residues so that additional pairwise interaction terms may improve prediction accuracy.
both of these modifications will significantly increase the number of model parameters, making regularization even more crucial for prediction accuracy.
authors' contributions
ajb conceived of the study, collected the data sets, and analyzed the prediction results.
hdm performed the numerical optimization to determine the model parameters.
both authors participated in drafting the manuscript and read and approved the final manuscript.
supplementary material
additional file 1
detailed description of the algorithm for creating cross-validation data sets with minimal peptide sequence overlap.
click here for file
additional file 2
this zip archive contains tab-separated files with the 5-fold cross-validation sets for  <dig> different mhc allotypes created for this study and also a separate text file describing the table format.
click here for file
additional file 3
tables s1-s <dig> present the prediction results for the iedb benchmark sets from el-manzawaly et al.
<dig>  <cit> .
prediction accuracies are given for rta as well as the three prediction methods studied in that paper.
click here for file
acknowledgements
the work of ajb was supported by the mayo clinic.

 class ii mhc is an essential protein complex in the adaptive immune system that is involved in activating helper t cells. it is found on the surfaces of specialized immune system cells, such as dendritic cells and macrophages, where it binds fragments of extracellular peptides and presents them to cd4+ helper t cells. class ii mhc is highly polymorphic and each allotype generally has different peptide binding preferences. the computational prediction of such binding affinities can be applied to the important problem of finding promiscuous epitopes that bind to multiple allotypes for use in rational vaccine design. also these predictions may aid in understanding autoimmune diseases, many of which have been linked to particular class ii mhc alleles. characteristic properties of peptide binding to class ii mhc molecules can be inferred from available high-resolution x-ray structures of peptide-mhc complexes. the binding cleft in class ii mhc is open at both ends so that it generally binds longer peptides, 15- residues in length, with the peptide n- and c-terminii extending out of the cleft. this is in contrast to class i mhc, which has a peptide binding cleft that is closed at both ends so that it only accommodates shorter peptides of 8- residues. an examination of x-ray crystal structures of peptide bound class ii mhc molecules reveals that the peptides bind in an extended polyproline ii helix with the conserved backbone structure maintained by conserved hydrogen bonds with the mhc. peptide side chains within a contacting nine residue core segment of the peptide bind into pockets on the mhc surface. many of the polymorphic residues occur in these pockets and so determine the distinctive peptide binding preferences of each mhc allotype. although the core residues of the bound peptide assume similar backbone conformations, each peptide can generally bind in multiple registers. this makes the prediction of peptide binding affinities for class ii mhc particularly challenging since, unlike class i mhc, all possible binding registers must be considered. prediction methods for class ii epitopes remains an active area of research, motivated both by their biomedical importance and the difficulty of achieving high accuracy. early prediction methods fit the total peptide binding energy, binding motif, geometric average binding affinity, or sequence alignment profile in a particular register to a linear combination of contributions from individual residues, and represented them as binding profile matrices. the scores for all possible peptide binding registers were calculated and either the maximum value or sum were used as a total peptide binding score. later methods employed various machine learning and data fitting approaches to prediction including partial least squares, gibbs sampling, linear programming, support vector machines, and kernel methods. one prediction method, called smm-align, successfully combined two methods, binding profile matrices and gibbs sampling. we introduce a new model for predicting binding affinities of peptides to class ii mhc called the regularized thermodynamic average, or rta, model. three features are incorporated into the model: independent binding of individual peptide residues so that the total binding energy is a sum of the binding energies of each core peptide residue, a thermodynamic average over all possible peptide binding registers, and regularization of the model to prevent overfitting. feature # the approximate independence of peptide binding on peptide residue interactions, is supported experimentally and is a basic assumption in most prediction models to date. indeed, it is the primary basis for binding profile methods. the importance of contributions of suboptimal peptide binding registers to the total binding affinity motivates feature # experimental studies have found such cases in which a peptide binds with significant affinity to class ii mhc in multiple registers. in fact, such peptide binding in multiple registers may be involved in the origin of autoimmune diseases. the thermodynamic average in the rta model insures that contributions from suboptimal binding registers, which may be significant, also contribute to the total binding affinity. regularization is especially important in prediction models that contain a relatively large number of parameters, such as ours. regularization is implemented in rta using an l penalty function, like in lasso regression. this has the advantage over alternative penalty functions, such as the l penalty employed in smm-align, that many coefficients vanish at lower cutoff values so that the number of model parameters is reduced. in addition, unlike in many previous prediction methods, the p residue for peptides binding to hla-dr is not restricted to be hydrophobic in the rta model. finally, we use deterministic local optimization methods to fit the model parameters, which we found to converge to better solutions than commonly employed stochastic optimization methods. the rta prediction method was tested using cross-validation sets from an smm-align study, for comparison with results from different prediction methods examined in that study, and also using larger cross-validation sets created from binding data from the iedb database. the rta method was also evaluated using additional peptide binding data sets from lin et al. and el-manzalawy et al.. a novel method to efficiently create data sets with minimal sequence overlap was developed to create the new cross-validation sets. the performance of rta in predicting peptide binding affinities and in predicting the correct peptide core binding residues were tested. finally, we examined the relative contribution of suboptimal peptide binding registers to the overall binding affinity. methods peptide-mhc binding data set peptide binding data for the mhc allotypes in table were downloaded from the iedb database http://www.immuneepitope.org. all ic values from either radioactively or fluorescently labeled peptide competition binding assays were collected. the median value was used for the few cases in which multiple experimental values were measured. the total numbers of data for each allotype are given in table these datasets were then divided into minimally overlapping 5-fold cross-validation sets using the method described in "creating independent cross-validation sets". prediction model as mentioned above, the prediction of peptide - class ii mhc binding affinity involves not only estimating the affinity of the bound core portion of the peptide but also accounting for multiple binding registers. the model begins with the approximation that the contributions to the total binding free energy of binding individual peptide residues to each mhc pocket are approximately independent. this assumption has been tested experimentally and is the basis of matrix-based prediction methods in which the total peptide-mhc binding energy is a sum of pocket specific residue binding energies. furthermore, x-ray structures of peptide - class ii mhc complexes show that the 9-residue core bound portions of the peptides adopt similar backbone conformations so that the total score over an ungapped 9-residue segment of the peptide sequence is appropriate. the peptide-mhc binding data from competition experiments are given as ic values. ic is not a direct measure of binding affinity but it can be related to the inhibition constant, ki, through the cheng-prusoff equation is the concentration of the free labeled peptide and kd is its binding affinity to the mhc molecule. it is assumed that the labeled peptide concentration is sufficiently low, i.e. = kd, so that ki ≈ ic this is an implicit assumption in other prediction methods as well. according to this approximation, experimental binding affinities were calculated as Δgexp = - kt log, in which kt is the boltzmann constant times the absolute temperature, which is approximately kcal/mol at room temperature. each peptide can bind the mhc molecule in different registers defined by which 9-residue core peptide segment contacts the mhc. the length of peptide k is denoted by l. the register will be denoted by the index m, which varies from for the first residues as the core segment to l - for the last residues as the core segment. the sequence of peptide k is represented by matrix elements, which are equal to if residue type j appears at position i in the peptide sequence and are equal to otherwise. the affinity of peptide k binding to the mhc molecule in register m is then the linear function in which the coefficients βij are to be fit using the experimental data. the apparent total binding affinity for peptide k can be calculated as in which is the fraction of bound peptides binding in register m and is the affinity for binding in register m, calculated in the model using eq. an entropic contribution to the free energy in eq. was omitted because it did not improve prediction performance. the fraction of peptide bound in a particular register is related to the binding affinities in eq. by,, and are the concentrations of unbound peptide, unbound mhc and peptide bound to the mhc in register m, respectively. is the dissociation constant for peptide k binding to the mhc in register m. by combining eqs. and one arrives at an expression for the total binding affinity as an exponentially-weighted average of the affinities for the peptide binding in different registers many prediction methods for peptide-mhc binding make the further approximation that the largest binding affinity dominates the sum in eq. so that only the largest binding affinity over all possible registers is included in the prediction model. because suboptimal terms are comparable in magnitude to the optimal term in many cases, we will retain all terms in eq. the final component in the prediction model is regularization. this is needed to avoid overfitting, since the model contains parameters, which is relatively large compared to the number of binding data for each mhc type. regularization is accomplished by imposing an upper bound on the l norm of the coefficient amplitudes with the cutoff parameter t. this type of constraint is also used in lasso regression and has the desirable property that as the cutoff t is lowered, an increasing number of coefficients become zero so that the number of nonzero parameters in the model is reduced. by combining these components, the binding affinity for peptide k, Δg, is calculated as the parameters in this model, {βij, i = k, j = k, 20}, are fit to the experimental binding affinities,, by minimizing the mean square error subject to the constraint in eq. optimization method to solve the optimization problem minimizing eq. under the constraint in eq. one first splits the variable βij into the difference of two nonnegative variables and, so that βij = -. the constraint eq. then has as the left side the sum over both new variables. overall one has a linearly constrained nonlinear and nonconvex optimization problem. in principle, global optimization methods would have to be applied to find the global minimizer. due to the dimensions of the problems these would be non-deterministic methods such as some of the many metaheuristics, which include simulated annealing and genetic algorithms. these methods require a large number of evaluations and due to their stochastic character would have to be run several times in order to increase the likelihood of finding the global optimum, although a guarantee for that is impossible. after initial tests we decided to instead use local solvers. these are very efficient and through various measures the chances of getting very good local minima can be increased substantially. using local solvers with "multistart" or several often randomly generated starting guesses is another way of solving global optimization problems. it is also implemented in several software packages. we instead used simple constant starting guesses. all of our results were obtained using the constant starting values of or for all components of β+ and β-. we did, however, use a sequence of different t values in the regularization constraint eq. and then chose the best minimum found for the solvers used. it is significant that the best solutions were obtained with values of the regularization parameter t that were small enough to restrict the fit βij values. this shows that the additional l constraint in eq. helps alleviate overfitting. in order to be able to easily call a variety of solvers we phrased the problem in the modeling language ampl. in order not to have to list many separate citations we state that we used the applicable solvers installed at neos but run locally, not through this free service in which we are also heavily involved. in the way described above we generated the values in tables and for table it was possible to just run the solver that had best solved the corresponding case in table on the data set and evaluate for the test sets. we did not use the optimal coefficients found but applied our method to the full training set in each case. the rta results were calculated using the cross-validation sets from ref.. the smm-align results were taken from the same reference. "✔" indicates a correct prediction and "✘" indicates an incorrect one. "na" denotes that a prediction is not available either because the method does not have a prediction model for that mhc type or the method cannot be applied to a peptide shorter than residues. creating independent cross-validation sets the most relevant measure of a prediction model's accuracy is how well it performs on novel peptide sequences that are dissimilar to those used to fit the model. a model with a large number of parameters, such as most peptide-mhc binding prediction models, can potentially overfit the data so that its accuracy is good on the training set but poor on unrelated peptide sequences. cross-validation is one procedure that can be used to estimate the accuracy of the prediction model for novel data. we have used 5-fold cross-validation in which the data is divided into approximately equal size subsets and predictions are made in turn on each subset using a model fit to the remaining data in the other subsets. because many of the peptide sequences in the experimental binding data are similar it is important to create the cross-validation sets with the fewest possible sequence similarities between subsets. the binding data from the iedb database was divided into 5-fold cross-validation data sets using a method that minimizes the number of common 9-mer subsequences shared between different cross-validation sets while keeping the size of the sets approximately the same size. the method was able to achieve perfect separation, i.e. no common 9-mers between sets, for all mhc allotypes except for hla-drb1* which had only common subsequences. we briefly outline the method here; details are given in additional file the peptide sequences are first divided into sets of similar sequences such that every sequence within a given set shares a common 9-mer with at least one other sequence in the same set but does not share any common 9-mers with sequences in other sets. if these sets are sufficiently small then perfectly separated 5-fold cross-validation sets can be created simply by partitioning these sequence sets into parts, each part containing sets with approximately 1/ of the total number of sequences. because there are no common 9-mers between sequence sets and because each sequence set is entirely contained within a cross-validation set there are no common 9-mers between cross-validation sets. the only impediment to generating perfectly separated cross-validation sets, without any inter-set sequence similarities, is if some of the sets contain more than 1/ of the sequences and so must be split. the peptide sequences are divided into sets sharing common 9-mers using a graph representation of their similarities. the nodes in the graph represent sequences and edges connect nodes for sequences sharing a common 9-mer. the sequences corresponding to the nodes in each connected component of this graph are then the initial sets of similar sequences referred to above. if one of these connected subgraphs is too large then it needs to be split in two in a way that minimizes the number of similar sequence similarities between the two new parts. since edges represent sequence similarities, this is accomplished by finding a graph cut that has few edges between the two new subgraphs. this graph cut problem is formulated as the sparsest cut problem in which the goal is to find a set of nodes s that minimizes the sparsity sp defined by in which the numerator is the number of edges between nodes in s and nodes outside of s and the denominator is the minimum number of nodes in either subset. this measure insures that not only does the cut cross few edges, it also divides the graph approximately in half. it would be easy, but not very productive, to find a cut that separates only one node, for example. furthermore we chose to solve this np-hard problem using an approximate algorithm based on spectral graph theory. although there are alternative algorithms that achieve tighter worst-case approximation bounds, the spectral method has worked well in applications. also it is much easier to implement, since it only requires calculating matrix eigenvalues and eigenvectors, for which we used svd. the approximate solution of the sparsest cut problem is based on the bound on the sparsity in which dmax is the maximum node degree and λ is the second smallest eigenvalue of the graph laplacian matrix l defined by with d the degree of node i. this bound implies that an efficient cut, with low sparsity, is guaranteed if λ is small. furthermore, the cut can be found using the associated eigenvector v or so-called fiedler vector. the cut is defined by dividing the graph into subgraphs s and according to the eigenvector components with for some cutoff c. the approximate solution to the sparsest cut problem is found by calculating the sparsity for all |v|- such cuts, defined by sc with cutoff c equal to each vector component, and choosing the one with the minimimum sp. the overall algorithm to find optimal cross-validation sets begins by grouping the peptide sequences by similarity according to connected components or, in cases in which the groups are too large, connected components split using the spectral partitioning algorithm just described. the groups only needed to be split for hla-drb1* and only a single iteration of splitting was necessary. each group is then iteratively added to the cross-validation set with the fewest sequences starting with the largest groups. the total number of pairs of similar sequences between the training and test sets, which ideally should be small, is equal to the total number of edges removed in the sparsest cut portion of the algorithm. this means that perfect separation is achieved for all mhc allotypes except for hla-drb1* the generated cross-validation sets are provided in additional file  binding affinity prediction the rta prediction method was evaluated both with the 5-fold cross-validation sets in the smm-align reference, so that its performance could be compared with the results for the different prediction methods examined in that paper, and also with the generally larger 5-fold cross-validation sets created using the method described above, which minimizes the overlap between training and test set peptide sequences. the prediction performance was evaluated using the area under the receiver operating characteristic curve, or auc, the root mean squared error, and the pearson correlation coefficient. the latter two statistics were not calculated for the smm-align data sets because the data for low affinity binders in those data sets were assigned a fixed cutoff value, which renders those statistics less meaningful. the prediction results for the new data sets are given in table and the results for the smm-align data sets are in table it may be seen from table that rta had better performance than smm-align on the same data sets for all mhc allotypes examined. the improvement is especially large for hla-drb3* with an auc for rta of compared with only for smm-align. however this smm-align data set is exceptional since it contains only strong binders, with ic < nm, and has the highest percentage of thresholded low-affinity values. ref. also provides auc values obtained by other prediction methods applied to the same smm-align data sets. the methods examined were a gibbs sampler, tepitope, svrmhc, mhcpred, arb, netmhcii, predbalb/c, and a variant of smm-align that includes information on the peptide length and flanking residues called smm-prf. a total of of the methods could make predictions for all hla-dr allotypes, while other methods could only be used for a subset of these allotypes. only netmhcii could be used for all murine allotypes while two other methods, arb and predbalb/c, could be applied to some of these allotypes. overall, rta was more accurate than the other prediction methods. rta yielded the same auc value as tepitope for hla-drb1* and the same value as arb for hla-drb1* it gave a lower auc value than one of other the other methods only for hla-drb1* hla-drb1* and h2-ias. significantly, rta yielded the highest auc values for all methods examined in that study for out of the allotypes, even though of the methods were evaluated without cross-validation so that some of the training data may have been included in the test sets. in general, the auc values for the iedb data sets are slightly lower than for the corresponding smm-align data sets for the same mhc allotypes. one factor that likely contributes to this difference is that the composition of the iedb and smm-align data sets are quite different. whereas the median percentage of low affinity peptides, with ic > μm, in the iedb data set is only %, it is 52% for the smm-align sets. one common trend between the results for the iedb and smm-align data sets is that the auc for the murine mhc allotypes is generally higher than for the human ones. the rta prediction results for the iedb data sets were not compared with results using other prediction methods because these other methods were fit to peptide-mhc binding data that likely significantly overlap the iedb data, especially considering that iedb is one of the most comprehensive databases. the occurrence of similar test set data in the training data can significantly increase the apparent prediction accuracy, leading to overly optimistic accuracy estimates. in the extreme case in which the rta model was trained on both the training and test set data, close to perfect auc values of > were obtained. this also emphasizes the importance of using minimally similar training and test data sets in order to reliably estimate the prediction accuracy for novel peptide sequences. the iedb cross-validation sets, created using the algorithm described above, had no overlap at all, defined by identical 9-residue segments, for all but one mhc allotype. likewise, the smm-align cross-validation data sets used for comparison with that method had low overlap of about -2%. finally, we point out that using more cross-validation sets, will allow a larger fraction of the limited data to be used for training but it will also generally increase the overlap between corresponding training and test sets. in the extreme case of leave-one-out cross-validation this overlap will be large, leading to overly optimistic performance estimates, unless all similar peptide sequences are removed. the rta method was also evaluated by making predictions for the peptide binding data in lin et al.. that study evaluated class ii mhc binding prediction servers using binding affinity data for peptides binding to seven different allotypes. the peptides comprise overlapping segments of four allergens and were chosen as an unbiased test set for comparing different binding prediction methods. rta predictions were made for all allotypes except for hla-drb1* for which insufficient binding data was available for training. the prediction performance, as evaluated by auc values, is shown in table a comparison with the results for the prediction servers in ref. shows that rta yields among the most accurate results for hla-drb1* drb1* drb1* and drb1* and about average results for hla-drb1* and hla-drb1* one potential uncertainty in this evaluation is how many of the peptides in this test set were used to train the individual prediction methods. as mentioned above, significant overlap between the test and training sets is expected to give an overly optimistic assessment of accuracy. this problem can be reduced through cross-validation, however this procedure is infeasible for evaluating prediction servers on the internet. finally, the rta method was evaluated using the iedb similarity reduced data sets from el-manzalawy et al.. these data sets were compiled from binding data downloaded from the iedb database by removing similar peptides at different similarity thresholds: unique sequences, no common 9-mer subsequences, and pairwise sequence identity < 80%. the prediction results for rta as well as the three prediction methods evaluated in ref. are given in tables s1-s in additional file these results show that the rta method performed significantly better than the other three prediction methods. rta had the highest auc out of all four methods for 9/ mhc allotypes for the upds data set and the highest auc for all mhc allotypes for both the srds and srds sets. peptide binding register prediction the peptide-mhc binding prediction model is based on the observation in x-ray structures of complexes that only a portion of the peptide binds to the mhc. the model predicts peptide binding at ambient temperature by summing the contributions of all possible 9-mer segments binding to the mhc. the difference between the binding affinity of the strongest binding 9-mer segment and that of other possible binding segments will be larger at the low temperatures used to determine x-ray structures and presumably is represented in the structure of the complex. the available x-ray structures of peptide-mhc complexes can be used to check this assumption by comparing the prediction of the binding register, defined by the largest term in the boltzmann-weighted average in eq. with the actual binding register observed in the x-ray crystal structure. the peptide binding registers inferred from the x-ray structures provide additional experimental data that is unrelated to the binding affinity data and so allow an independent estimate of the model's accuracy. finally, we note that there has been some confusion in the literature since the sequence of the peptide appearing in the x-ray crystal structure may in fact be different from the sequence of the actual co-crystallized peptide. this is because the peptide n- and c-termini may be disordered and so missing from the structure. in addition, some peptide residues with disordered side chains are modeled as alanines in the structure. because of this potential discrepancy the peptide sequences in table were obtained from papers reporting the pdb structures rather than the structure itself. the 9-residue core segment of the peptide that binds to the mhc molecule, which begins with the p residue, is highlighted in boldface in the peptide sequence. these structures show which peptide segment binds to the mhc molecule and so are used to evaluate how well the rta method predicts the correct peptide binding register. peptide binding register predictions were made for all unique peptide-mhc complexes with available x-ray structures. these are shown in table the 9-mer segment in the structure begins with the residue binding in the p pocket and the predicted 9-mer segment extends from residue mmax + to mmax + in which mmax is the m index of the largest term in the boltzmann sum appearing in eq. all peptides similar to those for which the predictions were made were removed from the training set. the results showed that the rta method predicted the peptide binding registers accurately. the predictions for all peptides in table were correct except for one, the peptide rgisqavhaahaeineagr binding to h2-iad from pdb entry 1iao. interestingly, the core segment for the clip peptide is the same for binding to both hla-drb1* and h2-iab, as has been observed for most other human and murine mhc allotypes. clip binds to all class ii mhc molecules as an intermediate step in mhc processing and peptide loading. for comparison, peptide binding register predictions were made for the peptide-mhc complexes listed in table using seven other prediction methods implemented as web servers: arb, mhc2pred, nn-align, propred, rankpep, smm-align, and svrmhc. each of these methods uses different computational approaches to predicting peptide-mhc binding affinities. arb, propred, and rankpep use position specific scoring matrices, svrmhc and mhc2pred use support vector machines, nn-align uses artificial neural networks, and smm-align uses a combination of scoring matrices and a novel gibbs sampler method. the overall prediction results for rta and these other methods are given in table the performance of the different prediction methods could not be directly compared since most methods could not be applied to all peptide-mhc complexes. however, rta is clearly among the most accurate methods for this data set, with out of the peptide binding cores correctly predicted. a recent study by wang et al. performed a similar analysis using an overlapping but slightly different set of peptide-mhc complexes and obtained qualitatively similar results for the same prediction methods examined in that study. that study concluded that the tepitope method, as implemented by the propred web server, achieved the best overall performance. we found similar results, with only rta and propred correctly predicting all hla-dr peptide binding registers. nn-align achieved similar performance with correctly predicted hla-dr peptide binding cores. significantly, some prediction methods, which generally attain accuracies that are comparable to or better than the earlier tepitope method in predicting peptide binding affinities, had significantly worse accuracy than tepitope in predicting the peptide binding cores. this demonstrates that a prediction model that achieves high accuracy for peptide binding affinity does not necessarily achieve corresponding high accuracy for predicting the peptide binding core. the fact that our rta method does perform well for both prediction tasks suggests that the underlying approximate physical model for the contribution of each binding register to the total binding affinity is correct and accurate. importance of suboptimal terms in the boltzmann sum a common approximation in other prediction models, such as tepitope, is to use only the peptide binding register with the highest affinity and so neglect the contributions of suboptimal binding registers. under this approximation eq. would become this approximation is accurate if only one term dominates the boltzmann average in eq. in order to estimate the importance of the suboptimal terms included in our model, we calculated the ratios and |t2/t1|, in which ti are the terms in the boltzmann sum in eq. ordered so that |t1| ≥ |t2| ≥ l ≥ |tl-8|. t are then the terms with the largest absolute values. both ratios vary from to. these ratios were calculated for the largest iedb data set, containing peptides binding to hla-drb1* a low value of indicates that suboptimal terms make a significant contribution to the sum so that the approximation using only the optimal term is inaccurate. the cumulative distribution of this quantity is shown in figure for example, suboptimal terms contribute at least 50% to the boltzmann sum for about 20% of the peptide sequences. likewise figure shows the cumulative distribution for the relative contribution of the leading suboptimal term relative to the optimal term, |t2/t1|. for example, it can be seen that this ratio is greater than for about 40% of the sequences. the ratios for other data sets are expected to have a similar distribution. to the extent that these predicted values are accurate, the distributions of these ratios indicate that suboptimal terms make significant contributions to the total binding affinity for many peptides so that the assumption that only one binding register dominates is incorrect for these peptides. the error in the predicted binding affinity for these peptides in a prediction model including only the optimal term may not be as large as expected since the model parameters are fit directly to experimental binding affinities and so the model may compensate to some degree for the contribution of missing suboptimal terms. in any case, to the extent that these estimates of the contributions of suboptimal terms are correct, prediction models, such as rta, that incorporate these terms are expected to be more accurate.  the prediction results for the rta method show that suboptimal peptide binding registers make a significant contribution to the total peptide binding affinity for many peptides, regularization with the l parameter constraint improves prediction accuracy for novel data, the method achieves good accuracy in predicting peptide-mhc binding affinities with generally higher auc values than other prediction methods evaluated on the same data sets, and the method accurately predicts the 9-residue binding core for all but one out of peptide-mhc complexes examined. while peptide-mhc binding affinities can be quickly calculated using eq. the efficient solution of the difficult global optimization problem for fitting parameters using local optimization methods contributed to the success of the method. the accuracy of the method is expected to continually improve as more experimental binding data become available, particularly for certain mhc allotypes with limited current data. one direction for future work is to extend the rta model to a multi-allotype model by including another parameter dimension reflecting the mhc peptide binding determinants, either through unique binding pockets, as in tepitope, or through polymorphic pocket residues, as in netmhciipan. another possible extension is to include interaction terms for nearby peptide residues. the approximation that contributions of individual peptide residues to the total binding free energy are largely independent is supported by side chain scanning experiments and by the accuracy of prediction models, such as rta, that are based on this approximation. however, recent experimental studies have found evidence for cooperative effects between peptide residues so that additional pairwise interaction terms may improve prediction accuracy. both of these modifications will significantly increase the number of model parameters, making regularization even more crucial for prediction accuracy. authors' contributions ajb conceived of the study, collected the data sets, and analyzed the prediction results. hdm performed the numerical optimization to determine the model parameters. both authors participated in drafting the manuscript and read and approved the final manuscript. supplementary material additional file 1 detailed description of the algorithm for creating cross-validation data sets with minimal peptide sequence overlap. click here for file additional file 2 this zip archive contains tab-separated files with the 5-fold cross-validation sets for different mhc allotypes created for this study and also a separate text file describing the table format. click here for file additional file 3 tables s1-s present the prediction results for the iedb benchmark sets from el-manzawaly et al.. prediction accuracies are given for rta as well as the three prediction methods studied in that paper. click here for file acknowledgements the work of ajb was supported by the mayo clinic.
we have developed the regularized thermodynamic average, or rta, method for predicting the affinities of peptides binding to class ii mhc. rta accounts for all possible peptide binding conformations using a thermodynamic average and includes a parameter constraint for regularization to improve accuracy on novel data. finally, we found that suboptimal peptide binding registers, which are often ignored in other prediction methods, made significant contributions of at least 50% of the total binding energy for approximately 20% of the peptides. the rta method accurately predicts peptide binding affinities to class ii mhc and accounts for multiple peptide binding registers while reducing overfitting through regularization. each mhc allotype generally binds a distinct subset of peptides and the enormous number of possible peptide epitopes prevents their complete experimental characterization. rta also gave the highest accuracy on all but three allotypes when compared with results from different prediction methods applied to the same data. rta was shown to achieve higher accuracy, as measured by auc, than smm-align on the same data for all mhc allotypes examined. a web server implementing the rta prediction method is available at http://bordnerlab.org/rta/. computational methods can utilize the limited experimental data to predict the binding affinities of peptides to class ii mhc. the binding of peptide fragments of extracellular peptides to class ii mhc is a crucial event in the adaptive immune response.
 the binding of peptide fragments of extracellular peptides to class ii mhc is a crucial event in the adaptive immune response. each mhc allotype generally binds a distinct subset of peptides and the enormous number of possible peptide epitopes prevents their complete experimental characterization. computational methods can utilize the limited experimental data to predict the binding affinities of peptides to class ii mhc.  we have developed the regularized thermodynamic average, or rta, method for predicting the affinities of peptides binding to class ii mhc. rta accounts for all possible peptide binding conformations using a thermodynamic average and includes a parameter constraint for regularization to improve accuracy on novel data. rta was shown to achieve higher accuracy, as measured by auc, than smm-align on the same data for all mhc allotypes examined. rta also gave the highest accuracy on all but three allotypes when compared with results from different prediction methods applied to the same data. in addition, the method correctly predicted the peptide binding register of out of peptide-mhc complexes. finally, we found that suboptimal peptide binding registers, which are often ignored in other prediction methods, made significant contributions of at least 50% of the total binding energy for approximately 20% of the peptides.  the rta method accurately predicts peptide binding affinities to class ii mhc and accounts for multiple peptide binding registers while reducing overfitting through regularization. the method has potential applications in vaccine design and in understanding autoimmune disorders. a web server implementing the rta prediction method is available at http://bordnerlab.org/rta/. 