The molecular machinery by which an antigen presenting cell (APC) presents T cell epitopes for recognition by T cell receptors (TR) and subsequent activation of T cells followed by the immune response cascade is fascinating. T cell epitopes are short antigenic peptide sequences (p) that are bound to and presented by the major histocompatibility complexes (MHC) for recognition by the TR [1]. These epitopes are essential subunit peptides that are required in order to stimulate cellular immune responses, especially the adaptive immune responses. Peptide epitopes can be of endogenous (processed within the cell) or exogenous (processed outside the cell) origins, which are presented for surveillance and recognition by the TR in an MHC allele and supertype dependant manner. Broadly classified into two types, MHC class I (MHC-I) complexes bind and present endogenous peptides whereas MHC class II (MHC-II) complexes prefer exogenous peptides. Typically, MHC-I proteins are heterodimers, consisting of a heavy α chain (I-ALPHA) of about 45 kDa, and a light chain, β2-microglobulin (Β2M) of about 12 kDa [2, 3]. The α chain consists of α1 (G-ALPHA1), α2 (G-ALPHA2) and α3 (C-LIKE) domains where G-ALPHA1 and G-ALPHA2 domains form the peptide binding groove or ‘cleft’ [4]. MHC-II proteins are also heterodimeric proteins consisting of an α chain (II-APLHA; 34 kDa) and a β chain (II-BETA; 29 kDa) with very similar overall quaternary structure to that of MHC-I proteins [5, 6, 7, 8, 9, 10]. However, their peptide binding groove is formed by the α1 and β1 domains of the two chains.

Peptides presented by MHC-I are generally between 8-11 amino acids in length. These peptides are ‘chopped’ within the cytosol of the cell by cytosolic proteases and are transported to the MHC binding groove within the endoplasmic reticulum by the transporters associated with antigen processing (TAP) proteins in an ATP dependant manner. Following which, the peptides bind to the MHC to form the peptide-MHC (pMHC) complex which is then transported to the APC cell surface and presented for recognition by the TR of CD8⁺ cytotoxic T cells (CTLs). Similarly, the peptides presented by MHC-II are usually 12-25 amino acids in length and are endocytosed into the cell by the lysosomes where they bind the MHC-II proteins by displacing the original MHC-II ligand known as the ‘CLIP’ peptide to form the pMHC complex. And again, they are transported to the APC cell surface for recognition by the TR of the CD4⁺ T helper cells. Identification of true T cell epitopes from the repertoires of immunologically significant antigenic peptide sequences is a vital prerequisite in the process of conventional molecular vaccine design for prevention and treatment of infectious, autoimmune, allergic and graft vs. host diseases. The key step in TR-mediated immune response is thus the binding and presentation of antigenic endogenous or exogenous peptide epitopes, which can be reasonably well predicted using sequence-based methods for alleles with large datasets of known binding peptides, as reviewed earlier [11, 12].

Experimental identification of T cell epitopes is a tedious, time consuming and expensive process owing to the large number and diversity of both MHC alleles and the antigenic peptides. Not to mention, is the extremely low chance of immunogenicity (1 in 2000 peptides) even amongst the peptides that bind strongly to the MHC (50%) [13]. Recently developed computational methods have proven to be vastly time and cost efficient in screening the vast oceans of peptides and MHC repertoires [14]. Current computational methods can be broadly classified into: 1. Sequence-based approaches which use sequence motifs [15], matrix models [16, 17], Artificial Neural Network [18, 19, 20], Hidden Markov Model [21] and Support Vector Machine [22, 23, 24] for large-scale screening of potential T cell epitopes from protein sequence databanks and 2. Structure-based approaches such as protein threading [25, 26], homology modeling [27, 28], rigid docking [29] and flexible docking [2, 3] which utilize three-dimensional data for detailed structural analysis of interactions between the MHC and bound segmental antigenic peptides. The former are more suitable for large-scale screening of potential T cell epitopes, while the latter work better for detailed analysis of short immunogenic regions of antigens [2]. Although sequence-based methods are well established, a major limitation of such techniques is the heavy reliance on the availability of large comprehensive training sets of peptides. Thus, these approaches are not appropriate for accurate prediction of peptides in circumstances where the data available is insufficient. Therefore, the coverage of sequence-based techniques is limited to subsets of binding peptides that belong to the most numerous groups and cannot generate reliable data for peptides that are least represented in the dataset [2], leaving structural immunoinformatics as the only option for such peptides [3, 5, 6, 7].

Antigenic peptides that bind strongly to MHC alleles are known to elicit T cell responses [1, 2, 3, 5, 6, 7, 11]. Hence, their identification is a vital first step in the process of structure-based immune epitope prediction. The usual approach adopted to address this important issue is to utilize a powerful concept, based on the principle of structure-based drug design called “docking”, where peptides are computationally placed in MHC grooves in the best orientation, reflecting steric and electrostatic complimentarity, using structure-based docking techniques. The accuracy with which the peptides are docked is measured in terms of Root Mean Square Deviation (RMSD) values obtained by comparing the docked conformations of the peptides to their original bound conformations in the respective X-ray crystal structures. With the development of new structural modeling and docking techniques and an increase in the number of protein structures deposited in the Protein Data Bank (PDB) [30] and the IMGT/3Dstructure-DB [31, 32], structure-based approaches are being more commonly used to predict potential T cell epitopes [33], often producing modeled structures accurate to within 2.00Å RMSD from the experimental crystal structure, providing a wealth of information for structural analysis and the development of prediction methods.

The development of an accurate protocol for flexible docking has helped us to successfully carry out quantitative predictions for both MHC-I and MHC-II alleles even with limited binding peptide data [3, 5, 6, 7]. Our earlier docking protocol consisted of three steps (extended to four for pMHC-II complexes for incorporating the flanking residues on either side of the nonameric core, which is the 9-mer anchored to the MHC molecule): (1) rigid docking of the peptide nonamer termini into the MHC binding groove; (2) loop closure of central residues by satisfaction of spatial constraints; (3) followed by iterative ab initio refinements of ligand backbone and; (4) extension of flanking peptide residues by satisfaction of spatial constraints [2, 3] (only for MHC-II related peptides). While accurate, this approach has multiple steps, resulting in suboptimal computational speeds. Therefore, the efficiency of this protocol for peptide docking to MHC needs to be improved for large-scale screening of T cell epitopes. A grid-based docking methodology has earlier been reported [34] to be highly accurate in pMHC docking over a limited MHC-I data. Hence, we have developed a grid-based peptide docking method (pDOCK) and have extensively tested it on both MHC-I and MHC-II peptides. The motivation behind the development of a faster and more accurate peptide docking methodology was to eventually improve the qualitative and quantitative efficacy of structure-based T cell epitope prediction.

In this study, we present pDOCK: a new computational technique for rapid and accurate docking of flexible peptides to the MHC receptors and primarily apply it to re-dock a non-redundant dataset of 186 (149 MHC-I and 37 MHC-II related) peptides, from MPID-T2 (http://biolinfo.org/mpid-t2) database for which X-ray crystal structures are available in the PDB and the IMGT/3Dstructure-DB, back into their respective MHC grooves. pDOCK comprises of two input preparatory steps followed by a single consolidated docking and refinement step as depicted in Figure 1. The pDOCK protocol involves: Preparatory step 1: receptor modeling and positioning; Preparatory step 2: determining the docking grid by defining the grid dimensions (length x breadth x height) for ligand placement and grid map generation within the vicinity of the receptor’s binding site and; Final docking and refinement step: ligand positioning within the grid, flexible docking of the peptide into the peptide binding groove and refinement of all ligand and binding site residues using the Internal Coordinate Mechanics (ICM) global optimization, docking algorithm [35] and a biased Monte Carlo procedure (see Methods section for more details). Our preliminary analysis of all pMHC complexes from the MPID-T2 database has provided us with standardized dimensions for the 3-D docking grids for both class I and class II pMHC structures. These standardized values were used to set the dimensions of the docking grids in all our experiments. Unlike the previously reported grid-based docking method [34], homology model building for MHC receptors has not been used in the development of pDOCK, instead using only experimentally determined X-ray crystallographic structures. The pDOCK method, however, is generic and is applicable to high quality homology models of alleles when experimental structures are not available. Here, the receptor modeling sub-step mentioned in the preparatory step 1 (Figure 1) can be used in the absence of structural data for the MHC proteins. Thus, the direct use of X-ray crystal structures in our docking simulations ensures accurate results.

The first experiment that we conducted was to ensure that an extended peptide bound to its cognate MHC receptor preferentially selecting the same nonameric core peptides as in the crystal structure and then to evaluate the accuracy of the docked peptide. Hence, we performed re-docking of 186 peptides back to their cognate MHC receptors to check for conformational accuracy of the predicted binding registers and their Cα RMSD against their respective crystal structures. We have then benchmarked pDOCK with our earlier docking protocol [2, 3] for a dataset of 50 selected (35 MHC-I and 15 MHC-II) pMHC complexes to verify the speed and accuracy of pDOCK against our earlier method. This was followed by validation and accuracy checks for pDOCK against available flexible peptide docking results obtained from the literature for a dataset of 15 peptides.

In the process of selecting immunogenic peptides for vaccine design, the two main aspects are to determine: (1) multiple peptides that bind to the same allele or MHC molecule and; (2) promiscuous or same peptides that bind multiple alleles. Therefore, as a secondary experiment, we have pursued to test the efficacy and robustness of our docking protocol in modeling the bound conformations of novel peptides to specific MHC alleles by carrying out docking of multiple peptides to a single MHC template structure (same MHC allele), suitable for immune epitope prediction from an antigenic protein, using a moving window of 9-mers along the entire sequence [3, 5, 6, 7]. Our third experiment was to dock a single peptide from particular PDB structures onto multiple MHC templates (multiple alleles) from other crystal structures, suitable for determining promiscuous peptides capable of binding to a set of related alleles and therefore, important for vaccine design.

The Cα RMSD values have been calculated only for the nonameric core of the peptide (for both MHC-I and MHC-II related peptides) which is a contiguous immunogenic segment that forms the “binding register” within the MHC peptide binding cleft, as reported earlier by our group [3]. For the peptides with nine and less number of amino acid residues the entire peptide was used for Cα RMSD calculation. pDOCK accurately detected all 186 binding registers, i.e., the nonameric cores of the peptides are identical to their respective crystal structures. pDOCK generated 85.5% of all the peptides with Cα RMSD of less than 1.00 Å compared to their respective X-ray crystal structures. Our benchmarking results imply up to 2.5 fold improvement in the accuracy of the new peptide docking methodology. The validation results represent a sevenfold improvement in the accuracy of our technique compared to that of the existing methodologies in flexible docking and modeling of peptides into MHC grooves. Amongst the 21 peptides docked in the second experiment, the Cα RMSD values for docked peptides compared to their respective crystal structures were below 1.00 Å for 20 peptides (details in Results and discussion section). The third experiment accounted for all 4 peptides docked with less than 1.00 Å Cα RMSD compared to the same peptides from the corresponding template crystal structures (details in Results and discussion section). Overall, pDOCK is up to 60% faster than our earlier protocol and hence provides a rapid and accurate docking method to evaluate pMHC binding for large scale immune-epitope prediction.

The fact that our earlier method was comparatively slower and that it involved rigid-docking of the peptide termini, acted as the platform for us to ‘revisit’ our pMHC docking methodology. Based on these requirements, we have developed a single step pMHC docking protocol (details in Methods section) as shown in Figure 1, which allows flexibility over the entire length of the peptide antigen and can be used as a generic method to obtain the conformations of bound peptide ligands to MHC binding grooves of both class I and class II MHC proteins. A systematic evaluation of pDOCK is performed as three separate tests: (1) exhaustive re-docking of all non-redundant peptides to their respective MHC grooves as a test case, benchmarking and validation; we then address two very significant practical problems faced by immunologists during the process of allele-specific peptide vaccine design: (2) the docking of multiple peptides that bind to same MHC allele, for immunogenic epitope scanning of antigenic sequences and; (3) docking of promiscuous peptides or same peptides binding to multiple MHC alleles for vaccine design, based on groups of disease-implicated alleles. A correctly docked structure is defined as the peptide with at most 2.50 Å Cα RMSD from the respective experimental X-ray crystal structure [2]. pDOCK has also been benchmarked against our previous docking protocol and validated on published peptide modeling and docking results from the literature. Bordner and Abagyan [34] suggested that while grid-based docking could be applied for pMHC-II, it was a more difficult problem. pDOCK has been successfully applied for MHC-II peptide docking as well with excellent results.

Experiment 1

Re-docking bound peptides to their cognate MHC grooves

pDOCK has been applied on a non-redundant dataset of 186 (149 MHC-I and 37 MHC-II) pMHC complexes from the MPID-T2 database (details in Methods section, data and docking results in Additional File 1 – Table S1). Initially, the peptides were extracted from the experimental pMHC complexes, randomized and set to extended conformations. This was followed by optimization of the peptide ligands and re-docking of the separated peptides back to their respective MHC grooves. As depicted in Figure 2, our technique generated 159 out of 186 peptides with Cα RMSD values less than 1.00 Å: 124 out of 149 peptides (83%) and 35 out of 37 peptides (~95%) for class I and class II MHC proteins, respectively. ~15% (22/149) and ~1% (2/149) of the peptides have their Cα RMSD values within the ranges 1.01-2.00 Å and 2.01-2.50 Å, respectively amongst the MHC-I peptides docked (Figure 2a). Similarly, ~5% of the peptides have their Cα RMSD values within a range of 1.01-1.30 Å amongst the MHC-II related peptides that were docked using pDOCK (Figure 2b). On an average, pDOCK resulted in a Cα RMSD value of about 0.56 Å for re-docking of peptides into their respective MHC grooves over the entire dataset of 186 pMHC complexes.

Our best results are shown in Figure 3, with structural comparison between the lowest energy docked conformation and the native conformation of the bound peptides for MHC-I (PDB code 1s7q) and MHC-II (PDB code 1d5x) structures. These docked conformations of peptide structures have the best Cα RMSD values of 0.09 Å and 0.11 Å respectively, obtained over the entire dataset. The MHC-II peptide in Figure 3b has 5 out of 6 amino acid residues replaced by amino acid analogues (chemical mimics) in the crystal structure. Nonetheless, it has the best Cα RMSD value among all the MHC-II related peptides used in this study, supporting pDOCK’s applicability to peptide or peptide analogues (containing amino acid mimics in structure-based drug design). pDOCK also generated the least energy docked orientations for all the peptides with accurate determination of their respective binding registers, i.e. having the exact nonameric core in the binding grooves, with respect to their native bound conformations in the X-ray crystal structures. All peptides except one from the class I pMHC crystal structure (PDB code 2gtw; Cα RMSD of 3.08 Å) were within the acceptable 2.50 Å Cα RMSD from their respective native conformations (Figure 2a). Also, none of the MHC-II related peptides showed any deviation from the acceptable 2.50 Å Cα RMSD threshold (Figure 2b).

We carefully examined the re-docked conformation of the peptide LAGIGILTV in the MHC groove of the complex 2gtw, with the X-ray structure. In 2gtw, peptide residues 1 to 5 interact with a formic acid molecule, which was not explicitly introduced into the docking simulation. When the formic acid molecule was included in the docking simulation, the predicted orientation of the peptide using pDOCK is energetically more favourable for pMHC complex formation than the predicted conformation when the formic acid molecule is omitted. The improvement in accuracy by the inclusion of the formic acid molecule is ~13 folds. This is portrayed in Figure 4 which clearly indicates that the peptide residues Leu 1, Ala 2, Gly 3, Ile 4 and Gly 5 that are not correctly predicted in the absence of the formic acid molecule (Figure 4a), are accurately docked when the formic acid molecule is introduced into the docking simulation (Figure 4b), resulting in an improvement in the Cα RMSD value from 3.08 Å to 0.24 Å.

Although water molecules and other common biological ions such as phosphate and chloride may mediate pMHC interactions in some cases, they were omitted from our experiments because the significance and contributions of these molecules towards pMHC binding vary immensely between different peptides and specific alleles over a large dataset like the one used in this study (186 complexes). Our previous protocol achieves a Cα RMSD of 1.53 Å for the bound structure of the peptide from pMHC complex 1jf1, due to the presence of a water molecule positioned around the peptide residues 5 to 7 in the crystal structure leading to erroneous prediction of the loop formed, which resulted in incorrect positioning of interacting residues [2]. However, pDOCK successfully overcomes this restriction to accurately predict the least energy bound conformation of this peptide with a Cα RMSD value of 0.30 Å. The enhancement in accuracy of docking is a direct consequence of the improved sampling of available conformational space in pDOCK. This preliminary experiment is a critical first step as it establishes the validity of our approach and helps us test the ability of our technique to accurately dock cognate peptides into their respective MHC receptors, using the proposed single-step docking procedure.

Benchmarking with our previous methodology

In order to ascertain the improvement in speed and accuracy of pDOCK compared to the old technique, we have benchmarked our peptide docking methodology with our earlier pMHC docking protocol [2, 3] over a subset of 50 pMHC complexes (35 MHC-I and 15 MHC-II) from the complete non-redundant dataset (listed in Additional File 1 – Table S1) and the results are presented in Table 1. pDOCK results are consistently better than our earlier docking methodology in terms of accuracy (Cα RMSD) of the modelled or docked peptide compared to their X-ray crystal structures after docking back into their respective MHC grooves. The new protocol also generates the least energy docked conformations for all 50 peptides with Cα RMSD values less than 1.00 Å, compared to eight peptides, docked using the earlier method, having Cα RMSD values above 1.00 Å (graphically shown in Additional File 2 – Figure S1). The new procedure outmatches the old protocol particularly well for complexes 1s9y, 1hhh, 1jf1, 1e27, 1jpf, 1qo3, 1wbz and 1g7p (Table 1) amongst the MHC-I structures and for structures 1uvq and 1aqd (Table 1) amongst the MHC-II structures (highlighted in yellow in Additional File 2 – Figure S1).

S. No.	Allele	PDB	Peptide Length	Peptide Sequence	Cα RMSD (Å)
					Previous method	pDOCK
MHC-I
1	HLA-A*1101	1qvo	10	QVPLRPMTYK	0.53	0.24
2	HLA-A*0201	1qr1	9	IISAVVGIL	0.46	0.29
3	HLA-A*0201	1akj	9	ILKEPVHGV	0.87	0.39
4	HLA-A*0201	1i1y	9	YLKEPVHGV	0.70	0.66
5	HLA-A*0201	1i7r	9	FAPGFFPYL	0.59	0.47
6	HLA-A*0201	1i7u	9	ALWGFVPVL	0.32	0.29
7	HLA-A*0201	1oga	9	GILGFVFTL	0.32	0.16
8	HLA-A*0201	1qsf	9	LLFGYPVAV	0.54	0.34
9	HLA-A*0201	1lp9	9	ALWGFFPVL	0.58	0.26
10	HLA-A*0201	1s9y	9	SLLMWITQS	1.09	0.39
11	HLA-A*0201	1hhh	10	FLPSDFFPSV	1.10	0.49
12	HLA-A*0201	1jf1	10	ELAGIGILTV	1.53	0.30
13	HLA-B*0801	1agc	8	GGKKKYQL	0.28	0.23
14	HLA-B*0801	1mi5	9	FLRGRAYGL	0.42	0.37
15	HLA-B*2705	1ogt	9	RRKWRRWHL	0.51	0.18
16	HLA-B*2705	2a83	9	RRRWHRWRL	0.55	0.18
17	HLA-B*3501	2cik	9	KPIVVLHGY	0.74	0.26
18	HLA-B*3508	3bwa	8	FPTKDVAL	0.56	0.26
19	HLA-B*5101	1e27	9	LPPVVAKEI	1.27	0.18
20	HLA-B*5301	1a1m	9	TPYDINQML	0.59	0.28
21	HLA-Cw*0401	1im9	9	QYDDAVYKL	0.49	0.34
22	HLA-G*0101	2dyp	9	RIIPRHLQL	0.43	0.16
23	H2-Db	1fg2	9	KAVYNFATC	0.25	0.19
24	H2-Db	3buy	9	LSLRNPILV	0.63	0.23
25	H2-Db	1yn7	10	SSLENFAAYV	0.62	0.14
26	H2-Db	1jpf	11	SGVENPGGYCL	1.14	0.36
27	H2-Dd	1qo3	10	RGPGRAFVTI	1.49	0.17
28	H2-Kb	1t0m	8	SSIEFARL	0.66	0.21
29	H2-Kb	1vac	8	SIINFEKL	0.32	0.22
30	H2-Kb	1wbz	9	SSYRRPVGI	0.89	0.19
31	H2-Kb	1s7q	9	KAVYNFATM	0.20	0.09
32	H2-Kb	1g7p	9	SRDHSRTPM	0.97	0.17
33	H2-Kd	1vgk	9	SYVNTNMGL	0.86	0.25
34	H2-Kk	1zt1	8	FEANGNLI	0.57	0.45
35	H2-Ld	2e7l	9	QLSPFPFDL	0.37	0.35
MHC-II
36	HLA-DQB1*0602	1uvq	20	MNLPSTKVSWAAVGGGGSLV	1.09	0.23
37	HLA-DRB1*0301	1a6a	15	PVSKMRMATPLLMQA	0.38	0.30
38	HLA-DRB1*0101	1aqd	14	GSDWRFLRGYHQYA	1.08	0.28
39	HLA-DRB1*0101	1fyt	13	PKYVKQNTLKLAT	0.68	0.23
40	HLA-DRB1*0101	2iam	15	GELIGILNAAKVPAD	0.56	0.24
41	HLA-DRB1*0401	1d5x	6	XXRXXX	0.23	0.11
42	HLA-DRB1*0401	1d5z	7	XXRAXSX	0.33	0.22
43	HLA-DRB1*0401	1d6e	8	XXRXMASX	0.32	0.14
44	HLA-DRB1*0401	1j8h	13	PKYVKQNTLKLAT	0.59	0.20
45	HLA-DRB3*0101	2q6w	11	AWRSDEALPLG	0.54	0.30
46	HLA-DRB5*0101	1fv1	20	NPVVHFFKNIVTPRTPPPSQ	0.88	0.59
47	I-Ad	1iao	14	RGISQAVHAAHAEI	0.81	0.27
48	I-Ak	1iak	13	STDYGILQINSRW	0.42	0.23
49	I-Au	2pxy	11	RGGASQYRPSQ	0.78	0.28
50	I-Ek	1r5v	13	ADLIAYPKAATKF	0.82	0.28

Cα RMSD values are calculated only for the nonamer binding cores (shown in bold) for peptides with more than 9 residues in the X-ray crystal structures. X: Amino acid analogues (chemical mimics). A graphical representation of the results is available in Additional File 2 – Figure S1.

These results suggest that some of the conformational limitations of our previous methodology, such as the presence of water molecules in and around the peptide and within the peptide binding groove in the original PDB structure, have been addressed in our new docking protocol making it highly accurate. Besides an improvement in the accuracy, pDOCK is also able to accurately model docked conformations for some peptides especially for MHC-II related peptides with more than 9 amino acid residues, thereby improving the coverage over the entire length of the peptides. Peptides from the pMHC complexes 1uvq and 2iam were among the highest coverage (20 and 15 residues respectively) obtained in this experiment with Cα RMSD values 0.42 Å and 0.46 Å respectively over the length of the entire peptide (results not shown). The reliability for the accurate prediction of flanking residues (especially for MHC-II peptides) depends upon their interactions with the MHC residues outside the peptide binding groove and therefore, have not been included in the calculation of Cα RMSD values reported.

In terms of the computational time to complete a single docking experiment, pDOCK is up to 60% faster (on an average) than the earlier method as summarized in Table 2. The average time taken by pDOCK is approximately 10 min. (the preparatory receptor positioning step of ~3 sec. {0.50%}, determining the docking grid taking ~42.6 sec. {7.10%} and the single docking and refinement step of ~9.24 min. {92.4%}), compared to 23.50 to 24.50 min (Step 1 taking ~5 min., Step 2 of ~30 sec., Step 3 taking ~18 min. and Step 4, which was only applicable to MHC-II related peptides, of ~1 min.) using the old protocol on a 2 CPU 3.20 GHz 3 GB RAM workstation. The average time taken for each of the steps using either of the methodologies is calculated over the entire non-redundant dataset of 186 pMHC complexes catalogued in additional file 1 – table S1. The mean Cα RMSD value for the least energy docked conformations of peptides, from the dataset of 50 peptides used for benchmarking, was 0.27 Å using pDOCK compared to 0.65 Å for the old procedure. This denotes almost two and a half fold improvement in the accuracy of our novel docking strategy over a larger dataset (50 peptides) than that used previously (40 peptides) [2].

Previous method	pDOCK
Step 1: ~ 5 min
Step 2: ~ 30 s	Preparatory Step 1: ~ 3*s
Step 3: ~ 18 min	Preparatory Step 2: ~ 42.6 s
Step 4#: ~ 1 min	Single docking and refinement step: ~ 9.24 min
Total: ~ 23.50 – 24.50 min	Total: ~ 10 min

Both methodologies were applied using a 2 CPU 3.20 GHz 3 GB RAM workstation. *Only for X-ray crystal structures of MHC proteins. The time taken for this step would increase if homology modeling needs to be carried out. ^#Applicable only to MHC-II related peptides.

Validation against previously published studies

Keeping in mind the essence of improving the accuracy and robustness of the proposed strategy, we have validated pDOCK with seven studies involving MHC-I peptide docking/modeling and one study involving MHC-II peptide docking, covering 15 pMHC structures and compared the results by re-running our earlier method. The results of our validation experiments are compiled into Table 3. Peptides 1, 2, 3, 4 and 15 (Table 3) are new in this study and are collated from recent publications [34, 36, 37], whereas the remaining 10 were from the validation studies reported for our earlier methodology [2]. To the best of our knowledge, these results represent a sevenfold increase in the accuracy of pDOCK compared to available flexible docking techniques in the remodeling of pMHC complexes. Interestingly, the validation criteria for almost all of the previously published studies [34, 36, 37, 38, 39, 40] involved either docking or remodeling of peptides back into their original crystal structure. Although the Cα RMSD values (0.29 Å and 0.30 Å, respectively) for peptides 2 and 3 (Table 3) were slightly higher, they are still comparable with the Cα RMSD values reported earlier (0.23 Å and 0.22 Å, respectively) [34]. Peptide 1 (Table 3) however, was generated with a better Cα RMSD (0.31 Å) compared to the Cα RMSD (0.76 Å) reported in the same earlier grid-based docking study [34]. The enhancement in the accuracy for peptide 1 could be a direct implication of more conformational sampling space in a flexible environment resulting from a relatively larger docking grid (35.36 Å x 35.52 Å x 35.79 Å) for MHC-I peptides and a lower temperature (300 K) used in pDOCK compared to the grid dimensions (34 Å x 34 Å x 25 Å) and temperature (700 K) used in the previous grid-based docking study [34]. Thus, pDOCK is not only comparable to but also surpasses the available techniques in flexible docking and remodeling of peptides with regards to the accuracy (Cα RMSD) with which it predicts the bound structure of a peptide to its respective MHC groove. By and large, our results illustrate the advantages of using grid-based flexible docking over conventional docking protocols.

S.No	Technique	Peptide Sequence	MHC class	PDB	RMSD (Å)
					Published	Previousmethod	pDOCK
1	Grid-based Flexible docking [34]	RGYVYQGL	I	1kpu#	0.76	0.59	0.31
2	Grid-based Flexible docking [34]	ALWGFVPVL	I	1i7u	0.23	0.32	0.29
3	Grid-based Flexible docking [34]	ELAGIGILTV	I	1jf1	0.22	1.53	0.30
4	Monte Carlo annealing [37]	LLFGYPVYV	I	1duz#	3.01	0.33	0.33
5	Simulated annealing [38]	FLPSDFFPSV	I	1hhh	1.59	1.10	0.48
6	Simulated annealing [38]	GILGFVFTL	I	1hhi#	0.46	0.32	0.16
7	Simulated annealing [38]	ILKEPVHGV	I	1hhj#	0.87	0.87	0.55
8	Simulated annealing [38]	LLFGYPVYV	I	1hhk#	0.78	0.33	0.33
9	Combinatorial buildup algorithm [39]	RGYVYQGL	I	2vaa#	0.56	0.32	0.22
10	Combinatorial buildup algorithm [40]	LLFGYPVYV	I	1hhk#	1.40	0.33	0.33
11	Combinatorial buildup algorithm [40]	ILKEPVHGV	I	1hhj#	1.30	0.87	0.55
12	Combinatorial buildup algorithm [40]	GILGFVFTL	I	1hhi#	1.60	0.32	0.16
13	Multiple copy algorithm [41]	FAPGNYPAL	I	2vab#	2.70	0.40	0.25
14	Multiple copy algorithm [42]	GILGFVFTL	I	1hhi#	1.40	0.32	0.16
15	GOLD/GLIDE Flexible docking [36]	XXRXMASX	II	1d6e	1.24/3.06	0.32	0.14

X: Amino acid analogues (chemical mimics). ^#These structures are not listed in Additional File 1- Table S1 due to redundancy in MPID-T2.

Figure 5 provides a pictorial representation of an example of the above discussed accuracy. This structural comparison between the least energy docked conformation generated using pDOCK and that of the native conformation of the cognate peptide in the complex 1duz portrays not only the highly accurate predicted conformation of the peptide, Cα RMSD of 0.33 Å compared to that of 3.01 Å reported earlier [37], but also highlights the fact that the peptide’s N-terminal residues (Leu 1, Leu 2 and Phe 3) were better modeled and structurally well aligned to that of its native conformation when compared to the lowest energy docked conformation reported earlier [37]. Notably, the least energy docked conformations generated for a common murine MHC (H2-Kb) related Sendai virus nucleocapsid peptide FAPGNYPAL and a very familiar human HLA (A*0201) related Influenza A virus matrix peptide GILGFVFTL have significantly lower Cα RMSD values of 0.25 Å and 0.16 Å respectively (Table 3) than those reported in earlier studies (2.70 Å and 0.46 Å, 1.60 Å, 1.40 Å respectively) [38, 40, 41, 42] and those obtained using our previous protocol (0.40 Å and 0.32 Å). These observations establish the efficacy of pDOCK to dock highly accurate multi-species related peptide structures permitting conformational sampling of the peptide in the binding groove during flexible docking.

We have developed pDOCK as a fast, accurate and robust method for flexible docking of peptides to MHC-I and MHC-II proteins. Our results provide evidence of overcoming limitations pertaining to the application of our previous methodology, such as the presence of water molecules in and around the peptide and within the peptide binding groove in the template and relatively longer computational time required. Benchmarking with our previous method for a dataset of 50 non-redundant pMHC complexes consistently produced least energy docked conformations of peptides below 1.00 Å Cα RMSD from their respective native orientations for all 50 peptides. The Cα RMSD range for the same dataset was 0.09 Å (1s7q) to 0.66 Å (1i1y) using pDOCK compared to a Cα RMSD range from 0.20 Å (1s7q) to 1.53 Å (1jf1) applying our previous protocol. These observations imply an improvement in the accuracy by upto two and a half folds compared to our previous protocol. The outcomes of our validation experiments suggest a seven-fold improvement in the accuracy of the pDOCK docking protocol. pDOCK can therefore be successfully applied as a generalized, efficient protocol for docking of peptides to MHC-I and MHC-II receptors with improved accuracy, greater coverage of peptide residues and vastly reduced computational time (up to 60% compared to our earlier method).

The average time taken to perform each step using pDOCK has also improved drastically compared to our old technique on a 2 CPU 3.20 GHz 3 GB RAM workstation. This is mainly due to the consolidation of the docking and refinement protocols into a single step docking and refinement procedure. Our results establish the efficacy of pDOCK to model highly accurate pMHC complex structures permitting conformational sampling of the peptide in MHC binding groove. The current study thus presents one of the most accurate pMHC docking protocols developed to date. pDOCK targets a more generic approach to generation of docked conformations of peptides using a single template for each allele. For some pMHC complexes however, appropriate addition of mediating molecules or considerations of solvent effects may lead to a possible improvement in docking accuracy. Rapid and large scale docking and scanning for identification of potential candidates for immunogenicity from repertoires of immunologically significant antigenic peptide sequences is possible by automating all the steps. No requirement for experimental data to be trained and the need of only a suitable template for a particular allele give pDOCK a prominent edge over other sequence-based techniques such as Artificial Neural Networks, Support Vector Machines, and Hidden Markov Models.

pDOCK is also highly efficient in accurately predicting the docked conformations of amino acid analogues or chemical components within the peptide ligand suggesting its possible use as a docking and evaluation tool in structure-based drug design protocols and chemoinformatics. The single and variant template docking experiments along with the validation experiments also serve as strong benchmarks for pDOCK against our old method. pDOCK can correctly predict the conformation of residues that extend into the MHC binding cleft and therefore could help identify essential contacts with the MHC receptor, responsible for reducing the half life of the pMHC complex such that the peptide is held long enough within the MHC groove for presentation at the APC cell surface leading to surveillance and recognition by the TR molecules which in turn results in the activation of T cells and triggers the adaptive immune response cascade. Another significant improvement in this study is that the peptide ligand is allowed full flexibility within the peptide biding groove of the MHC proteins, unlike our previous method where the peptide termini were docked rigidly to the MHC groove. This aspect of pDOCK has helped us carry out fully flexible peptide docking to the MHC proteins. Our results also indicate the successful application of this protocol for easy in silico identification of promiscuous peptide epitopes that are applicable to higher proportions of human population with greater propensity to bind to MHC proteins and consequently activate T cells making them key targets for the design of vaccines and immunotherapies.