Skip to main content
  • Research article
  • Open access
  • Published:

Statistical deconvolution of enthalpic energetic contributions to MHC-peptide binding affinity



MHC Class I molecules present antigenic peptides to cytotoxic T cells, which forms an integral part of the adaptive immune response. Peptides are bound within a groove formed by the MHC heavy chain. Previous approaches to MHC Class I-peptide binding prediction have largely concentrated on the peptide anchor residues located at the P2 and C-terminus positions.


A large dataset comprising MHC-peptide structural complexes was created by re-modelling pre-determined x-ray crystallographic structures. Static energetic analysis, following energy minimisation, was performed on the dataset in order to characterise interactions between bound peptides and the MHC Class I molecule, partitioning the interactions within the groove into van der Waals, electrostatic and total non-bonded energy contributions.


The QSAR techniques of Genetic Function Approximation (GFA) and Genetic Partial Least Squares (G/PLS) algorithms were used to identify key interactions between the two molecules by comparing the calculated energy values with experimentally-determined BL50 data. Although the peptide termini binding interactions help ensure the stability of the MHC Class I-peptide complex, the central region of the peptide is also important in defining the specificity of the interaction. As thermodynamic studies indicate that peptide association and dissociation may be driven entropically, it may be necessary to incorporate entropic contributions into future calculations.


Major Histocompatibility Complex (MHC) molecules play a central role in the adaptive immune system, forming a complex with foreign antigenic peptides and displaying them on the cell surface for inspection by receptors of cytotoxic and helper T cells. In consort with the T cell receptor, the repertoire of MHC-peptide complexes shapes the specificity of the T cell response. Two different types of MHC molecules – Class I and Class II – are recognized by distinct sets of T-cells: CD8+ and CD4+ respectively. MHC Class I molecules, which predominantly present antigenic material derived from the cytosol, are composed of an α heavy chain, a light chain (β2-microglobulin or β2m) and a peptide between 8–15 amino acids in length. The peptide is bound within a groove formed by the MHC heavy chain.

Analysis of MHC Class I-peptide binding interactions together with specific peptide-MHC crystal structures reveals several common trends. A conserved cluster of tyrosine residues are present at the N terminus of the peptide while a different set of MHC Class I residues form a series of hydrogen bonds and ionic interactions with the C terminus and peptide backbone. Previous attempts to predict potential epitopes have concentrated on binding motifs. Such motifs greatly emphasize so-called 'anchor' residues at peptide position P2 and the peptide C-terminus, which occupy the B and F pocket of the MHC Class I binding groove respectively [1, 2]. The B pocket is composed of eight α heavy chain residues (9, 24, 45, 66, 67, 70 and 99) and the F pocket is composed of four (77, 80, 81 and 116). Previous analysis of HLA-A*0201 [3] found that residues such as leucine, isoleucine, valine and methionine are favoured for the P2 position while the optimal residues for the C-terminus are leucine and valine (although methionine, isoleucine, alanine and threonine can also occupy this position). In some cases, contributions made by the peptide residues at positions 1, 3, 6 and 7 can compensate for the presence of unfavoured residues at the P2 and C-terminal positions). In an attempt to discover novel T cell epitopes, many motif-based techniques have been used to predict the affinity of MHC-peptide interactions. Diverse methods, such as EpiVax [4], Artificial Neural Networks [5], Hidden Markov Models [6], Support Vector Machines [7] and Profiles [8], have all been used to create T-cell epitope prediction programs. Typically, prediction techniques have been calibrated by using experimentally-determined IC50 and BL50 data [9], at least indirectly. Purely motif-based analysis of binding, however, may be too simplistic to accurately detect the effect of non-anchor residues within the groove. This is because changes in the peptide sequence may introduce subtle structural variations into the binding which can greatly influence peptide binding affinity [10]. These interactions can be partially accounted for by the additive effect of the binding. The calculation of the additive effect is based upon the method first introduced by Free and Wilson whereby each substitute makes additive and constant contributions to biological activity regardless of substituent variation. This assumes that peptide binding affinity can be decomposed into the sum of contributions from residues at each position and the interaction between those residues. The QSAR (Quantitative Structure-Activity Relationship)-based additive technique [1016] is based on a more sophisticated approach than motifs and uses additive technique in the calculation. The additive QSAR model therefore accounts for the sum of amino acid interactions, the sum of adjacent peptide side chain interactions, and the sum of every second side chain interaction, and so on.

A possible way to full incorporate the additive effects into the prediction is to derive the affinities based on structural data. The COMBINE (COMparative BINding Energy) method does this by creating a model which estimates binding free energy differences based on energy-minimised structures of ligand-receptor complexes [1719]. The model is fitted to a set of experimental binding parameters and is used to predict the affinity of novel ligands. This technique assumes that the free energy of binding correlates with a subset of energy components determined from the structure. In this paper, a technique similar to COMBINE is developed which incorporates both structural and data-driven methods in order to create a novel method for the analysis of peptide-MHC binding. This new method seeks to better understand the atomic interactions that determine binding affinities: static energy analysis calculations combined with QSAR analysis are used to calculate the relative contribution of interactions between the MHC Class I allele A*0201 and bound peptide residues, thus generating an accurate and predictive system for the A2 supertype.


QSAR analysis was carried out on three individual models: total/electrostatic/van der Waals energy descriptors; electrostatic/van der Waals energy descriptors and total non-bonded energy descriptors. The matrix consisted of 118 peptide rows and 4860 descriptor (energy, electrostatic and van der Waals) columns. To make this analysis tractable, any columns consisting entirely of values of ~0 kcal mol-1 (signifying that no contributions were made between the peptide and binding site) were omitted from the dataset, yielding a final matrix of 1504 descriptor columns (454 total energy descriptors, 331 electrostatic energy descriptors and 719 van der Waals energy descriptors). The resulting regression equation modeled the binding affinities of a training set of 83 peptides and was used to predict the binding affinities of an independent test set of 35 peptides.

The GFA and G/PLS algorithms were used to select the optimum number of descriptors for use in the regression analysis. Although the methods are similar, the application of G/PLS often allows the construction of larger QSAR equations while still avoiding over-fitting and eliminating most variables. Within the Cerius2software, the crucial variable is the initial number of descriptors chosen in the regression equation although this can be varied during the process by including a mutation probability for increasing or decreasing the equation length. The statistical indices for the best GFA and G/PLS algorithms for all 3 models are shown in Table 1. For all the QSAR models, the number of terms in the initial and final equations was varied between 15–20 descriptors. It was found that raising or lowering the number of descriptors outside of this range caused a steady increase in the PRESS value and a lowering of the r2. Models were created with 10, 25 and 30 descriptors and all exhibited a decrease in the quality of the model. The GFA Total + Electrostatic + VDW model is shown in Table 2 as example. It is likely that the use of too many descriptors or too few will cause overfitting or underfitting of the model respectively. 15–20 descriptors may represent a range where there is sufficient information to build a model without overemphasising specific interactions that do not represent general trends within the dataset. The best model was chosen based on both the r2 value of the test set and also on the (CV)r2 value from the training set. All selected models had the highest available test set r2 value, except for the G/PLS Electrostatic + VDW model where the 15 descriptor model was chosen as it has a (CV)r2 value significantly higher than the 16 descriptor model. Marked equations had a few anomalous outliers removed to consolidate the model.

Table 1 QSAR Results Based on full binding groove for GFA and G/PLS Total, Electrostatic + VDW and Total + Electrostatic + VDW models.
Table 2 QSAR Results for GFA Total + Electrostatic + VDW model based on 10, 20 25, 30 descriptors.

Further analysis was also carried out limiting the potential descriptors to those amino acids comprising previously-defined MHC anchor pockets: P2 and the C terminus [1]. In that model, a total of 138 residues were used (46 energy terms, 46 electrostatic terms and 46 van der Waals terms). Despite the reduced population size, the optimal number of descriptors still occurred within the 10–20 descriptor range (see Table 3). The reduced number of descriptors generally seems to decrease the predictivity of the models significantly. However, one exception (the G/PLS Total + Electrostatic + VDW model) did produce a test set r2 value of 0.711, which is higher than any model produced using the full range of descriptors. The GFA and G/PLS Total models produced extremely poor results with a test set r2 value of 0.044 and 0.263 respectively. The best individual models for both the GFA and G/PLS algorithms taken from both Tables 1 and 3 are presented in Table 4 and 5 respectively. Descriptor sets taken from the best models generated by the GFA and G/PLS calculations were used to validate the QSAR analysis using the PLS module in SYBYL 6.9. The predictivity of GFA and G/PLS models compares well with the PLS method, indicating that the predictive ability of the models is consistent in terms of descriptor choice, despite differences in the form of the predictive equation in each case.

Table 3 QSAR Results Based on Binding Pockets for GFA and G/PLS Total, Electrostatic + VDW and Total + Electrostatic + VDW models.
Table 4 Best QSAR models for GFA and G/PLS based on full binding groove in comparison with the PLS method. *PLS Models based on Cerius GFA & G/PLS Models respectively.
Table 5 Best QSAR models for GFA and G/PLS based on binding Pockets in comparison with the PLS method. * PLS Models based on Cerius GFA & G/PLS Models respectively.

Table 6 shows the residues that make up the descriptors for each of the six Cerius2QSAR models presented in Table 3. Each model's residue distributions around the peptide are shown in Figures 1, 2, 3, 4, 5, 6. The selected coefficients are generated from those residues which make a consistent contribution over the entire training set. An overview of the QSAR models would suggest that although specific MHC-peptide interactions are difficult to identify, certain regions of the α heavy chain seem to relate to specific peptide residues. P2 is associated with the α60–70 region while the α160–170 region is associated with P2-P4. The α150–160 and α140–145 regions relate to P5 and P7 specifically. These associations may be fundamental to the binding and are maintained irrespective of the sequence of the peptide. In several cases, the MHC residues interact with the main chain atoms of the peptide and therefore make energetic contributions independent of the peptide sequence. Conversely, regions α1–4, α26–33, α46–58, α101–111 and α168–180 all make little or no energetic contribution due to their distance from the peptide residues. It is difficult to identify specific interactions from MHC residues, however, the following associations between specific MHC and peptide residues are all represented in three out of the six models: Met 5 interacting with P3, Asp 77 interacting with P9, Thr 143 interacting with P9, Val 152 interacting with P6 and Leu 156 interacting with P5. This might suggest that those five interactions all have an effect in determining the specificity of the interaction. Of the B pocket residues, Tyr 99 is the only residue to appear in two models while B pocket residues Ala 24, Met 45, Lys 66 and Tyr 116 only appear in one model each. Asp 77 appears in three models as is the only one of the four residues lining the F pocket which is used as a descriptor. Both Lys 66 and Asp 77 have been identified as forming conserved hydrogen bonds within the groove, respectively with the oxygen of the P2 residue and the nitrogen of the P9 residue [1]. Lys 66 is of particular interest as it has been identified as being a functional "hotspot" due to the formation of an ion pair with Glu 63 and an apparent dual role in both peptide binding and T cell receptor interaction [2022]. In spite of this the residue is only used as a descriptor by G/PLS Total + Electrostatic + VDW model. The conserved bonds, Glu 63- P2 amide nitrogen and Tyr 99-P3 amide nitrogen, are not represented at all. Both residues, however, are used as descriptors in relation to P1 and in relation to P2 and P4 respectively in separate models. Previous work [10] has identified residues 9, 97, 114 and 116 as key to defining the A2-supertype, hence representative of vital points of interaction within the groove. They are poorly represented by the QSAR models, the only equation that incorporates one of them is the G/PLS Electrostatic + VDW model, where 116 interacts with P5. However, they do all occur within regions of the groove identified as important to the binding, again suggesting that it is regions rather than specific residues that determine the stability of the complex.

Table 6 Residue descriptors for the Best QSAR Results based on full binding groove.
Figure 1
figure 1

GFA Total + Electrostatic + VDW model. The bound peptide (red) is surrounded by the descriptor peptides (bubble opaque). Despite the apparent distance in some instances, it is established that long range interactions can be as important to the specificity of the binding as short range ones.

Figure 2
figure 2

G/PLS Total + Electrostatic + VDW model. The bound peptide (red) is surrounded by the descriptor peptides (bubble opaque). Despite the apparent distance in some instances, it is established that long range interactions can be as important to the specificity of the binding as short range ones.

Figure 3
figure 3

GFA Electrostatic + VDW model. The bound peptide (red) is surrounded by the descriptor peptides (bubble opaque). Despite the apparent distance in some instances, it is established that long range interactions can be as important to the specificity of the binding as short range ones.

Figure 4
figure 4

G/PLS Electrostatic + VDW model. The bound peptide (red) is surrounded by the descriptor peptides (bubble opaque). Despite the apparent distance in some instances, it is established that long range interactions can be as important to the specificity of the binding as short range ones.

Figure 5
figure 5

GFA Total model. The bound peptide (red) is surrounded by the descriptor peptides (bubble opaque). Despite the apparent distance in some instances, it is established that long range interactions can be as important to the specificity of the binding as short range ones.

Figure 6
figure 6

G/PLS Total model. The bound peptide (red) is surrounded by the descriptor peptides (bubble opaque). Despite the apparent distance in some instances, it is established that long range interactions can be as important to the specificity of the binding as short range ones.


Overall, within our models, the distribution of the descriptors does not lend support to the view that burial of anchor residues, within defined cavities, are mainly responsible for stabilising the peptide within the groove. Rather it indicates a fairly even distribution of MHC residue interactions along the length of the peptide. Attempts to improve the QSAR model by focusing on the MHC residues associated with the anchor residues were largely unsuccessful (see Table 3 and Table 5), in spite of the predictive accuracy of the G/PLS Total + Electrostatics + VDW model. Despite the fact that the anchor models generate a similar number of descriptors to those of the full groove, focusing the interaction towards the peptide termini at the expense of the variable region generally reduces the algorithm's predictivity. It is possible that this reflects the comparative dominance of electrostatic interactions in determining binding affinity within both of the anchor pockets. Model predictivity also deteriorates after removal of each descriptor type compared to the full residue model, particularly in the case of models based on Total energy. Moreover, attempts to split the interaction energy into long range electrostatic and short range van der Waals forces also decreases correlation in the full residue model, suggesting that a single residue's total energy contribution is primarily enthalpic in character. It should be noted that several of the residues highlighted by the models are distant from the peptide (see Figures 1, 2, 3, 4, 5, 6), demonstrating the importance of long-range electrostatic interactions, which can be up to 8Å length, in stabilizing the complex. This could also be an effect of superadditivity correlation (SAC) [23] where correlation between variables between two sites does not imply that the two sites directly interact. Mutational data has indicated that long range non-causal interactions can be stronger in magnitude than causal correlations induced by direct interactions [24, 25].

The Total energy models showed the highest correlation in the test set for both the GFA and G/PLS model although the (CV)r2 value of the G/PLS model dropped sharply in the training set. Both models placed great emphasis on the conserved P1-P3 region of the MHC Class I groove at the expense of the central and C-termini regions. The key interaction of Asp 77 with the P9 residue was incorporated into both the GFA and G/PLS model, indicating that it may be a vital descriptor for a successful predictive algorithm. Both sets of models where the energetic contributions were broken down into the electrostatic and van der Waals contributions showed a lower degree of correlation (although the PRESS values for the Test Set compared favourably with the Total Energy Models). Building comparative models using the program SYBYL (see Table 4 and 5) produced a consistently lower level of correlation compared to the Cerius2model but the apparent trends are comparable. The QSAR models also show a similar degree of correlation to models generated by the additive method [10] which produced a higher r2value of 0.954 but an equivalent (CV)r2 value of 0.602 from a model containing only 6 components.

The limitations of motif-based methods and the potential of structural analysis may be in part due to the variable nature of the binding with the central region of the peptide. Where as the termini regions (P1-P3 and P8-P9) are heavily constrained by interactions with conserved MHC residues, the central region (P4-P7) typically exhibits a far more variable conformation. Analysis of the MHC Class I structure [1] has shown that the central region is less constrained within the groove than the termini and can dramatically alter its conformation in response to the peptide sequence. It is also the region which is most exposed to the T cell receptor and because of this may arch away from the floor of the groove. The conformational variability of the central region means that definite contacts with specific MHC residues are not obvious. This makes predicting the contribution to affinity of that region of the ligand for the MHC molecule difficult. Nevertheless, the models demonstrate that the central region does make an important contribution to the specificity of the interaction.

It is also possible that static energy analysis overlooks significant energetic interactions occurring during binding. The thermodynamic property, which we are trying to estimate with our scoring function, is the energy released when ligand and receptor bind; the so-called Gibbs free energy of binding: ΔGbind. It can be represented as ΔGbind = ΔH - TΔS where ΔH is the enthalpic (internal) energy and TΔS is an entropy term which is indicative of the relative gain or loss of disorder upon binding [26]. The binding groove is hydrophobic in character, particularly in the region of the anchor residues, and the calculated energy interactions reflect this, particularly with the high incidence of charge-charge interactions within the groove. It is necessary for the solvent entropy from the burial of hydrophobic groups to offset the reduction in peptide conformational entropy that occurs upon binding. The measured interactions therefore reflect the hydrogen bonds and salt bridges that are formed when the peptide moves from the loosely packed, partially bound state to a being fully bound into a complex. Water is stripped from partial and full charges on the MHC and peptide residues and there is a reduction in the favourable hydrogen bond enthalpy associated from interactions with bound water molecules. In taking our work forward, it is therefore necessary to incorporate the entropic contribution to the energetic calculation into the descriptors in order to create a more accurate representation of the free energy change between the MHC molecules bound and unbound states. This may be done by use of the computationally efficient MM_PBSA/GBSA method which combines the molecular mechanical energies with the continuum solvent approaches [27, 28]. It is hoped that the inclusion of these entropic terms may significantly increase the predictive capability of future QSAR analysis. However, the local interactions have always proved to be more relevant to determining the affinity of a given peptide than the general properties of the ligand-receptor interaction [29].

Previous work [1016] has already indicated that motif-based analysis of the anchor residues within the MHC Class I binding groove can only provide limited accuracy in determining potential T-cell epitopes. Thus the QSAR models presented here strongly challenge the assumption that only the P2 and P9 residues of a nonameric peptide make a significant contribution to binding. While the majority of the descriptors do concentrate on the termini regions of the peptide, there is a clear, and substantial, energetic contribution made by the structurally variable central region. It is possible that more detailed free energy calculations of the MHC Class I-peptide binding interactions, whereby the entropic contribution to the association of the peptide to the MHC Class I molecule is incorporated into the model, will, in the future, lead to the generation of more accurate predictive algorithms.


Energy minimization of HLA-A*0201 structure

An HLA-A*0201 dataset, consisting of 118 peptides of known sequence, was obtained from Doytchinova et al. [10]. The structure of the MHC Class I allele HLA-A*0201 bound to a nonameric melanocyte-melanoma tumor-antigen peptide elucidated by X-ray diffraction at a resolution of 2.15Å [2] (pdb code: 1JHT) was used as the template. Residues of the bound peptide were remodeled using the crystallographic modeling program 'O' [29] to generate 118 MHC Class I-peptide structures equivalent to the data set. Energy minimization was carried out on the structures using AMBER Version 8 [30]. The interatomic potential energy function uses the AMBER forcefield 03 for parameterisation. Hydrogen atoms were added to the structure, and the system was fully solvated using TIP3 waters [31]. This function was performed by the LEaP program [30]. The full structure of the HLA-A*0201-peptide complex was explicitly represented within the simulation. The energy of the solvated molecular complex was minimized using a steepest descent method that continued for 20,000 1 fs time steps or until the root mean square deviation between successive time steps had fallen below 0.01Å. The peptides were then remodeled and annealed by raising the temperature of the system from 0 to 500 K for a period of 8 picoseconds and maintaining the system at that temperature for a further 6 picoseconds. The system was cooled to 0.2 K over a period of 33 picoseconds before being rested at 1 K for a further 30 picoseconds. Both the minimization and the annealing were performed using the sander program and static energy analysis was then carried out on the minimised structures using the anal program [30]. Anal calculates the group-group interaction energies between different parts of the system based upon the position of their composite atoms. The interaction energies, which are measured in kcal mol-1, between the nine peptide residues and the first 180 residues of the MHC Class I heavy chain (which incorporates the entire MHC Class I binding groove) were calculated. Three values were generated for each interaction; the electrostatic interaction energy, the van der Waals interaction energy and the total non-bonded interaction energy. Total energies are the sum of electrostatic and van der Waals terms; since these values can have different signs and be significantly different in magnitude, their sum gives rise to a non-trivial and informative descriptor.

QSAR regression

The BL50 binding affinities for the dataset, were extracted from the AntiJen database [3234]. The half-maximal binding level (BL50) is the peptide concentration yielding the half-maximal Fluorescence Index (FI) of the reference peptide in each assay. These values were converted to pBL50 (-logBL50) and used as the dependent variables in the QSAR regression. Peptides were partitioned between training and test datasets in a ratio of 70:30. A program was developed which transformed the nine amino acid peptide sequence into rows of a table. When choosing a test set, the program omits rows that contain the least contributing amino acid starting at position 1 and working systematically to positions 2, 3, through to position 9 until 30% of the dataset remains. Three different QSAR techniques were employed during the study: the Genetic Function Approximation (GFA) algorithm [35], Genetic Partial Least Squares (G/PLS)[36], which were both carried out using Cerius2[37], and Partial Least Squares (PLS), carried out using SYBYL 6.9 [37]. The genetic algorithm first selects equations that contain a randomly chosen subset of dependent variables (or descriptors) and, subsequently, improved equations are 'evolved' using a genetic crossover operation, which selects the best descriptors using multiple linear regression, and finally uses a least squares technique to produce the final model. G/PLS, an alternative to the GFA algorithm, uses GFA to select appropriate basis functions and PLS regression as the fitting technique to weight the basis functions' relative contributions in the final model. In the case of both algorithms, the population size of the initial equation set was fixed at 100 while the number of generations of equations (genetic crossovers) was set at 5000. Once a QSAR calculation has been completed, a list of outliers is generated by the GFA and G/PLS algorithms. Outliers are removed from the training set, the QSAR is repeated again and new equations are generated, until, at most, 10–15% of the training set outliers are removed.

Within Cerius2the validation terms used were correlation coefficient (r2), Friedman's Lack of Fit (LOF)[38] and the Cross-Validation correlation coefficient ((CV)r2) using the Leave-One-Out (LOO) procedure [39]. Adding more terms to the regression model always reduces the normal least squares error (LSE) but does not necessarily reduce the LOF measure. Thus, as a new term reduces the LSE, it also tends to increase the value of LOF. Within SYBYL, the predictive power of the PLS Cross-Validation Leave-One-Out (CV-LOO) model is assessed using the following parameters: cross-validated coefficient (q2) and the Standard Error of Prediction (SEP) [1016]. The optimal number of components (NC) resulting from the CV-LOO is then used in the non-cross validated model which was assessed using standard Multiple Linear Regression (MLR) validation terms: explained variance (r2) and Standard Error of Estimate (SEE) [1016]. Increasing the number of components improves the fit between target and explanatory properties and corresponds to the best q2.



Genetic Function Approximation


Genetic Partial Least Squares


Major Histocompatibility Complex


COMparative BINding Energy


Superadditivity Correlation


Partial Least Squares


Lack of Fit


Cross-Validation correlation coefficient




Standard Error of Prediction


Number of Components


Standard Error of Estimate


Multiple Linear Regression


  1. Saper MA, Bjorkman PJ, Wiley DC: Refined structure of the human histocompatibility antigen HLA-A2 at 2.6Å resolution. J Mol Biol 1991, 219: 277–319. 10.1016/0022-2836(91)90567-P

    Article  CAS  PubMed  Google Scholar 

  2. Sliz P, Michielin O, Cerottini JC, Luescher I, Romero P, Karplus M, Wiley DC: Crystal structures of two closely related but antigenically distinct HLA-A2/melanocyte-melanoma tumor-antigen peptide complexes. J Immunol 2001, 167: 3276–3284.

    Article  CAS  PubMed  Google Scholar 

  3. Sidney J, Southwood S, Mann DL, Fernandez-Vina MA, Newman MJ, Sette A: Majority of peptides binding HLA-A*0201 with high affinity crossreact with other A2-supertype molecules. Hum Immunol 2001, 62: 1200–1216. 10.1016/S0198-8859(01)00319-6

    Article  CAS  PubMed  Google Scholar 

  4. De Groot AS, Sbai H, Aubin CS, McMurry J, Martin W: Immunoinformatics: mining genomes for vaccine components. Immunol Cell Biol 2002, 80: 255–269. 10.1046/j.1440-1711.2002.01092.x

    Article  CAS  PubMed  Google Scholar 

  5. Brusic V, Rudy G, Honeyman MC, Hammer J, Harrison LC: Prediction of MHC class-II binding peptides using an evolutionary algorithm and artificial neural network. Bioinformatics 1998, 14: 121–130. 10.1093/bioinformatics/14.2.121

    Article  CAS  PubMed  Google Scholar 

  6. Udaka K, Mamitsuka H, Nakaseko Y, Abe N: Prediction of MHC class I binding peptides by a query learning algorithm based on hidden Markov models. J Biol Phys 2002, 28: 183–194. 10.1023/A:1019931731519

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  7. Donnes P, Elofsson A: Prediction of MHC class I binding peptides, using SVMHC. BMC Bioinformatics 2002, 3: 25. 10.1186/1471-2105-3-25

    Article  PubMed Central  PubMed  Google Scholar 

  8. Reche PA, Glutting JP, Reinherz EL: Prediction of MHC class I binding peptides using profile motifs. Hum Immunol 2002, 63: 701–709. 10.1016/S0198-8859(02)00432-9

    Article  CAS  PubMed  Google Scholar 

  9. Sette A, Sidney J: Nine major HLA class I supertypes account for the vast preponderance of HLA-A and -B polymorphism. Immunogenetics 1999, 50: 201–212. 10.1007/s002510050594

    Article  CAS  PubMed  Google Scholar 

  10. Doytchinova IA, Walshe VA, Jones NA, Gloster SE, Borrow P, Flower DR: Coupling in silico and in vitro analysis of peptide-MHC binding: a bioinformatic approach enabling prediction of superbinding peptides and anchorless epitopes. J Immunol 2004, 172: 7495–7502.

    Article  CAS  PubMed  Google Scholar 

  11. Doytchinova IA, Guan P, Flower DR: Coupling identifiying human MHC supertypes using bioinformatic methods. J Immunol 2004, 172: 4314–4323.

    Article  CAS  PubMed  Google Scholar 

  12. Hattotuwagama CK, Toseland CP, Guan P, Taylor PJ, Hemsley SL, Doytchinova IA, Flower DR: Class II Mouse Major Histocompatibility Complex Peptide Binding Affinity: In Silico bioinformatic prediction using robust multivariate statistics. Bioinformatics, in press.

  13. Hattotuwagama CK, Doytchinova IA, Flower DR: In Silico prediction of peptide binding affinity to class I mouse major histocompatibility complexes: A Comparative Molecular Similarity Index Analysis (CoMSIA) study. J Chem Inf Mod, in press.

  14. Hattotuwagama CK, Guan P, Doytchinova IA, Flower DR: In Silico QSAR-Based Predictions of Class I and Class II MHC Epitopes. Immunoinformatics: Opportunities and Challenges of Bridging Immunology with Computer and Information Sciences. , in press.

  15. Hattotuwagama CK, Guan P, Doytchinova IA, Flower DR: New Horizons in Mouse Immunoinformatics: Reliable In Silico Prediction of Mouse Class I Histocompatibility Major Complex Peptide Binding Affinity. Org Biomolec Chem 2004, 2: 3274–3283. 10.1039/b409656h

    Article  CAS  Google Scholar 

  16. Hattotuwagama CK, Guan P, Doytchinova IA, Zygouri C, Flower DR: Quantitative online prediction of peptide binding to the major histocompatibility complex. J Mol Graph Model 2004, 22: 195–207. 10.1016/S1093-3263(03)00160-8

    Article  CAS  PubMed  Google Scholar 

  17. Wang R, Wade R: Comparative binding energy (COMBINE) analysis of OppA-peptide complexes to relate structure to binding thermodynamics. J Med Chem 2002, 45: 4828–4837. 10.1021/jm020900l

    Article  CAS  PubMed  Google Scholar 

  18. Wang R, Wade R: Comparative binding energy (COMBINE) analysis of influenza neuraminidase-inhibitor complexes. J Med Chem 2002, 6: 961–971.

    Google Scholar 

  19. Tokarski JS, Hopfinger AJ: Prediction of ligand-receptor binding thermodynamics by free energy force field (FEFF) 3D-QSAR analysis: application to a set of peptidometic renin inhibitors. J Chem Inf Comput Sci 1997, 37: 792–811. 10.1021/ci970006g

    Article  CAS  PubMed  Google Scholar 

  20. Madden DR, Garboczi DN, Wiley DC: The antigenic identity of peptide-MHC complexes: a comparison of the conformations of five viral peptides presented by HLA-A2. Cell 1993, 75: 693–708. 10.1016/0092-8674(93)90490-H

    Article  CAS  PubMed  Google Scholar 

  21. Baker BM, Turner RV, Gagnon SJ, Wiley DC, Biddison WE: Identification of a crucial energetic footprint on the alpha1 helix of human histocompatibility leukocyte antigen (HLA)-A2 that provides functional interactions for recognition by tax peptide/HLA-A2-specific T cell receptors. J Exp Med 2001, 193: 551–562. 10.1084/jem.193.5.551

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  22. Baxter TK, Gagnon SJ, Davis-Harrison RL, Beck JC, Binz AK, Turner RV, Biddison WE, Baker BM: Strategic mutations in the class I major histocompatibility complex HLA-A2 independently affect both peptide binding and T cell receptor recognition. J Biol Chem 2004, 279: 29175–29184. Giraud BG, Heumann JM, Lapedes AS: Superadditive correlation. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 1999, 59:4983–4991 10.1074/jbc.M403372200

    Article  CAS  PubMed  Google Scholar 

  23. Oelschlaeger P, Mayo SL, Pleiss J: Impact of remote mutations on metallo-{beta}-lactamase substrate specificity: Implications for the evolution of antibiotic resistance. Protein Sci 2005, 14: 765–774. 10.1110/ps.041093405

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Marvin JS, Hellinga HW: Manipulation of ligand binding affinity by exploitation of conformational coupling. Nat Struct Biol 2001, 8: 795–798. 10.1038/nsb0901-795

    Article  CAS  PubMed  Google Scholar 

  25. Binz AK, Rodriguez RC, Biddison WE, Baker BM: Thermodynamic and kinetic analysis of a peptide-class I MHC interaction highlights the noncovalent nature and conformational dynamics of the class I heterotrimer. Biochemistry 2003, 42: 4954–4961. 10.1021/bi034077m

    Article  CAS  PubMed  Google Scholar 

  26. Honig B, Nicholls A: Classical electrostatics in biology and chemistry. Science 1995, 268: 1144–1149.

    Article  CAS  PubMed  Google Scholar 

  27. Sitkoff D, Sharp KA, Honig B: Accurate calculation of hydration free energies using macroscopic solvent models. J Phys Chem 1994, 98: 1978–1988. 10.1021/j100058a043

    Article  CAS  Google Scholar 

  28. Doytchinova IA, Walshe VA, Borrow P, Flower DR: Towards the chemometric dissection of peptide – HLA-A*0201 binding affinity: comparison of local and global QSAR models. J Comput Aided Mol Des 2005, 19: 203–212. 10.1007/s10822-005-3993-x

    Article  CAS  PubMed  Google Scholar 

  29. Kleywegt GJ, Jones TA: Model-building and refinement practice. Methods Enzymol 1997, 277: 208–230.

    Article  CAS  PubMed  Google Scholar 

  30. Case DA, Pearlman DA, Caldwell JW, Cheatham III TE, Ross WS, Simmerling C, Darden T, Merz KM, Stanton RV, Cheng A, Vincent JJ, Crowley M, Ferguson DM, Radmer R, Seibel GL, Singh UC, Weiner P, Kollman PA: AMBER 6. Univ of California, San Francisco;

  31. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML: Comparison of simple potential functions for simulating liquid water. J Chem Phys 1983, 79: 926–935. 10.1063/1.445869

    Article  CAS  Google Scholar 

  32. Blythe MJ, Doytchinova IA, Flower DR: JenPep: a database of quantitative functional peptide data for immunology. Bioinformatics 2002, 18: 434–439. 10.1093/bioinformatics/18.3.434

    Article  CAS  PubMed  Google Scholar 

  33. Sette A, Sidney J, del Guercio M-F, Southwood S, Ruppert J, Dalberg C, Grey HM, Kubo RT: Peptide binding to the most frequent HLA-A class I alleles measured by quantitative molecular binding assays. Mol Immunol 1994, 31: 813–822. 10.1016/0161-5890(94)90019-1

    Article  CAS  PubMed  Google Scholar 

  34. McSparron H, Blythe MJ, Zygouri C, Doytchinova IA, Flower DR: JenPep: A novel computational information resource for immunology and vaccinology. J Chem Inf Comput Sci 2003, 43: 1276–1287. 10.1021/ci030461e

    Article  CAS  PubMed  Google Scholar 

  35. Cerius2 © Molecular Simulations/Biosym inc. San Diego, USA; 1998.

  36. Rogers D, Hopfinger AJ: Application of Genetic Function Approximation to Quantitative Structure-Activity Relationships and Quantitative Structure-Property Relationships. Chem Inf Comput Sci 1994, 34: 854–866. 10.1021/ci00020a020

    Article  CAS  Google Scholar 

  37. SYBYL 6.7 Tripos Inc., 1699 Hanley Road, St. Louis, MO 63144; 1988.

  38. Friedman J: Multivariate Adaptive Regression. In Spline Technical Report No. 102. Stanford University, Stanford, CA; 1988.

    Google Scholar 

  39. Cramer RD, Patterson DE, Bunce DE: Crossvalidation, Bootstrapping, and Partial Least Squares Compared with Multiple Regression in Conventional QSAR Studies. Quant Struct-Act Relat 1988, 7: 18–25.

    Article  Google Scholar 

Download references


The authors would like to thank Dr. Mark T. Storr and Roger Robinson at the University of Reading for general help, advice, time and use of computer facilities in Computational Group as well as Dr Andy Purkiss and Dr David Houldershaw of Birkbeck Crystallography Department for their assistance with the energy minimisation simulations.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Matthew N Davies.

Additional information

Authors' contributions

Matthew Davies remodelled the MHC Class I structures, carried out the energy minimisation simulations and calculated the static energy interactions between the MHC molecule and its bound peptide

Channa Hattotuwagama carried out the QSAR analysis of the static energy interactions data using of Genetic Function Approximation (GFA) and Genetic Partial Least Squares (G/PLS) algorithms

David Moss supervised the molecular modeling work at Birkbeck College

Michael Drew supervised the QSAR analysis at the University of Reading

Darren Flower co-ordinated the analysis of the static energy interactions by QSAR methods and helped to interpret the data in an immunological context as well as overseeing the whole project

Matthew N Davies, Channa K Hattotuwagama contributed equally to this work.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Davies, M.N., Hattotuwagama, C.K., Moss, D.S. et al. Statistical deconvolution of enthalpic energetic contributions to MHC-peptide binding affinity. BMC Struct Biol 6, 5 (2006).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: