Peptide-MHC Surface Morphology Predictions Protocol
[Park et. al. Molecular Immunology 2013] Read Here
1. Peptide Sequence + MHC Sequence
The first step to predicting peptide-MHC structures is inputting a peptide sequence and a specific MHC class I molecule.
The first step to predicting peptide-MHC structures is inputting a peptide sequence and a specific MHC class I molecule.
2. Homology Modeling
Next, we build a homology-based model using MODELLER (1) software. MODELLER (1) requires a template with a previously available X-ray crystal structure. Of all available crystal structures, one with a sequence homologous to our input pMHC complex is preferentially selected as a template.
3. 100 Runs of Simulated Annealing (1ns)
Starting from the predicted structure from MODELLER (1), we next employ all-atom molecular dynamics simulations with the AMBER force field to generate an ensemble of pMHC conformations using a simulated annealing scheme (2): The pMHC, initially configured by the MODELLER (1), is heated from 300K to 1500K during 80ps. Next, it is equilibrated at 1500K for 80ps. Next we slowly cool the system to 300K for 800ps. Finally, we cool again to 283K for 40ps. This scheme is expected to guide the complex to an energy minimized structure. Here we perform Langevin dynamics using the sander module of the AMBER9 with a time step of 2fs. We implement the SHAKE algorithm for proton covalent bonds and introduce a cutoff distance of 10 Å for Lennard-Jones interactions. We treat solvent effects implicitly by adapting the pairwise Generalized Born model as parameterized by Tsui and Case (3). We impose two kinds of restraints: First, a 20kcal/mol harmonic energy barrier is assigned to each of the receptor atoms of the α1 and α2 domains to limit the motion of the receptor. Second, we maintain the distances of the four hydrogen bonds between the peptide and the receptor during the simulation. These four conserved hydrogen bonds, identified by LIGPLOT analysis (4), are: (1) OG1 atom of THR 143 and C-terminal oxygen of the peptide, (2) NE1 atom of TRP 147 and O atom of the peptide at position 8, (3) OH atom of TYR 159 and O atom of the peptide at position 1, and (4) OH atom of TYR 171 and N-terminal nitrogen of the peptide. We repeat this simulated annealing cycle for each pMHC complex until we obtain an ensemble of 100 pMHC conformation candidates.
4. Bayesian Clustering
We use Bayesian clustering (5) to correctly select the most abundant configuration from the ensemble of 100 simulated annealed structures (2). Using pairwise all-atom root mean square deviation (RMSD) values among the 100 structures, we choose the conformation closest to the centroid of the major cluster.
We use Bayesian clustering (5) to correctly select the most abundant configuration from the ensemble of 100 simulated annealed structures (2). Using pairwise all-atom root mean square deviation (RMSD) values among the 100 structures, we choose the conformation closest to the centroid of the major cluster.
5. MD Simulations on 100 Structures (10ns) (100 conformations)
We perform a 10 nano-second all-atom MD simulation (6) using the AMBER force field on each of the 100 structures generated by the simulated annealing: We perform Langevin dynamics simulations at 283K during 10ns using the pmemd module of the AMBER9. The parm99 force field parameters are used with the pairwise Generalized Born implicit solvent model as parameterized by Tsui and Case (3). Prior to the dynamics run, we relax the receptor by steepest descent energy minimization. We employ the same restraints used during the simulated annealing step. All bonds involving hydrogen atoms are constrained using the SHAKE algorithm. A nonbonded cutoff distance of 12 Å and a time step of 2fs is used.
We perform a 10 nano-second all-atom MD simulation (6) using the AMBER force field on each of the 100 structures generated by the simulated annealing: We perform Langevin dynamics simulations at 283K during 10ns using the pmemd module of the AMBER9. The parm99 force field parameters are used with the pairwise Generalized Born implicit solvent model as parameterized by Tsui and Case (3). Prior to the dynamics run, we relax the receptor by steepest descent energy minimization. We employ the same restraints used during the simulated annealing step. All bonds involving hydrogen atoms are constrained using the SHAKE algorithm. A nonbonded cutoff distance of 12 Å and a time step of 2fs is used.
6. Monitoring Conformational Transitions
During the MD simulations6, we monitor the simulation trajectory to determine whether conformational transitions occur by comparing the peptide structure with the starting conformation. We flag all conformations with a shift in all-atom RMSD greater than 0.5 Å as possible conformational transitions. This threshold enables us to focus only on major changes rather than minor motions of peptide side chains.
During the MD simulations6, we monitor the simulation trajectory to determine whether conformational transitions occur by comparing the peptide structure with the starting conformation. We flag all conformations with a shift in all-atom RMSD greater than 0.5 Å as possible conformational transitions. This threshold enables us to focus only on major changes rather than minor motions of peptide side chains.
7. Yes: If there are conformational transitions, we collect the post- transition structures and designate the most prevalent structure as our final prediction of the target pMHC’s surface morphology.
8. No: In the absence of any conformational transitions, we designate the most abundant configuration from the ensemble of 100 simulated annealed structures as our prediction. This structure is defined as the closest structure to the centroid of the major cluster of the ensemble designated by the Bayesian clustering method5.
References
1. Marti-Renom, M.A., Stuart, A.C., Fiser, A., Sanchez, R., Melo, F., Sali, A., 2000. Comparative protein structure modeling of genes and genomes. Annual Review of Biophysics and Biomolecular Structure 29, 291–325.
2. Frenkel, D., Smit, B., 2002. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. Academic Press, San Diego.
3. Tsui, V., Case, D.A., 2000. Molecular dynamics simulations of nucleic acids with a generalized Born solvation model. Journal of the American Chemical Society 122, 2489–2498.
4. Wallace, A.C., Laskowski, R.A., Thornton, J.M., 1995. LIGPLOT: a program to generate schematic diagrams of protein–ligand interactions. Protein Engineering 8, 127–134.
5. Shao, J.Y., Tanner, S.W., Thompson, N., Cheatham, T.E., 2007. Clustering molecular dynamics trajectories: characterizing the performance of different clustering algorithms. Journal of Chemical Theory and Computation 3, 2312–2334.
6. Fagerberg, T., Cerottini, J.C., Michielin, O., 2006. Structural prediction of peptides bound to MHC class I. Journal of Molecular Biology 356, 521–546.
2. Frenkel, D., Smit, B., 2002. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. Academic Press, San Diego.
3. Tsui, V., Case, D.A., 2000. Molecular dynamics simulations of nucleic acids with a generalized Born solvation model. Journal of the American Chemical Society 122, 2489–2498.
4. Wallace, A.C., Laskowski, R.A., Thornton, J.M., 1995. LIGPLOT: a program to generate schematic diagrams of protein–ligand interactions. Protein Engineering 8, 127–134.
5. Shao, J.Y., Tanner, S.W., Thompson, N., Cheatham, T.E., 2007. Clustering molecular dynamics trajectories: characterizing the performance of different clustering algorithms. Journal of Chemical Theory and Computation 3, 2312–2334.
6. Fagerberg, T., Cerottini, J.C., Michielin, O., 2006. Structural prediction of peptides bound to MHC class I. Journal of Molecular Biology 356, 521–546.