In-Silico Design of Peptides for Inhibition of HLA-A*03-KLIETYFSK Complex as a New Drug Design for Treatment of Multiples Sclerosis Disease

Multiple sclerosis is recognized as a chronic inammatory disease. Human leukocyte antigen (HLA) plays an important role in initiating adaptive immune responses. HLA class I is present in almost all nucleated cells and presents the cleaved endogenous peptide antigens to cytotoxic T cells. HLA-A*03 is one of the HLA class I alleles, which is reported as substantially related HLA to MS disease. In 2011, structure of the HLA-A*03 in complex was identied with an immunodominant proteolipid protein (PLP) epitope (KLIETYFSK). This complex has been reported as an important autoantigen-presenting complex in MS pathogenesis. In this study, new peptides were designed to bind to this complex that may prevent specic pathogenic cytotoxic T cell binding to this autoantigen-presenting complex and CNS demyelination. Herein, 14 new helical peptides containing 19 amino acids were designed and their structures were predicted using the PEP-FOLD server. Binding of each designed peptide to the mentioned complex was then performed. A mutation approach was used by the BeAtMuSiC server to improve binding anity of the designed peptide. In each position, amino acid substitutions leading to an increase in binding anity of the peptide to the mentioned complex were determined. Finally, the resulting complexes were simulated for 40 ns using AMBER18 software. The results revealed that out of 14 designed peptides, “WRYWWKDWAKQFRQFYRWF” peptide exhibited the highest anity for binding to the mentioned complex. This peptide can be considered as a potential drug to control multiple sclerosis disease in patients carrying the HLA-A*03 allele. Web-based servers used their characteristic calculations and automatically predicted the hot-spots in inter-subunit interactions of the p-HLA-TCR selected complexes. These computational approaches allowed us to predict possible points of interaction with TCR located on the α 1 and α 2 -helix of α chain of the HLA-A*03 molecule. Interface residues indicated in this step with 9 residues of PLP epitope were selected as binding site for all peptides(cid:0) docking (active residue) in ClusPro and HADDOCK servers. properties


Introduction
Approximately 2.5 million individuals suffer from multiple sclerosis (MS) worldwide, as a persistent in ammatory neurological disorder of the brain and spinal cord, which is a prevalent cause of signi cant physical impairment, fatigue, and pain in young adults, especially women [1].
Differences in clinical symptoms are correlated with spatiotemporal distribution within pathological sites of lesions in the central nervous system (CNS) [2]. These lesions are characteristic of MS and are caused by the blood-brain-barrier (BBB) crossing of immunogenic cells, which causes in ammation, demyelination, gliosis, and neuroaxonal degeneration, nally leading to neuronal death [3]. Macrophages dominate the in ltrate, followed by CD8 + T cells, while fewer CD4 + T, B, and plasma cells can also be found [4]. Many overlapping genetic and environmental risk factors contribute to an immune response to CNS myelin antigens [5]. Several alleles of the human leukocyte antigen (HLA) or major histocompatibility complex (MHC) on the short arm of chromosome 6 [6p21.3] were identi ed as basis of genetic susceptibility to MS [6]. MHC class I (HLA A, B, and C) molecules present endogenous peptides to CD8 + cytotoxic cells. The CD4 + T cells recognize the processed exogenous antigens located on MHC class II (HLA DRB1, DPB1, and DQB1) alleles existing on antigen-presenting cells, such as B cells to stimulate their antibody production [7]. The in uence of both HLA classes of I and II genes on MS susceptibility may re ect important role regarding antigen recognition of HLA class I-dependent CD8 + T cells and HLA class II-dependent CD4 + T cells in pathophysiology of MS [8].
CD8 + T cells have a major role in pathogenesis of autoimmune diseases like MS. CD8 + T cells recognize antigen peptides on MHC class I via the clonotypically-expressed αβ TCR and the lineage-speci c CD8 co-receptor [9]. These cells can then kill target cells presenting the peptide antigens. TCR-HLA crystal structures of the peptide complex, combined with functional experiments have demonstrated that the complementarity-determining region (CDR) loops of TCR engage both the presented peptide and the HLA presentation platform in a generally conserved diagonal orientation, with the TCR-α chain position over the MHC-α 2 helix and the TCR-β chain over the MHC-α 1 helix (Fig. S1) [10,11].
HLA-A3 (HLA-A*03) is one of the class I HLA alleles, which is rmly related to MS. The presence of this allele has been reported as a predisposing risk factor for patients with MS in Khuzestan Province, Iran [12]. Supplementary effect of class I HLA loci (e.g., HLA-A*03) has been described to double the risk of MS [13,14].
In 2011, structure of the HLA-A*03 in complex with an immunodominant proteolipid protein (PLP) epitope (KLIETYFSK) was determined [15]. PLP makes up 50% of whole myelin sheath protein [16]. Protein-protein interactions (PPIs) play pivotal roles in majority of biological processes. Therefore, the improved approaches to target and disrupt PPIs would provide tools for chemical biology and lead to therapeutic development [17]. Peptide-based inhibitors have more advantages than small-molecule drugs. Firstly, targeting protein interactions using small molecules is considered to be di cult due to size of proteins and the lack of well-de ned binding pockets [18,19]. Secondly, the peptide has high speci city towards a target and may not accumulate in organs. This would minimize the side effects and toxicity [20]. Peptides modulating PPIs can be directly derived from crystallographic interface, or through screening of peptide sequences that do not originate from natural proteins using phage display library screening protocol [21].
Besides, computational methods combining simulations of molecular dynamics and calculations of binding free energy could present both structural and energetic outlook of peptide inhibitors to evaluate binding a nity of peptides [21]. In the present study, it is attempted to prevent formation of trimolecular complexes (HLA-A*03-PLP( 45-53 )-TCR) and then, inhibit proliferation of the activated T cells. In other words, herein, rational design of binding peptides towards HLA-A*03-PLP( 45-53 ) complex is studied.
Identi cation of the HLA Hot-Spot Residues for Binding to TCR The lack of available triple molecular crystal structures is a serious obstacle to determine the most important residues of the mentioned p-HLA complex for binding to TCR. Initially, sequence similarity of HLA-A*03 alleles with other class I, HLA proteins was investigated. Sequence similarity was determined using pBLAST (https://blast.ncbi.nlm.nih.gov/pBlast) and Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) tools. Sequence identity of HLA-A*03 alleles with HLA-A*11 and HLA-A*02 was equal to 97 and 93%, respectively. The used TCR databanks included TCR structural repertoire database (https://tcrdb.ibbr.umd.edu/) [22], VDJdb web browser (https://vdjdb.cdr3.net) [23], STCRDab database (http://opig.stats. ox.ac.uk/webapps/stcrdab/) [24], and immune epitope database and analysis resource (http://www.iedb.org/) [25]. The selected complexes for this study consisted of an MHC molecule (HLA-A*02 or HLA-A*11), a 9-or 10-mer peptide, and a TCR. Both TCR subunits were different in the TCR beta variable gene (TRβV) and TCR alpha variable gene (TRαV). Only 12 structures that were different in CDR1, CDR2, and CDR3 sequences were retained) (Table S1). Intersubunit interactions of p-HLA-TCR selected complexes were analyzed by three different computational approaches. Firstly, Robetta, a fast computational alanine scanning method was used to predict energetically important amino acids in the selected p-HLA-TCR complex (http://robetta.bakerlab.org/ alaninescan) [26]. Secondly, knowledge-based FADE and Contacts2 (KFC2) server was used to predict hotspots by analyzing chemical and physical properties in surrounding p-HLA-TCR PPIs interface (https://mitchellweb.ornl.gov/KFC_Server/index.php) [27]. Finally, the HotRegion database was used to predict the hot-spots by accessible surface area and knowledge-based pair energies of each amino acid residue in the target complex (http://prism.ccbb.ku.edu.tr/hotregion/) [28]. Webbased servers used their characteristic calculations and automatically predicted the hot-spots in inter-subunit interactions of the p-HLA-TCR selected complexes. These computational approaches allowed us to predict possible points of interaction with TCR located on the α 1 and α 2 -helix of α chain of the HLA-A*03 molecule. Interface residues indicated in this step with 9 residues of PLP epitope were selected as binding site for all peptides docking (active residue) in ClusPro and HADDOCK servers.

In-Silico Design of Peptide
The rst peptide for inhibition of HLA-A*03-PLP( 45-53 ) was designed innovatively by the above-described crystal structures interfacial contacts. TCR-HLA interfacial contacts were analyzed by the PISA server [29], using a cut-off value of 4 Å for contacts. Considering length of α 1 and α 2 -helix of HLA-A*03 molecule on which the hot/warm spots were distributed, some helical 19-mer peptides were designed in this research. Before docking, structure of the rst peptide was modeled using the PEP-FOLD server (https://bioserv.rpbs.univ-parisdiderot.fr/services/PEP-FOLD/) [30]. The best 3D model was selected according to the PEP-FOLD server, considering the lowest energy model that indicates peptide stability [31]. The binding a nity obtained from the docking results suggested that these peptides were not suitable and they required to be improved. To this end, several steps of a nity maturation were performed by the BeAtMuSiC server (http://babylone.ulb.ac.be/beatmusic/) [32]. This server employed a mutation strategy, and any residue was substituted with all other 19 amino acids separately. This server uses the known structure of a protein-protein complex to evaluate change in binding a nity between two proteins caused by single mutations in their sequence [33]. Each substitution was accompanied by an evaluation of peptides binding a nity and several physicochemical properties, using ClusPro docking server [34] and ProtParam (https://web.expasy.org/protparam/) web tool [35]. ClusPro assigns a weighted ranking to binding energy of the protein complex. In parallel, molecular docking was separately conducted using the HADDOCK server to validate interaction mode of each peptide-HLA-A*03-PLP( 45-53 ) complex. This server is a powerful docking platform that uses a data-driven approach supporting a wide variety of experimental data [36]. Molecular docking was used for obtaining the most possible binding pose for each designed peptide and determining probable binding site of each peptide into the mentioned p-HLA complex. Finally, at each position, amino acid substitutions leading to an increase in binding a nity and stability of the peptides to the p-HLA complex were distinguished and mutations were subsequently used in combinations to design the secondgeneration peptide libraries. Similarly, subsequent peptide libraries were generated until a peptide that had more binding a nity relative to reference peptide was identi ed.
Several physicochemical properties of peptides including molecular weight, net charge at pH=7, pI, aliphatic index, grand average of hydropathicity index (GRAVY) were measured using ProtParam and PepCalc methods [37]. Generally, a positive GRAVY value implies hydrophobicity while a negative GRAVY value suggests a protein's hydrophilicity degree [38]. Besides, antigenicity of the peptide sequence was predicted using the VaxiJen 2.0 server [39] and the predicted antigenic peptides server (http://imed.med.ucm.es/Tools/antigenic.pl). Finally, 14 proper peptides were selected for investigation via MD simulation (Table 1).

MD Simulations
The AMBER18 MD simulation package [40] was used to perform MD simulations of p-HLA complex bound to each of 14 different peptides ( Table 1). Calculations were conducted using the SANDER GPU-accelerated version [41] with Amber ff14SB force eld [42], in periodic boundary conditions. In all cases, the peptide was neutralized by adding sodium ions and was solved in a truncated octahedron box lled with TIP3P water model molecules [43], which led to formation of a solvent shell with a distance of at least 10 Å between each box face and the solute using the tleap module. About 7,000 water molecules were added in solvation stage. The total number of atoms was about 35,000 for each system. Structural minimization was performed in two steps to eliminate any steric con icts. At rst, water molecules and sodium were relaxed with harmonic constraint potential of 2.0 kcal/molÅ 2 on all atoms of both p-HLA protein and designed peptides. Then, the entire system was subsequently minimized by 5,000 steps of the steepest descent minimization followed by 5,000 steps of conjugate gradient minimization. Langevin thermostat was used [44] with coupling coe cient of 2 ps and force constant of 2 kcal/molÅ 2 on the complex. After minimization, MD simulations were initiated by gradually heating each system in canonical ensemble (NVT) from 0 to 300 K in three steps for 200 ps. Restraint constant of heating process was equal to 100, 10, and 1 kcal/molÅ 2 , respectively. Then, density equilibration of 100 ps was performed on the complex by removing all the restraints in isothermal-isobaric (NPT) ensemble. Berendsen barostat with a pressure-coupling constant of 2 ps was used to maintain pressure at 1 bar [45]. Ultimately, production runs were carried out for 40 ns of MD simulations in the NPT ensemble at constant temperature of 300 K for each complex. During MD simulation, particle mesh Ewald (PME) approach was used to treat long-range electrostatic interactions [46] and the SHAKE algorithm restricted all covalent bond lengths involving hydrogen atoms [47]. The cut-off distance for long-range van der Waals (VDW) interaction term was equal to 10 Å. Time step was set to 2 fs, and trajectory coordinates were recorded every 2 ps. Post-MD analyses were performed, using the CPPTRAJ module in Amber18 software [48]. Visualization of trajectories was performed using VMD software [49].

Free Energy Calculations and Per-Residue Contributions
For investigating interaction of each peptide from energetic perspective, binding free energy was calculated based on trajectories of MD simulations by molecular mechanics/generalized born surface area (MM/GBSA) and molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) methods using mmpbsa.py module [50]. These calculations were carried out following the "single trajectory method", where snapshots of binding partners were extracted from MD trajectories of pHLA-peptide complexes. The single trajectory method neglects energetic contributions due to conformational changes leading to a drastic reduction in statistical uncertainty of free energy components [51]. Binding free energy is a signi cant thermodynamic parameter that provides the detailed information about the peptideprotein interaction. It also offers a clear understanding of binding mechanism including their different energy contributions to molecular identi cation, such as enthalpy and entropy. For further exploring the detailed interaction information of each peptide, free energy decomposition was performed using MM/GBSA method to identify key residues responsible for binding energy. Free energy decomposition calculations were done during the last 10 ns of MD trajectories using a single trajectory approach averaged over 1,003 snapshots with 10 ps intervals. The residues contributing less than -2 kcal/mol were considered to be very signi cant for binding of p-HLA complex with peptide and these residues were designated as hot-spot amino acids.
During simulation process using the CPPTRAJ module, the number of hydrogen bonds formed versus simulation time was calculated to detect complex stability and con rm important active site residues interacting with each peptide. Maximum hydrogen-acceptor distance of 3.5 Å and a minimum hydrogen donor-acceptor angle of 135° were allocated to hydrogen bonds [52][53][54]. The frames adopted for calculation were taken from the entire 40 ns of trajectories at intervals of 2 ps. Furthermore, for exactly determining how hydrogen bonds play dominant roles in maintaining peptide stability in MD simulations, hydrogen bond occupations formed during this period were calculated [53]. The occupations vary from 1 to 100%, and a higher percentage represents a more stable hydrogen bond. Hydrogen bonds were only considered if their occupancies were more than 20%.
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 a nity than other peptides. This relatively high binding a nity exhibited by peptide 14 can be signi cantly 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 uctuation 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 uctuations 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 uctuation observed in Fig. 2(a) will be attributable to α 3 helix. In this study, the β 2 M chain was omitted to decrease the computational cost. The absence of the β 2 M domain in simulation led to high RMSD of the α 3 domain. In class I of HLA molecules, it was found that β 2 M interacts with all three domains of α 1 , α 2 , and α 3 whereas the α 3 domain is mostly responsible for binding to β 2 M [55]. As shown in Figs. S2(a-c), the three mentioned regions were distinguished by RMSD uctuations 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 uctuations at rst 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 oor of HLA, and peptides residues, there was hardly a space for conformational movements and uctuation in MD simulation ( Table 2).

Root Mean Square Fluctuation (RMSF)
For understanding uctuations 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 uctuations (F overall ) or total exibility for the three described regions, peptides, and p-HLA complexes were calculated. Binding peptides led to different F overal 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 F overal value relative to other peptides. In other words, peptide 14 induced the most structural uctuation in p-HLA complex. P-HLA complex in interaction with peptide 11 had the lowest F overall 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 uctuations 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 F overal values of PLP belonged to binding peptides of 7 and 12, respectively. In other words, binding peptide of 12 led to more uctuation 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 F overal of PLP peptide (Residues of 275-283 in Fig. 3(a)). It may re ect 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 oor are primarily in uenced by each peptide binding [58]. Our results showed that the RMSF values of α 1 and α 2 helices and binding groove (α 1 +α 2+ β-strand oor) 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 F overall 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 exibility 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 F overal value of peptides belonged to peptides 7 and 11. Then, these peptides had the least uctuation in binding to p-HLA complex. Also, the most uctuation 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 F overal value of p-HLA complex for complex 11 (Table 2). In other words, compaction can lead to low uctuation. 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 F overal 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 F overal 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).  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 e ciency 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 exible 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 identi cation 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 (ΔG Eff ) 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 (∆E ele +∆G GB ) and Van (Fig. 6). Electrostatic interaction and polar solvation energy (∆G GB ) 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 (∆E vdw ) and non-polar solvation energy (∆G SA ) (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 exibility (F overal ) 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 Hbond was de ned 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,  For studying origin of stronger binding a nity 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 ve 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 con rmed 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. 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 con rming low F overal value of PLP in binding of peptide 7 to p-HLA complex (Table 2) (Table 6). HLA class I supertypes have been classi ed 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. Speci city 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 identi ed 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 speci city is determined by electrostatic interactions and polar part of solvation free energy, as con rmed 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).

Conclusion
Pathophysiology of MS is extremely complicated including several immune system branches, and it has yet to be fully understood. Current MS therapies are often non-speci c, resulting in inhibition of the general immune response to combat against pathogenic infections. As a result, more antigen-speci c therapies are needed to avoid this general suppression. In this study, a therapeutic peptide design strategy was proposed to interrupt the TCR-pHLA interaction. Results of the post-MD analyses revealed that the designed peptide 14 was highly capable of binding to the TCR-binding sites of p-HLA complex in-silico, in a selective manner. This peptide could induce a local conformational change in p-HLA complex structure that led to an increase in total hydrogen bond occupancy and favorable protein-peptide interaction. It is hoped that ndings of this research introduce a new strategy for treatment of MS disease. Abbreviations MS; multiple sclerosis, p-HLA; HLA-A*03-PLP( 45-53 ) complex, PPI; protein-protein interaction; MD; molecular dynamics, MHC; major histocompatibility complex, SASA; solvent accessible surface area, RMSD; root mean square deviation, Rg; radius of gyration Declarations Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and material
All data generated or analyzed during this study are included in this published article and its supplementary information les.

Competing interests
The authors declare that they have no competing interests.   RMSD of the P-HLA complex (2a) and designed peptide (2b).

Figure 3
Per-residue backbone RMSF atoms for p-HLA complex (3a) and designed peptide (3b) during the last 10 ns of MD simulation  The number of hydrogen bonds during MD simulation between p-HLA complex and peptides.