Docking and Structural Properties of the Peptides
In this study, residues of 1-283 refer to the HLA-A*03 and residues of 55-85 refer to the α1-helix and residues of 139-175 represent the α2-helix of the HLA-A*03 allele, and residues of 275-283 show the PLP immunodominant epitope in HLA-A*03- PLP[45–53] complex (p-HLA complex) (Fig. 1). Also, residues of P1-P19 refer to positions 1-19 of each peptide. P-HLA complex in the presence of each peptide (triple complexes) is named as complexes of 1 - 14.
Docking of the peptide into the p-HLA binding site is one of commonly used molecular modeling methods for determining the desired orientation of each peptide. Totally, 14 designed peptides were docked into the p-HLA complex at the TCR-binding interface. The obtained active residues for HADDOCK server explained in Methodology Section were residues of 58, 62, 63, 64, 65, 66, 67, 69, 71, 72, 73, 74, 75, and 76 from α1 helix and residues of 146, 147, 152, 153, 154, 155, 157, 158, 159, 161, 163, and 167 from α2 helix and residues of 275-283 from PLP neuroantigen. Also, all 19 residues of the designed peptides were considered as active residues for the HADDOCK server. The predicted binding score is presented in Table 1 for each peptide calculated by the HADDOCK web server.
The results of docking showed that the peptide 14 (WRYWWKDWAKQFRQFYRWF) had the highest binding affinity than other peptides. This relatively high binding affinity exhibited by peptide 14 can be significantly attributed to replacement of Glu residues with Tyr at location 3 and Phe with Ala at position 9 of peptide 1 (WREWWKDWFKQFRQFYRWF). Of course, residue 3 is more important than residue 9, because it has more contribution in hydrogen binding formation between peptide and p-HLA complex. Selected peptides were assessed for grand average hydropathicity (GRAVY) and antigenic region. Positive and negative GRAVY values represent hydrophilicity and hydrophobicity, respectively. All peptides had a strong degree of hydrophilicity (negative GRAVY value) without any antigenic region in the peptides҆ sequence (Table 1). Increased peptide hydrophilicity may lead to passage across the blood-brain barrier (BBB) membrane.
Results of MD Simulations
Sometimes, docking results are unreliable because some of the best-docked conformations become separate from the protein within a short time of MD simulation. Thus, herein, MD simulations were performed to make sure that the docked peptide stays at the p-HLA binding site.
Changes in root mean square deviation (RMSD) relative to their initial structure during simulation showed stability of backbone atoms (Fig. 2(a)).
The average values of backbone RMSD of p-HLA are summarized in Table 2. P-HLA complex in the presence of each peptide was characterized with an initially high RMSD fluctuation followed by a decrease in values of backbone RMSD during the last 30 ns (Fig. 2(a)). The highest backbone RMSD of p-HLA belonged to system 14 (p-HLA in complex with peptide 14). However, the average backbone RMSD value of the p-HLA complex indicated strong fluctuations in certain regions. Then, RMSD of the residues forming the α1-helix (Fig. S2(a)), α2-helix (Fig. S2(b)), and PLP-bound neuroantigen (Fig. S2(c)) of the p-HLA complex was calculated separately (Table 2). Comparing Fig. 2(a) with Figs. S2(a-c) showed that RMSD values of these three regions were much lower than those of the p-HLA in the complex. Since, the alpha chain of class I HLA molecules consists of three domains of α1, α2, and α3, it was concluded that high fluctuation observed in Fig. 2(a) will be attributable to α3 helix. In this study, the β2M chain was omitted to decrease the computational cost. The absence of the β2M domain in simulation led to high RMSD of the α3 domain. In class I of HLA molecules, it was found that β2M interacts with all three domains of α1, α2, and α3 whereas the α3 domain is mostly responsible for binding to β2M [55]. As shown in Figs. S2(a-c), the three mentioned regions were distinguished by RMSD fluctuations of less than 2.5 Å, followed by a relative decrease in RMSD. After approximately 30 ns, RMSD values converged to about 2 Å.
Furthermore, RMSD of all peptides experienced various degrees of fluctuations at first but gradually tended to converge (Fig. 2(b)). Also, probability distribution of RMSD values over the last 10 ns of simulation was calculated for all RMSD values (data not shown) and the Gaussian shape of probability showed no statistical errors in simulation time and sampling.
Besides, in 14 simulated complexes, the RMSD values of the backbone atoms of PLP peptide varied between 0.2 - 1.3 Å (Fig. S2(c) and Table 2). Since, the PLP was surrounded by HLA residues from all directions, α1 and α2 helices and beta floor of HLA, and peptides҆ residues, there was hardly a space for conformational movements and fluctuation in MD simulation (Table 2).
Root Mean Square Fluctuation (RMSF)
For understanding fluctuations in residues and investigating the effect of peptide binding on p-HLA complex dynamics, RMSF of residues during the last 10 ns of MD trajectories was calculated and the results are presented in Fig. 3(a).
Also, using the formula mentioned in the literature [56], overall fluctuations (Foverall) or total flexibility for the three described regions, peptides, and p-HLA complexes were calculated. Binding peptides led to different Foveral values in the p-HLA complex and change in the overall structure of p-HLA (Fig. S3, Table 2). In the presence of peptide 14, the p-HLA complex had the highest Foveral value relative to other peptides. In other words, peptide 14 induced the most structural fluctuation in p-HLA complex. P-HLA complex in interaction with peptide 11 had the lowest Foverall value and there were slight changes in the p-HLA complex vibration.
Also, interaction of different peptides with PLP peptides in the HLA binding site led to different fluctuations in PLP peptides. In fact, each peptide represents particular new positioning of PLP (antigenic peptide presentation mode) in HLA-A*03 binding grooves. Fig. S3 provides a comparison of the PLP conformation in the HLA-A*03 binding groove in the presence of each peptide. The least and the most Foveral values of PLP belonged to binding peptides of 7 and 12, respectively. In other words, binding peptide of 12 led to more fluctuation or more structural change in PLP in the binding groove of HLA relative to other peptides. Madura in a study showed TCR-induced alteration in interaction of HLA with anchoring antigenic peptide (PLP). T cell antigen recognition required an anchor residue shift at N-terminus of the antigenic peptide [57]. Results of another study showed that latter portion of the antigenic peptide is “lifted” from the binding groove by 1–1.5 Å upon TCR binding, beginning at residue of 4 and continuing through the C-terminal residues [11]. Our results also indicated that binding of the designed peptides changes Foveral of PLP peptide (Residues of 275-283 in Fig. 3(a)). It may reflect weaker interactions of peptides with class I HLA proteins at their C-terminal ends [58]. Similarly, in the presence of peptides 3 and 11, the C-terminal of PLP detaches from their binding pockets and stretches outside the HLA-A3 binding groove because, these peptides make strong hydrogen bonds with C-terminal of PLP (Table S3).
Motions or dynamics of α1-α2 helices and β-strand floor are primarily influenced by each peptide binding [58]. Our results showed that the RMSF values of α1 and α2 helices and binding groove (α1+α2+β-strand floor) changed in binding of different peptides (Table 2). In the presence of peptides 11 and 13, HLA-A*03 binding groove showed the lowest and most Foverall values, respectively. Then, peptide 13 induced the most structural changes in binding groove of HLA relative to other peptides. Knapp et al., showed that only flexibility of α2-1 and α2-2 residues changed in interaction with TCR [59].
Also, RMSF of the designed peptides was calculated during the last 10 ns of MD simulation and the results are presented in Fig. 3(b) and Table 2. The least Foveral value of peptides belonged to peptides 7 and 11. Then, these peptides had the least fluctuation in binding to p-HLA complex. Also, the most fluctuation belonged to peptide 12.
Radius of Gyration (RG)
In the presence of each peptide, level of compaction of α1 and α2 helix of the HLA-A*03 molecule and PLP was determined by RG (Rg) so that, lower Rg indicates more compactness of the p-HLA complex. Fig. 4 and Table 3 show RG of p-HLA complex and each peptide measured over the last 10 ns of simulation. The results indicated partial structural changes in p-HLA complex in binding of different peptides. According to Fig. 4, it was concluded that the p-HLA complex was more compact while binding to peptide 11 relative to other peptides. Compaction of p-HLA complex in binding to peptide 11 was compatible with the lowest Foveral value of p-HLA complex for complex 11 (Table 2). In other words, compaction can lead to low fluctuation.
According to Table 3, the PLP was partially uncompact while binding of peptide 14 to p-HLA complex relative to other peptides. The most Rg value of PLP peptide was seen in triple complex 14 that was compatible with the highest Foveral value in PLP peptide in complex with peptides 12 and 14 (Fig. S4 and Tables 2 and 3).
Average values of Rg were less for α1-helix than α2-helix as the second α2-helix was not completely straight and contained kinks at residues of 150 and 165 such that, one can distinguish three α-helical sub-segments (Fig. 1). The most Rg and Foveral values of α2 helix belonged to complex containing peptide 13. Then, these parameters were compatible in α2 helix. (Fig. S4 and Tables 2 and 3).
Fig. S4 and Tables 2 and 3 show that binding of the designed peptides can lead to induction of different tertiary structures in p-HLA complex and PLP and α1 and α2 structures. Radius of the designed peptides was similar, but the least and most average Rg values of the designed peptides belonged to peptides 8 and 7 (Table 3).
Solvent Accessible Surface Area (SASA)
SASA shows the amount of exposure of a given region to solvent medium. Figs. S5(a-c) depict SASA values of α1 and α2 helices, and PLP versus simulation time, also the average values of SASA are summarized in Table 4. Here, decreased values of SASA denote to coverage of the mentioned regions by peptides. The results of analyzing SASA suggested that peptides 5, 11, and 1 decreased SASA of PLP more than other peptides, in other words, these peptides covered surface of PLP peptide better than other peptides and could inhibit TCR binding to p-HLA complex. In the presence of peptide 3, PLP peptide had the most SASA and the least coverage with PLP peptide.
The most coverage (the least SASA) and least coverage of α1- helix belonged to peptides 9 and 13, respectively.
Peptides 8 and 14 had the least and most SASA values during the last 10 ns of MD simulation, respectively. Similar to Rg results, the estimated SASA values for the α1-helix were less than α2-helix values. Optimal α2 helix coverage was not feasible due to kinks at residues of 150 and 165.
Center of Mass Distance
The distance between center of mass (COM) of the p-HLA complex and peptides was calculated. Fig. 5 demonstrates the distance based on all Cα atoms plotted versus simulation time. The average COM values of the last 10 ns of the simulation are summarized in Table 5. The least and most COM distance values belonged to peptides 9 and 2, respectively. Of course, the average COM is not a precise parameter, because it depends on peptide size and initial position (obtained from docking). Also, probability distribution of center of distances was obtained (the plot was not shown). The Gaussian shape of this plot showed no statistical error in MD simulation and sampling.
Binding Free Energy
The mentioned parameters, such as RMSD, RMSF, SASA, and COM are qualitative parameters for investigation of peptide binding to p-HLA complex. For investigating binding mode of peptides to p-HLA complex quantitatively, binding free energy was calculated during the last 10 ns of MD simulation after MD simulation. In a previous study, Hou et al., [60] showed that MM/GBSA calculation was more proper for indication of relative binding free energies. Considering computational efficiency of the MM/GBSA approach, it was decided to calculate binding free energy by this approach (Table 5). The best negative binding free energy belonged to peptide 14 with a value of about -75 kcal/mol indicating favorable protein-peptide interaction. The calculated values were overestimated because of the lack of entropic and energetic contributions. Entropic contributions were unfavorable in the case of flexible peptides. Energetic contributions or conformational changes occur due to using the single trajectory approach [61]. Binding free energy results were compatible with docking results.
Also, for investigating contributions of binding free energy of each residue, energies are decomposed into individual residue contributions using the MM/GBSA approach. This type of calculations not only rationalize molecular recognition processes but also can guide identification of residues that mimic determinants of binding in protein-protein complexes (hot-spots)[62]. Contribution of residues in the total binding free energy of each peptide is presented in Figs. S6-18. Only those residues whose contributions to effective energy (ΔGEff) were equal or less than −2 kcal/mol were considered. For understanding the main binding driving forces of peptides to the p-HLA complex, ∆G values were divided into two terms: electrostatic interaction plus polar solvation energies (∆Eele+∆GGB) and Van der Waals interaction energy plus non-polar solvation energies (∆Evdw+∆GSA)[63].
For instance, in system 14, the main energy contribution belonged to the residues of Arg 2, Tyr 3, Arg 17, Trp 8, Trp 4, Trp 18, Phe 12, Phe 15, and Trp 5 of peptide and residues of Arg 65, Phe 281, Glu 58, Gln 62, and Glu 152 of the complex (Fig. 6). Electrostatic interaction and polar solvation energy (∆GGB) played the main roles in binding of the residues of Arg 2 and Arg 17 of peptide, and Glu 152 of the complex. Binding energy contribution of the other residues was mainly due to van der Waals interaction (∆Evdw) and non-polar solvation energy (∆GSA) (Fig. 6).
Analysis of Hydrogen Bond
Hydrogen bonds play an important role in stabilizing peptides inside the HLA binding groove [59]. Loss of residue-residue hydrogen bond connections is a potential explanation of the increased residue flexibility (Foveral) as noted in Table 2. McMahon et al., showed a hydrogen bond network between PLP peptide and HLA-A3 peptide-binding groove in X-ray structure of the p-HLA complex [15]. Hydrogen bonding occupancy in the presence of each peptide between the PLP and HLA residues is shown in Table S2 of supplementary material, throughout 40 ns of MD simulation. Also, the average of number of hydrogen bonds was calculated during the last 30 ns of MD simulation in this study. During MD trajectory, these hydrogen bonds fully changed in the presence of peptides in most of them. For instance, in the presence of peptide 14, it was found that Glu 4 of PLP peptide made a hydrogen bond with Tyr 171 of HLA chain. The newly developed H-bond was defined with high occupancy of 85% and an average H-bond distance of 2.75 Å. McMahon et al., demonstrated a conventional hydrogen bond with Tyr 171 formed with N-terminal of PLP residues [13]. In our study, a pocket was occupied by the Glu 4 of PLP. As can be seen in Table S2, instead of Tyr 7, stated in the study by McMahon et al., N-terminal residue (Lys 275) of PLP preferred to form a hydrogen bond with Arg 6 and N atoms of the binding groove (HB occupancy= 53%).
For studying origin of stronger binding affinity of peptides to p-HLA complex in detail, hydrogen bond was analyzed for all peptides. The average number of hydrogen bonds was calculated between each peptide and the p-HLA complex during simulation (Fig. 7). Peptide 14 exhibited the highest total H-bond occupancy between the peptides and the p-HLA complex (547%) compared to other peptides (Table 5). Calculation of the average of number of hydrogen bonds between peptide 14 and p-HLA complex during simulation time indicated that tyrosine substitution at location 3 of peptide 14 provided appropriate distance and angle for growth of more hydrogen bonds over simulation time (Figs. 7 and 8). Within peptide 14, it was found that the formed hydrogen bonds between Arg 2 of peptide and Arg 65 from HLA’s α1- helix exhibited high H-bond total occupancy of 130 % with an average H-bond distance of 2.8 Å. As can be seen in Table 6, the stated hydrogen bonds were created by the backbone carbonyl group of Arg 2 whereas, Arg 2 and sidechain atom formed hydrogen bonds with Glu 58 atoms from HLA’s α1-helix (H-bond occupancy= 69 %). H-bond distance was responsible for high bond strength in peptide 14 complex. There were five hydrogen bonds in peptide 14 that reached more than 50% of their occupation while 3, 3, 2, 2, and 3 hydrogen bonds with the mentioned occupancy percentages were found in peptides of 11, 10, 12, 5, and 3, respectively.
H-bond distance and occupancy of 14 respective peptides throughout simulation time are presented in Table S3 for other peptides. In peptide 3, two of three strong hydrogen bonds resulted from interaction of C-terminal residues of PLP with Arg 2 of peptide 3.
Similarly, in peptide 11, hydrogen bonds between C-terminal of PLP (Lys 283) and Arg 2 of peptide 11 played an important role in hydrogen bond strength and increase in total occupancy (Table S3).
Strength of the mentioned hydrogen bonds was well confirmed by free energy decomposition analysis (Fig. S5-S17). In peptides of 3 and 11, the total binding energy contribution of Arg 2 was equal to -14.3 and -16.2 kcal/mol, respectively.
As shown in Table S2, in the presence of peptide 7, not only C-terminal of PLP was properly positioned inside the F-pocket (ASP 74 and ASP 77) but also position 3 of PLP was able to maintain the hydrogen bond described in the study by McMahon et al., and hydrogen bond was formed between Ile3 of PLP and Tyr 99 of C-pocket confirming low Foveral value of PLP in binding of peptide 7 to p-HLA complex (Table 2).
Without taking into account entropic contribution, the calculated MM/GBSA energies were equal to the available experimental enthalpies. MM/GBSA binding free energy of peptide 14 with p-HLA was about -75.1 kcal/mol, which is in fair agreement with the TCR-pHLA complex enthalpies, usually ranging from −10 to −30 kcal/mol, measured by van't Hoff analysis or calorimetric experiments [64] and theoretical procedures. Also, enthalpies values reported for TCR-HLA-A2-(MP58–66) and TCR-HLA-A2-SL9 ternary complexes reactions were equal to -23 and -10.4 kcal/mol, respectively [65]. Non-covalent interactions at PPI sites consist of van der Waals interactions, hydrogen bonds, and salt bridges. Interaction surfaces of both TCR and HLA–peptide complex are large so that, approximately a thousand square angstroms are buried when an interaction event occurs and most of the contacts are mediated by van der Waals interactions. Both hydrogen bonds and salt bridges, which are comparatively rare in TCR and p-HLA PPI interfaces, add specificity because of their dependence on stereochemistry and geometry of the interactions. Although, TCR-p-HLA interactions bury large amounts of surface area, they seem to have relatively less optimal shape complementarity (dissociation constants in the micromolar range) relative to cytokine-receptor interactions (dissociation constants in the nanomolar range); this low shape complementarity can lead to low average affinity of TCR-p-HLA interactions [66]. On the other hand, PPIs are extremely difficult to target because of their high surface area, and they are normally flat and shapeless, characterizing the interaction [57].
Murray described a conserved HLA motif (a.a., R65-X-X-K-A-X-S-Q72) and found that the TCR’s CDR3α loop hooks up to R65-joint and this position predisposed the HLA-A lineage to TCR alloreactivity. Furthermore, CDR3α loop hooks up to Arg 65-joint and maintains its CDR2α loop at a distance of about 4 Å from polymorphic amino acid positions of the α-2 helix [67]. As a result, blocking Arg 65 of the conserved motif of p-HLA with Arg 2 of peptide 14 (Table 6) led to inhibition of binding of TCR to the p-HLA complex. Besides, an overall H-bond occupancy of 194.7% was observed between Arg 17 of peptide 14 and Glu 152 of α2-helix (Table 6). HLA class I supertypes have been classified according to similar peptide-binding motifs. HLA-A*0301, -A*1101, -A*3101, -A*3301, and -A*6801 alleles belong to the A3-like superfamily [68]. Here, the HLA-A*03 allele sequence was compared with other A3 superfamily alleles and it was found that at position 152, only the HLA-A*03 allele has glutamine residue. Specificity of the peptide14 is explained based on association with sites that are distinctive to HLA-A*03. Furthermore, T cell maturation optimizes the α/β TCR repertoire for recognition of major histocompatibility complex (HLA) class-I polymorphism [67]. The interaction with α2-helix was accomplished via a remarkably compatible conjunction involving Gln155 as stated in the study conducted by Murray. Carbonyl group of Gln155 formed partially hydrogen bonds with Gln 11 and hydrogen atoms of peptide 14) H-bond occupancy: 20%(. Since, proximity of CDR2α loop over HLA (residues of 151–158) is essential for interaction with α2 helix, it appears that due to a strong hydrogen bond with Glu 152, the α2 helix is not recognizable by TCR. Furthermore, Murray stated that HLA-A24, -23 does not have Arg 65, but interestingly, has conjunction involving CDR3α contact with Glu 62 of the α1 helix. As shown in Table 6, a strong hydrogen bond formed between Trp 4 of peptide 14 and Gln 62 of p-HLA complex (H-bond occupancy: 61%). Also, Gln 62 was partially occupied with Trp 5 atoms. Therefore, peptide 14 can form strong hydrogen bonds with polymorphic and non-polymorphic contact sites of TCR. Interaction of HLA’s residue with peptide 14 is compatible with results of a former study on a TCR-pHLA triple complex model [11, 69, 70]. Tripathi identified hydrogen bond patterns and electrostatics in IG4TCR-HLA-A2-tumor epitope (NY-ESO) ternary complexes҆ interface with six hydrogen bonds. In the study by Tripathi, three hydrogen bonds with H-bond occupancy of 90%, one case with H-bond occupancy of 77%, and two cases with H-bond occupancy of 40-50% were observed [71]. The hydrogen bond developed between peptide 14 and the p-HLA complex in the present research showed higher total HB occupancies compared to those reported in the study by Tripathi (total HB occupancies: 547%). This is consistent with free energy decompositions, which revealed that, in general, peptide residues of P2, P17, P8, P4, P18, P12, P15, and P5 had the largest contribution to effective free energy, respectively (Fig. 6(.Studying free energy decomposition also revealed that van der Waals interactions and non-polar part of solvation free energy represent binding strength of the peptides, whereas binding specificity is determined by electrostatic interactions and polar part of solvation free energy, as confirmed in the literature [71]. Also, peptide 14 typically had contact with the exposed side chain of the PLP neuroantigen (Phe 281( and partly with C-terminal residues (Lys 283) (Fig. 6).