Athanassios Stavrakoudis

Ioannis G Tsoulos

Katalin Uray

Ferenc Hudecz

Vasso Apostolopoulos

Homology modeling and molecular dynamics simulations of MUC1-9/H-2K b complex suggest novel binding interactions

Keywords: Class I MHC, H-2K b, homology modeling, molecular dynamics, MUC1, tumor

Human MUC1 is over-expressed on human adenocarcinomas and has been used as a target for immunotherapy studies. The 9-mer MUC1-9 peptide has been identified as one of the peptides which binds to murine MHC class I H-2K b The structure of MUC1-9 in complex with H-2K b has been modeled and simulated with classical molecular dynamics, based on the x-ray structure of the SEV9 peptide/H-2K b complex. Two independent trajectories with the solvated complex (10 ns in length) were produced. Approximately 12 hydrogen bonds were identified during both trajectories to contribute to peptide/MHC complexation, as well as 1-2 water mediated hydrogen bonds. Stability of complex was also confirmed by buried surface area analysis, although the corresponding values were about 20% lower than those of the original x-ray structure. Interestingly, a bulged conformation of the peptide's central region, partially characterized as a β-turn, was found exposed form the binding groove. In addition, P1 and P9 residues remained bound in the A and F binding pockets, even though there was suggestion that P9 was more flexible. The complex lacked the numerous water mediated hydrogen bonds that were present in the reference peptide x-ray structure. Moreover, local displacements of residues Asp4, Thr5 and Pro9 resulted in loss of some key interactions with the MHC molecule. This might explain the reduced affinity of the MUC1-9 peptide, relatively to SEV9, for the MHC class I H-2K b .

Introduction

Major histocompatibility complex (MHC) proteins bind small peptide fragments derived from pathogenic proteins and form peptide/MHC (pMHC) complexes [START_REF] Raghavan | MHC class I assembly: out and about[END_REF]. MHC proteins are divided into two classes: class I (MHC-I) and class II (MHC-II). The MHC-I consists of a polymorphic transmembrane heavy chain and β2-microglobulin, which are non-covalently associated [START_REF] Zhang | Structural principles that govern the peptidebinding motifs of class I MHC molecules[END_REF]. The proteolysis of intracellular proteins by the proteasome produces the majority of peptides suitable for MHC-I binding. In most cases, peptides of 8-10 residues in length are found in the binding groove of MHC-I.

After the first crystal structures of pMHC complexes were available, [START_REF] Bjorkman | Structure of the human class I histocombatibility antigen, HLA-A2[END_REF][START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF][START_REF] Matsamura | Emerging principles for the recognition of peptide antigens by MHC class I molecules[END_REF] it was suggested that peptides bound to MHC-1 with a canonical extended strucure. MHC class I residues that form the binding groove are responsible for the specificity of the peptide selection. Six (out of 8-10) residues of the peptide sequence are accommodated within the A-F binding pockets of the MHC-I protein [START_REF] Saper | Refined Structure of the Human Histocompatibility Antigen HLA-A2 at 2.6 Å Resolution[END_REF]. Residues that do not participate directly in binding are believed to interact with the TCR.

Human mucin, MUC1, is a membrane-bound glycoprotein, expressed on the surface of epithelial cells. It is often overproduced and/or underglycosylated in adenocarcinomas (breast, ovary, colon, lung, kidney, etc) and is present in the serum of cancer patients. MUC1 is immunogenic in mice and in humans, with both humoral and cellular immune responses being induced by MUC1-based vaccine constructs (Tang et al., 2008b,a). MUC1 mucin partly consists of a variable number of tandem repeats region of the consensus sequence 1 PDTRPAPGSTAPPAHGVTSA 20 which is repeated 40-80 times [START_REF] Gendler | A highly immunogenic region of a human polymorphic epithelial mucin expressed by carcinomas is made up of tandem repeats[END_REF]. The majority of anti-MUC1 antibodies recognize sequences within the SA 1 PDTRPAP 7 region [START_REF] Price | Immunological and structural features of the protein core of human polymorphic epithelial mucin[END_REF][START_REF] Xing | Effect of variaions in peptide sequence on anti-human milk fat globule membrane antibody reactions[END_REF][START_REF] Xing | Epitope mapping of anti-breast and anti-ovarian mucin monoclonal antibodies[END_REF][START_REF] Burchell | A short sequence, within the amino acid tandem repeat of a cancer-associated mucin, contains immunodominant epitopes[END_REF]. The SAPDTRPAP (MUC1-9) 9-mer peptide was also found to be presented by MHC-1 H-2K b and to be immunogenic [START_REF] Apostolopoulos | MUC1 peptide epitopes associated with five different H-2 class I molecules[END_REF]. MUC1-9 binds with low affinity to H-2K b [START_REF] Apostolopoulos | MUC1 peptide epitopes associated with five different H-2 class I molecules[END_REF] via a noncanonical mode and it was suggested that the C-terminus of the peptide looped out of the peptide binding groove [START_REF] Apostolopoulos | Anti-MUC1 antibodies react directly with MUC1 peptides presented by class I H2 and HLA molecules[END_REF][START_REF] Apostolopoulos | Noncanonical peptides in complex with MHC class I[END_REF].

Computer simulation of molecular dynamics is a well established method for studying several aspects of biomolecular structure and function [START_REF] Hansson | Molecular dynamics simulations[END_REF][START_REF] Karplus | Molecular dynamics of biological macromolecules: a brief history and perspective[END_REF][START_REF] Aksimentiev | Computer modeling in biotechnology: a partner in development[END_REF][START_REF] Tantar | Docking and Biomolecular Simulations on Computer Grids: Status and Trends[END_REF]. In recent years such computational approaches have been increasingly incorporated in drug design [START_REF] Galeazzi | Molecular Dynamics as a Tool in Rational Drug Design: Current Status and Some Major Applications[END_REF], in immunological reasearch [START_REF] Morikis | Physical methods for structure, dynamics and binding in immunological research[END_REF][START_REF] Mallik | Applications of Molecular Dynamics Simulations in Immunology: A Useful Computational Method in Aiding Vaccine Design[END_REF]Stavrakoudis, 2010b) and also to peptide/MHC complexes [START_REF] Omasits | Analysis of key parameters for molecular dynamics of pMHC molecules[END_REF][START_REF] Knapp | A critical cross-validation of high throughput structural binding prediction methods for pMHC[END_REF]. Moreover, biomolecular modeling can complement experimental studies [START_REF] Van Gunsteren | Molecular simulation as an aid to experimentalists[END_REF] and can elucidate dynamics of immunological synapse [START_REF] Wan | Toward an atomistic understanding of the immune synapse: large-scale molecular dynamics simulation of a membrane-embedded TCR-pMHC-CD4 complex[END_REF], allows to study the dynamics of a peptide bound to antibody [START_REF] Tatsis | Molecular Dynamics Simulations of the TSSPSAD Peptide Antigen in Free and Bound with CAMPATH-1H Fab Antibody States: The Importance of the β-Turn Conformation[END_REF]Stavrakoudis, 2010a), could be used to model disulphide peptide complexed proteins (i.e. C8γ [START_REF] Stavrakoudis | A disulfide linked model of the complement protein C8γ complexed with C8α indel peptide[END_REF]) or even more excitingly to help in clinical decision making [START_REF] Sadiq | Patient-specific simulation as a basis for clinical decision-making[END_REF].

Modeling of the MUC1-9 peptide with both murine and human MHC class I, H-2K b and HLA-A2 respectively have been previously performed [START_REF] Apostolopoulos | Anti-MUC1 antibodies react directly with MUC1 peptides presented by class I H2 and HLA molecules[END_REF], based on a simulated annealing protocol and high temperature molecular dynamics [START_REF] Chelvanayagam | Structural comparison of major histocompatibility complex class i molecules and homology modelling of five distinct human leukocyte antigen-a alleles[END_REF]. That work was a considerable progress in our knowledge of peptide/MHC interactions in the MUC1-9 case and provided a possible structural explanation of the antibody binding of MUC1 peptides presented by the MHC molecules. However, modern progress in computational biophysics, accompanied with the big enhancement of available computer power, can be utilized to further improve the computer-generated model of the MUC1-9 peptide complexed the the MHC class I H-2K b .

Here, we present a homology modeling and molecular dynamics approach of MUC1-9 (SAPDTRPAP) in complex with MHC class I H-2K b . Since the initial conformation was modeled rather than taken from an x-ray structure, we chose to perform two indepent simulation runs, to obtain more robust results. Longrun dynamics, inclusion of the whole MHC molecule and explicit representation of solvent have been utilized in order to more accurate picture the MUC1-9 structure and interactions with the MHC molecule. Such approach has been suggested to give more reliable results in MD investigations [START_REF] Van Gunsteren | Molecular simulation as an aid to experimentalists[END_REF][START_REF] Omasits | Analysis of key parameters for molecular dynamics of pMHC molecules[END_REF]. Our results suggest that this was a beneficial approach in the current study, and has given insights into the peptide binding mode of the MUC1-9.

Methods

Initial coordinates for the SEV9/MHC complex were downloaded from Protein Data Bank [START_REF] Berman | The Protein Data Bank[END_REF], access code: 1kpv. The original peptide from sendaivirus, FAPGNYPAL was mutated to SAPDTRPAP, whilst MHC molecule remained untouched. The SEV9 peptide was selected from other candidates due to its homology with the MUC1-9 peptide. Pro residue homology in positions P3 and P7 was also crucial for selection. Since the backbone dihedral angle φ of Pro residue is restrained, it is preferable to choose a peptide that has the same residue in these positions. Ideally, it would be perfect to also have alignment for position P9, however there was no such option. Topology and force field parameters for all atoms were assigned from the CHARMM22-CMAP parameter set (MacKerell Jr et al., 2004b,a).

It has been noted that addition of cross terms with CMAP potential improves the system parametrization and helps to avoid undesired backbone helical transitions [START_REF] Buck | Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme[END_REF][START_REF] Stavrakoudis | Molecular dynamics simulations of an apoliprotein derived peptide[END_REF].

Hydrogen atoms were added with the VMD program [START_REF] Humphrey | VMD: Visual Molecular Dynamics[END_REF] and its autopsf utility. Protonation status of Histidine side chains were determined with the REDUCE program [START_REF] Word | Asparagine and glutamine: using hydrogen atom contacts in the choice of side-chain amide orientation[END_REF]. The peptide/MHC complex was centered 3 in a rectangular box with dimensions 95.7 ×88.3 ×102.9

. The box was filled with TIP3P water molecules and neutralized with the addition of 26 Na + and 20 Cl -ions respectively, to approximate a 0.1 mM ion concentration. Crystallographic water molecules (345) were also included in the model. The final system contained 24429 water molecules. Total number of atoms of the entire system were 80598.

Non-bonded van der Waals interactions were gradually turned off at a distance between 12 and 14 Å [START_REF] Yonetani | Liquid water simulation: a critical examination of cutoff length[END_REF]. Long range electrostatics were calculated with the PME method [START_REF] Darden | Particle mesh ewald: An N log(N ) method for ewald sums in large systems[END_REF]. Non-bonded forces and PME electrostatics were computed every second step. Pair list was updated every 10 steps.

Bonds to hydrogen atoms were constrained with the SHAKE method allowing a 2 fs time step for integration. The system was initially subjected to energy minimization with 5 000 steps. The temperature of the system was then gradually increased to 310 K, with Langevin dynamics using the NVT ensemble, during a period of 3 000 steps, by stepwise reassignment of velocities every 500 steps. The simulation was continued at 310 K for 100 000 steps (200 ps). During minimization and equilibration phases, protein backbone atoms (N, C α , C', O) and oxygen atoms of crystallographic waters were restrained to their initial positions with a force constant of 50 kcalmol -1 Å-2 .

The system was equilibrated for further 200 ps with the force constant reduced to 5 kcalmol -1 Å-2 . Finally, 400 ps of NVT simulation at 310 K was performed with total elimination of the positional restraints. The simulation was passed to the productive phase, by applying constant pressure with the Langevin piston method [START_REF] Feller | Constant pressure molecular dynamics simulation: The Langevin piston method[END_REF]. Velocities were re-initialized and two independent trajectories were produced (trA and trB). Pressure was maintained at 1 atm and temperature at 310 K. Results are based to a period of 10 ns of this isothermal-isobaric (NPT) runs. Snapshots were saved to disk at 1 ps interval for structural analysis.

The initial structure of the SEV9/MHC complex (PDB code 1kpv) were also simulated under identical conditions for comparative analysis (tr0 trajectory).

Trajectory analysis was performed with Eucb [START_REF] Tsoulos | Eucb: a C++ program for molecular dynamics trajectory analysis[END_REF][START_REF] Stavrakoudis | A disulfide linked model of the complement protein C8γ complexed with C8α indel peptide[END_REF] and Carma [START_REF] Glykos | Carma: a molecular dynamics analysis program[END_REF] software packages. Secondary structure analysis was performed with STRIDE [START_REF] Frishman | Knowledge-based protein secondary structure assignment[END_REF]. Circular l33 l34 l35 data statistics (dihedral angles, etc) were calculated with appropriate corrections [START_REF] Agostinelli | Circular Statistics with R[END_REF]Do ¨ker et al., 1999) β-turn classifications were based on geometrical characteristics of the backbone conformation [START_REF] Hutchinson | A revised set of potentials for beta-turn formation in proteins[END_REF]. Initially, a β-turn was accepted water molecules in the peptide/MHC interface the instantaneous water coordination number (N c ) approach [START_REF] Petrone | MHC-peptide binding is assisted by bound water molecules[END_REF]. This method counts the water oxygen atoms within a range (typically 3.5 A) of any water oxygen atom, which is actually the first hydration shell. The N c can be found between 0 and 15, depending on the local structure of water. In the bulk water this number is always greater than 3, while in the protein interior l35 l57 l58 l59 l60 is 0 to 2. This implies that a water molecule has no other water neighbours and it is inside the protein interior. The N c is measured for all the MD trajectory and isolated water molecules are identified if the N c value is small for a prolonged period of time. In the current study, a search of water molecule with N c < 1 for at least 70% of the MD time has been performed. l6l

l45 if d(G a -G a ) < 7 Å and |a G a -G a -G a -G a | < 90 o ,

MM-PBSA calculation of ∆G binding

The binding free energy of the association of two molecules (A+B→AB) can be estimated, according to the MM-PBSA approach [START_REF] Kollman | Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models[END_REF][START_REF] Wan | Peptide recognition by the T cell receptor: comparison of binding free energies from thermodynamic integration, Poisson-Boltzmann and linear interaction energy approximations[END_REF], as:

∆G binding = G AB -G A -G B , (2)
where: In the above equations, h.i denotes average value for a a set of snapshots along a molecular dynamics trajectory, while E M M is the molecular mechanics energy of the i th molecule in the gas phase, namely the sum of f internal bonded energy (comprising bond, angle and dihedral terms), van der Waals and electrostatic interactions. G solv is the solvation free energy of the i th molecule. This term can be estimated as the sum of the electrostatic solvation free energy calculated by the Poisson-Boltzmann equation and the non-polar solvation free energy calculated from the SASA.

∆G i = hE M M i + hG solv i -T S. (3)
Hence, the binding free energy is:

∆G i = h∆E M M i + h∆G solv i -T S. (4)
The average properties can be computed directly from the MD trajectory snapshots. In the current study, the last 5 ns were used, assuming that equilibrium was reached after the first 5 ns of the simulation. 5000 structures were utilized for the l72 solv

SASA and E M M calculations, while 50 structures (one every 100 frames) were used l73 for the calculation of the G elec with the APBS [START_REF] Baker | Electrostatics of nanosystems: application to microtubules and the ribosome[END_REF]Dolinsky et al., l74 2004) software. we take into consideration the simulation temperature (310 K) these values are considered small, indicating the stability of the complex. Moreover, there is strong evidence that the MHC molecule did not undertake significant conformational changes upon mutation of the peptide residues [START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF][START_REF] Matsamura | Emerging principles for the recognition of peptide antigens by MHC class I molecules[END_REF]. This is in accordance with other X-ray studies of the H-2K b MHC class-I molecule with different nonamer peptides in the binding groove. These observations corroborate our hypothesis that homology modeling coupled with molecular dynamics simulations produces a reliable model of the MUC1-9/H-2K b complex.

Peptide's RMSF values of C a atoms showed an interesting differentiation between 8 trA and trB trajectories. While values of 0.5-1.0 Å were recorded for residues 1-7 in both cases, trA trajectory showed increased values of 1.5 and 2.0 Å for residues 8 and 9 respectively. In trB trajectory, RMSF remained close to 1.0 Å for all residues. Values in the order of 2.0 Å are still considered relatively small, however, the differentiation is notable. Since this fact was observed in only one of the two trajectories, it could be considered as a relative random effect of the simulation.

On the other hand, it definitely indicates that the peptide binding to the MHC groove is not so tight at the C-terminal region, as previously has been suggested [START_REF] Apostolopoulos | Anti-MUC1 antibodies react directly with MUC1 peptides presented by class I H2 and HLA molecules[END_REF][START_REF] Apostolopoulos | Noncanonical peptides in complex with MHC class I[END_REF]).

Peptide's RMSD time series of backbone atoms were very similar in both cases.

RMSD values ranged between 0.99 and 2.1 Å and averaged at 1.46 (0.16) Å for trA, whilst the RMSD values ranged between 1.04 and 1.69 Å and averaged at 1.36 (0.09) for trB case. There is only a minor difference between these two profiles: Backbone conformations play an important role in peptide/MHC binding [START_REF] Barinaga | Getting some " backbone" : how MHC binds peptides[END_REF][START_REF] Matsamura | Emerging principles for the recognition of peptide antigens by MHC class I molecules[END_REF]. Here we present a detailed analysis of the peptide's backbone conformation. , 2004). However, this happens to the peptide's region that is outside of the binding group. In the current study, we have identified a very interesting case of β-turn in the central region of the peptide, covering residues Pro3 to Arg6. This sequence has been found in β-turn conformation for 50 and 77% of the simulation time, in the trA and trB trajectories respectively. We did not observed any intra-peptide hydrogen bond stabilizing this β-turn. Table 1 lists the values of backbone dihedral angles as calculated for the central residues of the β-turn, Asp4 and Thr5 respectively.

Both trA and trB trajectories showed very close values of backbone φ and ψ dihedral angles. These values differ from the initial values found in the crystal structure of the SEV9 peptide. However, the common finding from the two independent trajectories (trA and trB) corroborate the suggestion that a β-turn around the Asp4-Thr5

region exists, at least partially.

Bulged peptide structure have been identified in previous studies [START_REF] Speir | Two different, highly exposed, bulged structures for an unusually long peptide bound to rat MHC class I RT1-Aa[END_REF][START_REF] Tynan | T cell receptor recognition of a' super-bulged' major histocompatibility complex class I-bound peptide[END_REF], however for longer peptide sequences, like the ILF-PSSERLISNR peptide sequence [START_REF] Speir | Two different, highly exposed, bulged structures for an unusually long peptide bound to rat MHC class I RT1-Aa[END_REF]. The presence of a bulged (turn)

1 0 1 l101 28l
structure ln the nonamer peptldes bound ln the class I MHC molecules has been suggested to profoundly affect the recognltlon of the pMHC complex by the T-cell receptor [START_REF] Rudolph | How TCRs bind MHCs, peptides, and coreceptors[END_REF].

Interactions between the peptide and the MHC 29l 30l

The blndlng mode of nonamer peptldes wlth the H-2K et al., 1992;[START_REF] Saper | Refined Structure of the Human Histocompatibility Antigen HLA-A2 at 2.6 Å Resolution[END_REF]. Resldues P4 and P5 do not make dlrect contacts wlth the MHC molecule and protrude towards the solvent, hence thelr slde chalns are avallable for lnteractlon wlth the TCR. The charge groups of N-and C-termlnal resldues make strong lnteractlons wlth the MHC blndlng clefts (pockets A and F respectlvely).

A general vlew of the peptlde/MHC blndlng motlf ln shown ln Flgure 3, whllst the peptlde's orlentatlon lnslde the MHC's blndlng ls deplcted at Flgure 4.

Peptlde's Ser1 (P1) was found to form two stable hydrogen bonds wlth the MHC molecule. Its backbone atoms N and 0 were found ln hydrogen bond state wlth slde chalns of Glu63A and Tyr159A respectlvely. These hydrogen bonds were conserved, ln both trA and trB trajectorles, for approxlmately 91 to 95% of the slmulatlon tlme (Table 2). The dlstance between Ser1:N and Glu63A slde chaln oxygen atoms, ln the lnltlal structure, were found 4.6 and 5.8 Å for 0 c1 and 0 c2 respectlvely, whlch lndlcates that thls strong (charged) hydrogen bond between the N-termlnal group of the peptlde and the slde chaln of Glu63A was formed durlng the modellng proccess and was not present ln the lnltlal structure. Indeed, Glu63A's 8 slde chaln (atom 0 c1 ) actually was to form a hydrogen bond wlth Ala2:N atom, ln the structure of the orlglnal peptlde [START_REF] Matsamura | Emerging principles for the recognition of peptide antigens by MHC class I molecules[END_REF]. The hydrogen bond between Ser1:0 and Tyr159A:0 η , on the other hand, was well formed ln the lnltlal structure (dlstance 2.67 A) and very well conserved ln both MD trajectorles (Table 2). Another hydrogen bond lnteractlon between Ser1 and the MHC molecule was present between the slde chalns of Ser1 and Tyr7A (or Tyr171A for short perlods), for approxlmately 95% of the slmulatlon tlme. Thls ls very lnterestlng, slnce no slde-chaln lnteractlons have been observed ln the x-ray structure of SEV9 peptlde [START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF]. Thus, overall two to three hydrogen bonds contrlbuted to peptlde's blndlng. These results corroborate the lmportance of thls blndlng pocket ln the peptlde/MHC blndlng process.

Slde chaln of Glu63A (pocket B) accepted hydrogen bond from Ala2 Nltrogen atom (posltlon P2). Thls lnteractlon was conserved for 93.5% (trA) or 98.6% (trB) of the slmulatlon tlme, and lt was well formed ln the lnltlal structure (the dlstance between Ala2:N and Glu63A:0 δ1 was found 2.9 A). Thls findlng underllnes the lmportance of the Glu63A resldue, slnce lts negatlvely charged slde chaln formed two stable hydrogen bonds wlth the peptlde's backbone amlde groups. Slde chaln of Lys66A was found ln hydrogen bond state wlth Ala2:0 atom for over 90% of the slmulatlon tlme. The correspondlng dlstance between Lys66A:N ζ and Ala2:0 atoms ln the lnltlal structure was found 2.7 A, lndlcatlng the exlstense of the hydrogen bond. Moreover, slde chalns of Tyr7A and Tyr45A made hydrophoblc contacts wlth Ala2's allphatlc slde chaln. The above analysls ls for the Ala2 lnteractlons ls almost ldentlcal wlth the x-ray structure of the SEV9 peptlde [START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF], lndlcatlng the fact the preservatlon of the Ala2 resldue ln posltlon P2 (blndlng pocket B) contrlbuted to the retalng of the same peptlde/MHC lnteractlons.

Pro3 (P3) made lmportant hydrophoblc lnteractlons wlth Tyr159A's slde chaln.

Average dlstance of thelr slde chaln centers were found 4.0 A(0.6) or 4.2 A(0.6) durlng trA or trB MD trajectorles respectlvely. For approxlmately 25% of the tlme, 34l 35l 36l the two slde chalns were found ln parallel orlentatlon formlng a stacklng lnteractlon.

It ls noted that Tyr159A's slde chaln donated a hydrogen bond to Ser1:0, hence thls MHC resldue ls consldered to contrlbute slgnlficantly to peptlde's blndlng. The orlglnal hydrogen bond between Pro's backbone oxygen atom and Asn70A's slde chaln was found to be relatlvely weak durlng trA and trB MD trajectorles: 12.7 and 27.2% of the frames respectlvely satlsfied the hydrogen bond crlterla.

Central resldues Asp4 and Thr5 dld not show any slgnlficant lnteractlons wlth the MHC's resldues. 0nly Asp4's slde chaln was found hydrogen bonded to Arg62A's slde chaln for llmlted perlod of slmulatlon tlme, ≈ 15%. Both resldues were exposed outslde of the blndlng groove.

Blndlng pocket C plays an lmportant role ln peptlde recognltlon by MHC H-2K b molecules [START_REF] Molano | Peptide selection by an MHC H-2Kb class I molecule devoid of the central anchor (" C" ) pocket[END_REF]D. H. Fremont and E. A. Stura and M. Matsumura and P. A. Peterson and I. A. Wllson, 1995;[START_REF] Huard | The critical role of a solvent-exposed residue of an MHC class I-restricted peptide in MHC-peptide binding[END_REF]). Peptlde's resldue

Arg6 slde chaln at posltlon P6, was found to form a strong hydrogen bond wlth Glu24A slde chaln. Actually, these slde chalns remalned hydrogen bonded the entlre tlme tlme ln both trA and trB MD trajectorles. 0n the other slde, there was no backbone lnteractlon wlth the MHC molecule. However, the ablllty of the MHC molecule to blnd dlfferent peptlde sequences, slnce the orlglnal peptlde has Tyr ln thls posltlon [START_REF] Apostolopoulos | Crystal structure of a non-canonical high affinity peptide complexed with MHC class I: a novel use of alternative anchors[END_REF], whlch ls a canonlcal resldue at thls posltlon for MHC blndlng. Tyr6 (SEV9 peptlde) to Arg6 (MUC1-9 peptlde) mutatlon led to some loss of hydrophoblc lnteractlons between peptlde and MHC molecule, a fact that mlght explaln the reduced blndlng affinlty of the MUC1-9 peptlde, relatlve to SEV9 peptlde. However, the Arg6 remalned lnslde the canonlcal C-pocket, unllke the Arg6 resldue ln YEA9 peptlde (SRDNSRIPM) whlch utlllzed the non-canonlcal E blndlng pocket [START_REF] Apostolopoulos | Crystal structure of a non-canonical high affinity peptide complexed with MHC class I: a novel use of alternative anchors[END_REF].

Resldue Pro7, at peptlde's P7 posltlon, had a weak backbone hydrogen bond wlth Tyr117A's slde chaln. 0ccurence was found 28% ln trA and only 7% ln trB trajectorles respectlvely. Glven the fact that ln crystal structures of peptldes bound 37l 38l 39l in the H-2K b molecule, no such hydrogen bond exist (Table 2), the result is not so suprising. However, significant hydrophobic interactions with Trp147A and Trp133A side chains were found to contribute in peptide/MHC interactions. For example, side chain distances between Pro7 and Trp147A varied between 3 and 5 Å and averaged at 3.6 Å (0.2). To a lesser degree, Leu156A and Tyr116A also made hydrophobic contacts with side the chain of Pro7.

Position P8 was occupied by Ala8. The backbone carbonyl group of this residue was found to be in hydrogen bond state with Trp147A's side chain. This is a well expected interaction, as it has been found in the crystal structure of the original peptide. A relatively weak hydrogen bond was also formed for part of trA trajectory, between Ala8:N and Glu152A:Oc2. The corresponding distance in the initial structure was found to be 5.8 A.

Finally, residue Pro9 at position P9 (binding pocket F). The C-terminal carboxyl group was found to form two hydrogen bonds (Table 2) with Thr143A and Lys146A side chains, for almost all of the simulation time, in both trA and trB trajecoties. The same interactions were also present in the x-ray sructure that served as initial point for these calculations. However, the lack of amide hydrogen in proline's structure resulted to the abolishment of a backbone hydrogen bond between peptide and the MHC molecule. Thus, the Leu to Pro (SEV9 to MUC1-9 peptide) mutation resulted in a small shift of the position of this residue. These subtle changes in peptide's conformation have been shown [START_REF] Hoare | Subtle Changes in Peptide Conformation Profoundly Affect Recognition of the Non-Classical MHC Class I Molecule HLA-E by the CD94-NKG2 Natural Killer Cell Receptors[END_REF] to affect drastically the MHC recognition and might explain to some extent the reduced affinity of the MUC1-9 peptide when bound to class I H-2K b . Pro9's side chain also made hydrophobic contacts with Val76A, Leu81A and Trp147A side chains. For at least 90% of the simulation time, a pair of side chain heavy atoms from these residues were in close contact (distance less than 4.5 A) with a side chain heavy atom from Pro9.

These hydrophobic interactions further stabilized the peptide/MHC interactions, and along with the hydrogen bonds strengthen the anchoring role of Pro9. Overall, as it can be seen from Figure 6, there were approximately 12 hydrogen bonds between the peptide and the MHC molecule, during both MD trajectories.

This number approximates very well the number of the reported [START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF] hydrogen bonds (11) between the peptide SEV9 and the MHC molecule. Buried surface area (BSA) is a good indicator of the binding of a ligand into a protein [START_REF] Olsson | The thermodynamics of protein-ligand interaction and solvation: Insights for ligand design[END_REF]. Figure 7 shows the time evolution of BSA between the peptide and the MHC molecule. BSA fluctuated between 666.9 and 1005.6 Å2 400 and averaged at 848.7(47.5) 2 in the trA trajectory. In the trB case, BSA values 40l were found between 656.7 and 999.6 2 with mean value of 824.9(52.9) 2 . As it were expected as a measurement error rather than actual difference in BSA. These findings corroborate our statement that the peptide/MHC complex was stable and that the fluctuations in BSA time series are normal.

The BSA value in the x-ray structure of the SEV9 peptide was 1076 Å2 , while 4l2 the BSA value in the MUC1-9/MHC complex after restrained energy minimization 4l3 was found to be 937.5 2 . The loss of approximately 140 2 can be attributed to A main difference between MUC1-9's simulated structure and SEV9's x-ray structure was the φ angle of Thr5. A significant transition from +74 o to ≈ -150 o occurred. Since it is well known that backbone conformation plays a very important role in peptide/MHC recognition [START_REF] Barinaga | Getting some " backbone" : how MHC binds peptides[END_REF], it is expected that this conformational transition would alter the MHC's binding affinity for the peptide, most possibly downwards. Moreover, MHC H-2K b molecules prefer hydrophobic residues at position P6 (for nonamer peptides), even though MUC1-9 has Arg in this place.

This has resulted in a notable alteration of the backbone conformation of the central part of peptide and the enhancement of the exposure of the Asp4-Thr5 re-4 5 45l gion outside of the MHC's binding groove. For a considerable amount of simulation time this bulged region [START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF]) adopted a β-turn conformation, however without the presence of the characteristic hydrogen bond. This had not been noted in previous modeling studies [START_REF] Apostolopoulos | Anti-MUC1 antibodies react directly with MUC1 peptides presented by class I H2 and HLA molecules[END_REF] and provides a new framework for the peptide/MHC interactions.

Inclusion of explicit water molecules in the current study helped a lot to clarify the role of the solvent in peptide/MHC interactions. Water mediated hydrogen were found only sparingly and although existed, a clear contribution to the binding process can not be attributed to this kind of interaction.

Leu to Pro mutation at position P9 resulted in slight movement of this residue within the F binding pocket. However, this fact, along with the loss of a hydrogen bond interaction of the Leu amide hydrogen might be enough reason for observing the reduced affinity of the MUC1-9 peptide to H-2K b binding.

  fluctuations (RMSF) of the C a atoms of the MHC and peptide chains, as well as the time evolution of the root mean square deviation (RMSD) of the backbone atoms (N, C a , C') of the MHC and peptide chains, during both MD trajectories, trA and trB respectively, are shown in Figure 1.In both trA and trB cases, RMSF profiles of chains A and B from the MHC molecule were almost identical, which indicates the robustness of the study. RMSF values were between 0.5 and 2.0 A, which is quite common in similar MD studies of protein complexes around equilibrium. Similarly, RMSD time series were also very similar for chain A and B, with only a small exception of the trA trajectory: RMSD values escaped from stationarity around 4ns in trA, and a small peak of RMSD 2.8 Å was observed. In general, both trajectories were quite stable, Figure1. Time series of RMSD fluctuated around 1.5-2.0 Å for chain A and around 1.0 Å for chain B. If

  trB trajectory showed slightly smaller values with smaller standard deviation of the time series. This is possibly due to increased fluctuation at the C-terminal end in trA . However, as it was previously noted for chains A and B of the MHC molecule, RMSD profiles of the peptide corroborate the stability of the trajectories and the validity of the proposed model.

Figure 2

 2 Figure 2 displays the distribution (Ramachandran map) of the backbone dihedral angles φ, ψ of peptide residues in the region 2-8. It is evident that, for most of the residues, the backbone dihedral angles show very similar distributions in the trA and trB trajectories. The only exception comes from the Ala8 residue. As it has been noted, the C-terminal residues showed increased mobility (higher RMSF values),

  drawn from the graphical representation of BSA time evolution, and from basic statistical analysis, both trajectories showed similar profiles for the calculated BSA of peptide/MHC interface. The difference of approximately 25 Å2 (3%) in the 405 406 mean values is very small and could be considered to be within expected error. In a recent experimental re-investigation of BSA of protein x-ray structures (

  that occurred during MD run in order for the mutated peptide to adapt to the binding groove of the MHC molecule. Considered, however, that trajectories were obtained in 310 K, thus the spontaneous thermal moving of the atoms resulted in somewhat reduced BSA values. 2 4l8 Thus, the difference of more than 140 in the BSA of the peptide/MHC in-the SEV9 and MUC1-9 cases, is another indication of the lower binding affinity that the MUC1-9 has to the H-2K b molecule, relatively to the SEV9 peptide. molecular dynamics simulations have been used to assess the structure of the SAPDTRPAP/H-2K b complex. Results presented here indicate that a stable complex is formed, based on the analysis of two MD trajectories. MHC binding pockets A and F interacted closely with the N-and C-terminus of the peptide which played an important role in stabilizing the complex. The Buried Surface Area of the peptide/H-2K b interface remained constant during the simulation indicating the stability of the complex and its similarity to the initial peptide/MHC complex. Replacement of Leu with Pro at P9 position did not affect significantly the MHC's binding of the peptide. The C-terminal carboxyl group was found to form stable hydrogen bonds with the MHC molecule, and the non-polar side chain of Pro residue made a number of close contacts with hydrophobic residues of the MHC's F binding pocket. However, the peptide showed relatively increased mobility in the C-terminal region, that may affect the strength of the MHC binding.

  Figures

589Figure 1 :

 1 Figure 1: Root mean square fluctuation (left column) of C a atoms and root mean square deviation (right column) time series of backbone atoms (N, C a , C') of the pMHC complex after fitting the corresponding atom positions from MD trajectory to initial (X-ray) coordinates. Results from different trajectories (tr0, trA and trB) are indicated with different line colors. A) RMSF of MHC chain A, B) RMSF of MHC chain B, C) RMSF of MHC chain P (peptide), D) RMSD of MHC chain A, E) RMSD of MHC chain B and F) RMSD of MHC chain P (peptide).

Figure 2 :

 2 Figure2: Ramachandran plot of backbone dihedral angles of the peptide. Horizontal axis is for φ and vertical axis is for ψ angle respectively. The plots represent propability density maps, z-axis is the percentage of frames found within 10 o dihedral angle bin. The adjacent colour bar is used to identify regions of low (grey) versus high (blue) populations.

Figure 3 :Figure 4 :

 34 Figure 3: A) Ribbon representation of five selective structures of the MUC1-9/H-2K b complex (one frame every 2 ns) from trA trajectory, B) Ribbon representation of five selective structres of the MUC1-9/H-2K b complex (one frame every 2 ns) from trB trajectory, C) Stick representation of the peptide bound in the MHC groove from trA trajectory, D) Stick representation of the peptide bound in the MHC groove from trB trajectory, E) Important hydrogen bond interactions between the pepetide and MHC molecule in the trA trajectory and F) Important hydrogen bond interactions between the peptide and MHC molecule in the trB trajectory. Hydrogens were omitted from stick represantations. Structures have been fitted to the first frame using the backbone atoms.

Figure 5 :

 5 Figure5: A) and B) Backbone overlay of the Pro3-Arg6 region of the peptide from the trA and trB trajectories respectively. This fragment has been found in β-turn conformation for considerable amount of time. C) Backbone superimposition of SEV9 peptide (green) from the X-ray structure with representative structures from trA (cyan) and trB (orange) trajectories. The differentation of backbone conformation at fragment Asp4-Thr5 is well seen. Side chains of residues 2, 3, 6 and 7 share common orientation towards the MHC binding groove. Interestingly, conformations of residues at positions 1 and 9 deviate from the original structure.

Figure 6 :

 6 Figure 6: A) Total number of hydrogen bonds between the peptide and the MHC molecule, as evolved over simulation time. Data were averaged every 10 ps. B) Total number of water mediated hydrogen bonds between the pepetide and MHC molecule, as evolved over simulation time. Data were taken every 10 ps.

Figure 7 :

 7 Figure 7: Time series of apolar buried surface area (BSA) between the peptide and the MHC molecule in tr0, trA and trB trajectories.

Figure 1 Figure 6

 16 Figure 1Click here to download Manuscript: rms.eps

  b MHC class I molecule has been lnvestlgated ln the past. There are numerous studles ln the llterature (Mat-

samura et al., 1992;[START_REF] Fremont | Crystal structures of two viral peptides in complex with murine MHC class I H-2K b[END_REF][START_REF] Apostolopoulos | Crystal structure of a non-canonical high affinity peptide complexed with MHC class I: a novel use of alternative anchors[END_REF] Meljers et al., 2005; D. H. Fremont and E. A. Stura and M. Matsumura and P. A. Peterson and I. A. Wllson, 1995)

addresslng the prlnclples of peptlde anchorlng to MHC's blndlng groove. It ls generally assumed that H-2K b has slx blndlng pockets, A to F, that accomodate resldues P1,P2,P3,P6,P7 and P9 of nonamer peptldes

(Matsamura

Table 1

 1 All of the above observations reflected well in the reduction of the BSA between the peptide and the MHC molecule, where a loss of 140 Å2 has been measured. Backbone dihedral angles in the region Asp4-Thr5 of the MUC1-9 peptide, where a β-turn was found in the MUC1-9 peptide. Corresponding region of the SEV9 peptide is Gly4-Asn5. Column PDB lists the corresponding values from the crystal structure of the SEV9 peptide, with the residues Gly and Asn at positions 4 and 5 respectively. Averages values (and variances in parentheses) is fiven from trajectories tr0, trA and trB.

		Finally, it seems that while the MUC1-9 peptide forms stable complex with the H-
	46l	2K b molecule, it is clear that certain structural reorganization occurred and resulted
		in reduced binding affinity.

Table 2

 2 Hydrogen bond interactions between the SEV9 and MUC1-9 peptides and the H-2K b molecule. Percentage of frames is given, from trajectories tr0, trA and trB, that met the geometrical criteria for hydrogen bond interaction. Distance between donor-acceptor atoms are taken from the initial structure (PDB column).

	Donor	Acceptor	PDB (A)	tr0 (%)	trA (%)	trB (%)
	Phe 1P :N	Tyr 59A :O η	4.13	32.4
	Phe 1P :N	Glu 63A :O c 1,2	4.61	95.6
	Ser 1P :N	Glu 63A :O c 1,2	5.83		93.2	94.1
	Tyr 159A :O η	Ser 1P :O	2.67	59.6	91.4	94.8
	Tyr 59A :O η	Ser 1P :O γ	5.77		77.7	81.9
	Ser 1P :O γ	Tyr 7A :O η	5.82		93.2	97.9
	Ser 1P :O γ	Tyr 171A :O η	5.00			20.6
	Ala 2P :N	Glu 63A :O c 1,2	2.90	90.1	93.5	98.6
	Lys 66A :N ζ	Ala 2P :O	2.74	76.6	97.3	90.0
	Asn 70A :N δ	Pro 3P :O	3.63	63.4	12.7	27.2
	Arg 62A :N η 2	Asp 4P :O δ 1,2	6.78		19.0
	Arg 6P :N η 1,2	Glu 24A :O c 1,2	5.21		92.4	97.2
	Tyr 116A :O η	Pro 7P :O	4.11		28.2	7.2
	Ala 8P :N	Glu 152 :O c 1,2	5.87		38.5	5.5
	Trp 147A :N c	Ala 8P :O	2.86	11.7	70.4	98.5
	Leu 9P :N	Asp 77A :O δ 1,2	3.02	82.4
	Lys 146A :N ζ	Leu 9P :O τ 2	3.00	95.2
	Tyr 84A :O η	Leu 9P :O τ 2	2.84	39.3
	Thr 143A :O γ	Pro 9P :O τ 1,2	2.68		94.1	92.7
	Lys 146A :N ζ	Pro 9P :O τ 1,2	3.00		98.6	98.7

l73

Acknowledgements

NAMD parallel execution have been performed at the Research Center of Scientific Simulations (RCSS) of the University of Ioannina. The open source community is gratefully acknowledged for proving all the necessary tools (Linux, NAMD, GNU, etc) that made this work possible.

Figure 7

Click here to download attachment to manuscript: bsa.eps