Molecular Dynamics Simulation Reveals the Selective Binding of Human Leukocyte Antigen Alleles Associated with Behçet's Disease

Behçet’s disease (BD), a multi-organ inflammatory disorder, is associated with the presence of the human leukocyte antigen (HLA) HLA-B*51 allele in many ethnic groups. The possible antigen involvement of the major histocompatibility complex class I chain related gene A transmembrane (MICA-TM) nonapeptide (AAAAAIFVI) has been reported in BD symptomatic patients. This peptide has also been detected in HLA-A*26:01 positive patients. To investigate the link of BD with these two specific HLA alleles, molecular dynamics (MD) simulations were applied on the MICA-TM nonapeptide binding to the two BD-associated HLA alleles in comparison with the two non-BD-associated HLA alleles (B*35:01 and A*11:01). The MD simulations were applied on the four HLA/MICA-TM peptide complexes in aqueous solution. As a result, stabilization for the incoming MICA-TM was found to be predominantly contributed from van der Waals interactions. The P2/P3 residue close to the N-terminal and the P9 residue at the C-terminal of the MICA-TM nonapeptide served as the anchor for the peptide accommodated at the binding groove of the BD associated HLAs. The MM/PBSA free energy calculation predicted a stronger binding of the HLA/peptide complexes for the BD-associated HLA alleles than for the non-BD-associated ones, with a ranked binding strength of B*51:01 > B*35:01 and A*26:01 > A*11:01. Thus, the HLAs associated with BD pathogenesis expose the binding efficiency with the MICA-TM nonapeptide tighter than the non-associated HLA alleles. In addition, the residues 70, 73, 99, 146, 147 and 159 of the two BD-associated HLAs provided the conserved interaction for the MICA-TM peptide binding.


Introduction
Behçet's disease (BD) has caused recurrent inflammation of symptom complex in genital and oral ulceration, inflammatory lesion of central nervous systems, gastrointestinal tract and eyes leading to blindness [1]. This disease is usually reported in adults aged between 18 to 40 yearolds and presents a higher mortality risk in males than females [2]. A high epidemiology of BD is commonly found in Japan, Turkey, China and the Middle East and the Mediterranean countries [3]. The systemic inflammatory disease is characterized by recurrent exacerbations and spontaneous remissions with varying healing times among patients, and the disease course usually becomes several years or more. On the ocular disease, recurrent and often bilateral episodes of acute exacerbations of intraocular inflammation (ocular attacks) in the anterior chamber (iridocyclitis) normally occur with or without posterior involvement (retinal vasculitis, hemorrhages, exudates, retinal vein occlusion and optic neuritis). Recurrent ocular attacks of posterior involvement lead to poor visual prognosis [4]. Intestinal and neurological diseases are serious complications of BD and can severely deteriorate the patients' psychosomatic status. Because BD can induce severe damage in the eyes and other organs, it is important to clarify the pathogenesis of this disease and to develop a novel treatment based on the pathogenesis.
The exact etiology of this disease has not been clarified yet, but it is probably mediated by a combination of excessive immune reactions, genetics factors [5], immune reactions against infectious agents [6], heat shock proteins [7], oxidative stress [8] and environmental factors. Indeed, several inflammatory cytokines involved in acute inflammatory reactions, such as interleukin (IL)-6, tumor necrosis factor-α and IL-8 might play a role in the pathogenesis of this disease [9]. In addition, genetic factors, including the human leukocyte antigen (HLA)-B51 alleles and especially the HLA-B Ã 51:01 allele is strongly associated with BD across different susceptible ethnic groups [10,11]. Recently, a meta-analysis of 78 studies involving 4,800 cases of BD and 16,289 controls reported the pooled odds ratio (OR) of HLA B51/B5 carriers to develop BD compared to non-carriers (OR = 5.78 95%, CI = 5.00-6.67) [12].
The notion of an environmental trigger of BD in patients with genetic susceptibilities has long been advocated. Several infectious agents have been investigated, especially bacteria (Streptococcus, Mycoplasma and Helicobacter pylori) and viruses (Herpes simplex virus 1 and 2, Hepatitis virus and Parvovirus B19) [13]. The autoimmunity against bacterial and human heat-shock proteins (HSP) was also speculated as a trigger for BD, because autoantibodies against HSP65 and HSP60 have been reported in BD patients [14].
On the other hand, the genetic association of the major histocompatibility complex (MHC) class I chain related to gene A (MICA) with BD has been reported in various ethnic groups [15,16], and a strong association of the polymorphism in the transmembrane region of MICA (MICA Ã 009) with BD was observed [17,18]. The self nonapeptide of AAAAAIFVI located on the MICA transmembrane region (MICA-TM), induced autoreactive CD8+ cytotoxic T lymphocytes in BD patients with HLA-B Ã 51 [19]. Accordingly, the MICA-TM peptide is one of the candidate bechetogenic epitopes [17,19], where the self-antigen AAAAAIFVI could function as the antigen that triggers the T-cell receptor and develops to an autoimmune reaction [19][20][21]. The MICA is usually expressed in fibroblasts, monocytes, epithelial cells and endothelial cells in a stress-dependent manner, and plays a key role in initiating the antibody-dependent rejection of organ-transplants [22][23][24].
The HLA molecules are human equivalents of the MHC found in most vertebrates [25]. They are controlled by genes located at the terminal region at band p21.3 on the chromosome 6. Normally, the HLA molecules play a key role as a guard to protect cells in immune recognition. The HLAs are categorized into three different classes (I, II and II), with HLA class I further divided into three main genes(HLA-A,-B and-C) [26]. The HLA alleles associated with BD involve the HLA class I genes HLA-A and HLA-B that importantly function to bind and send short foreign antigens from within the cell to the T-cell receptor (TCR) of CD8+ glycoprotein cells (cytotoxic or killer T-cells) [27]. A digested antigen peptide consisting of 8-13 amino acid residues in length (most commonly 9 to 10 residues) is specifically bound within the peptide-binding groove located on the HLA surface [28,29]. The binding groove is constructed from the eight ß-strands and α1 and α2 helixes, while α3 and ß 2 -microglobulin (ß 2 m) are the floor and non-covalent supported units, respectively, as shown in Fig 1 [28]. The different diversity of residues present in the binding-cleft markedly affects the specific binding to the antigen peptide and pathogenesis. The residues at or close to the two ends of the peptide, P2/ P3 and P9 residues for the MICA-TM nonapeptide, generally interact with the binding groove of HLA class I and serve as an anchor, whilst the central peptide oozes up out of the groove for direct interaction with the TCR [30,31]. However, peptides longer than nine residues and up to 13 residues can be recognized by the HLA class I genes through "the middle-out loop shape form" [30], where the residues at or close to the peptide terminals are well occupied in their sub-sites [29].
The association of BD with the HLA-B Ã 51:01 allele in various ethnic groups in Japan and Korea has been frequently observed [3], while the HLA-A Ã 26:01 allele was detected in BD patients in Greece, Japan and Taiwan [32][33][34]. Some possible candidate genes for BD have been studied in various HLA-A and-B alleles [35]. Recently, in Japanese uveitis patients, BD was associated with the HLA-A Ã 26:01 allele at a 37.5% phenotype frequency more than the controls, and so HLA-A Ã 26:01 is a possible marker as a susceptible allele for ophthalmic BD in Japanese ethnics [36]. By meta-analysis and positive pathergy tests of the Japanese data, the BD clinical manifestations of uveitis, skin lesions and arthritis, and genital ulcers were found to be significantly associated with the HLA-A Ã 26:01, A Ã 02:07 and A Ã 30:04 alleles, respectively [37]. Currently there is no treatment (e.g. vaccine) against BD based upon the HLA recognition antigens despite the apparent restricted association with specific alleles, since the individual risk factors with HLA are unknown. Therefore, we aimed to investigate the link between BD and two specific HLA alleles associated with BD (HLA-A Ã 26:01 and HLA-B Ã 51:01) in terms of their binding affinity to the MICA-TM peptide using MD simulations, in comparison with two HLA alleles not associated with BD (HLA-A Ã 11:01 and HLA-B Ã 35:01). The molecular understanding of the specific binding of the MICA-TM peptide with the HLA alleles associated with the disease may be useful for therapeutic vaccine design.

Structure preparation
The starting structures of three of the HLA class I alleles were taken from the Protein Data Bank, being entries 1E27 [38] for B Ã 51:01 with the HIV-1 epitope at 2.20 Å, 1X7Q [39] for A Ã 11:01 with the SARS nucleocapsid at 1.45 Å, and 1A9E [40] for B Ã 35:01 with the EVB peptide at 2.5 Å (S1 Table). That for A Ã 26:01 was built from homology modeling using the above A Ã 11:01/SARS nucleocapsid structure as a template and the amino acid sequence of A Ã 26:01 from GenBank (accession no. AAA03720) [41], which has 94.3 and 96.6% amino acid identity and similarity to A Ã 11:01, respectively. The HLA/MICA-TM complexes were constructed by changing the original peptide within the X-ray structure to be consistent with the nine amino acids of the antigenic MICA-TM peptide (AAAAAIFVI) using the align sequence profiles module implemented in the Discovery studio 2.5 (Accelrys, Inc.).

MD simulation
The four studied HLA/MICA-TM complexes were investigated by molecular dynamics (MD) simulations using the AMBER 10 software package [42] with the ff03 force field. The missing atoms were added using the LEaP module [43]. The protonation state of all possible charged residues (arginine, lysine, histidine, aspartate and glutamate) in HLA-allele complexes was assigned at pH 7.0 by PROPKA server [44]. The total charges with negative value of the HLA/ MICA-TM complexes were randomly neutralized by Na + counterions (2, 1, 1 and 2 ions for the B Ã 51:01, B Ã 35:01, A Ã 26:01 and A Ã 11:01 systems, respectively). Afterwards, the individual complex was solvated by TIP3P [45] water molecules leading to approximately 65,000 atoms in total. The dimensions of the simulation box used for all the systems were 86 × 90 × 88 Å 3 . The periodic boundary condition with the NPT ensemble and a simulation time step of 2 fs was used. All energy minimizations and MD simulations were performed using the SANDER module of AMBER 10. All bonds and angles concerned to hydrogen atoms were constrained by algorithm of SHAKE [46]. The long-range electrostatic interactions were treated by particle mesh Ewald method and the non-bonded interactions with a cutoff distance of 12 Å were considered [47]. All MD simulations were run with a 12 Å residue-based cutoff for non-bonded interactions and the particle mesh Ewald method was applied for an adequate treatment of long-range electrostatic interactions [48]. Each system was subjected to the four stages of the restrained MD simulations at 298 K with force constants of 10, 7.5, 5 and 2.5 kcalÁmol -1 ÁÅ 2 for 500 ps in each stage accordingly. These subsequent steps could allow the peptide to adapt its geometry and orientation from the initial model to fit better within the peptide binding groove. Then, the constraints were completely removed and fully unrestrained MD simulations were performed until 50 ns. The convergences of energies, temperature, and global root meansquare displacement (RMSD) were used to verify the stability of the systems. The MD trajectories were collected every 0.2 ps from the production phase for further analysis.

Stability of the HLA/MICA-TM complexes
To study the structural stability of the four MD simulations, the RMSDs for all atoms of the four HLA alleles (B Ã 51:01, B Ã 35:01, A Ã 26:01 and A Ã 11:01 alleles) complexed with the MICA-TM peptide compared with those of the starting structures were monitored along 50 ns of simulation time using the PTRAJ module implemented in AMBER10 package [49]. The RMSD fluctuations of each complex and its structural components (binding groove and ß 2microglobulin) and MICA-TM peptide, were plotted in S1 Fig. The RMSD values quickly increased until~5 ns and then reached a plateau except for the HLA-B Ã 51:01/MICA-TM complex where the RMSD value continuously increased over the first~15 ns. The MICA-TM in HLA-B Ã 51:01 system had a low fluctuation around 1 Å, and suddenly at 23 ns, it jumped up to 2 Å until the end simulation. The RMSD increasing caused from the non-polar of MICA-TM peptide change orientation to hydrophobic region within pocket. Although the whole HLA/ MICA-TM complex for the two BD-associated HLAs (B Ã 51:01 and A Ã 26:01) fluctuated a great deal, their binding groove and the incoming MICA-TM peptide demonstrated a rather low level of fluctuation. All systems seemed to reach equilibrium after 25 ns with~1 Å of RMSD fluctuation, and so the MD trajectories from 25 to 50 ns (the production phase) were used for analysis.

Regional flexibility
The backbone flexibility of the HLA-B Ã 51:01, B Ã 35:01, A Ã 26:01 and A Ã 11:01 complexes with MICA-TM were investigated by B-factor calculation over the last 25 ns trajectories (Fig 2). Note that high flexibility of protein was displayed in red and vice versa in blue. Among the four systems, the HLA-A Ã 11:01 allele showed the highest degree of protein flexibility in the region of peptide binding groove and in particular at the S7-S9 sub-sites (defined in Fig 1D), whereas the other HLA alleles were likely to have similar degree of flexibility. This might be due to the

Per-residue decomposition (DC) energy
The decomposition (DC) energies were calculated and used to scan for potentially important residues for binding [50]. In order to seek the fingerprint of the HLA/MICA-TM interactions, the interaction energy between each HLA residue and the MICA-TM peptide and vice versa was calculated over 100 snapshots of the production phase. The obtained DC energies of the HLA and MICA-peptide residues are plotted in Figs 3 and 4, respectively. The HLA/MICA-TM interactions mostly occurred on the α1 and α2 helixes with additionally some residues of β-strands (Fig 3). Using the criterion of a total DC energy < -0.5 kcal/mol as an important residue, then 22, 18, 16 and 14 potentially important residues of the HLA-B Ã 51:01, B Ã 35:01, A Ã 26:01 and A Ã 11:01 alleles (S2 Table), respectively.
Van der Waal (vdW) interactions were found to play an important role in the complex, where the magnitude of the MICA-TM peptide binding with the BD-associated HLA alleles (B Ã 51:01 and A Ã 26:01) was greater than that for the corresponding matched non-BD-associated HLA alleles (B Ã 35:01 and A Ã 11:01). For the nonapeptide, as expected the non-polar residues interacted with HLA through vdW interactions (Fig 4). Based on the definition of a total DC energy of < -3 kcal/mol for strong binding, there were four peptide residues that firmly bound to HLA-B Ã 51:01 (P3, P5, P8 and P9) and A Ã 26:01 (P2 and P7-P9), while only two and three residues were found in B Ã 35:01 (P7 and P8) and A Ã 11:01 (P2, P7 and P8), respectively. Note that the P9 and P2 or its adjacent residue are known as the anchor for the peptide accommodation at the binding groove of HLA class I [30,31]. Loss of the P9 contribution may lead to an unbinding recognition of the non-BD-associated HLAs towards the MICA-TM peptide. More details of the interaction and orientation of this peptide at the binding groove of the four HLA alleles investigated are discussed in the hydrogen bond pattern section.

H-bond patterns between HLA and the MICA-TM nonapeptide
To investigate the intermolecular H-bond interactions between the HLA protein and the incoming short MICA-TM nonapeptide, the number of H-bonds was evaluated using the two acceptance criteria of (i) a distance between the proton donor (D) and acceptor (A) atoms of D−A 3.5 Å; and (ii) an angle of D−H . . . A > 120°, as previously reported [50,51]. The obtained results were compared between the two paired BD-associated and non-associated HLA alleles in each locus (A or B), and are shown in Fig 5. The peptide binding sub-sites (S1-S9) were classified by the peptide contact position in the binding groove of HLA, as seen in the HLA-B Ã 51:01/peptide complex (Fig 1D). The last snapshot of each complex was used to represent the protein-protein H-bond formation detected from the MD simulations (see in S2 Fig).
As shown in Fig 1D, the peptide binding groove of HLA class I is constructed by at least five sub-pockets with a large groove volume (S1 and S2 yellow; S3 magenta; S4-S6 dark blue; S7 orange; S8 and S9 light blue), and it is able to support a variety of incoming short peptides in the binding step. Moreover, the binding groove of the HLA protein provides a hydrophobic cavity to support the nine spanning residues (AAAAAIFVI) of the MICA-TM peptide. Due to the nonpolar nature of this peptide, H-bond formation with the HLA protein was expected through the peptide backbones. In Fig 5A and 5B, both HLA-B alleles similarly stabilized the P9 residue of the MICA-TM peptide by the formation of three H-bonds which two of them are the salt bridge interactions between the C-terminal carboxylate group and the ammonium group of K146. The C-terminal residue (P9) could bind significantly stronger with the BDassociated HLA allele than the non-associated allele (~60-90% H-bond occupation compared to only~50-70% in the non-associated allele), which could be due to the different residue 80 on the α1 helix (I80 in HLA-B Ã 51:01 and N80 in HLA-B Ã 35:01). Through the strong H-bond (70% occupancy) with the amide group of N80 (Fig 5B), this polar residue at the S9 sub-site of the peptide induced the C-terminal to change its orientation and move closer to the α1 helix of the non-associated BD HLA-B  peptide had no interaction with any HLA-B Ã 51:01 residues but instead the P3 residue was stabilized by two α1-helix residues, Y9 and N70, with *90 and *50% H-bond occupations, respectively (Fig 5A and S2 Fig). The rearrangement of the H-bonding network at the N-terminal was previously observed in the octamer and nonamer peptides binding to HLA-B Ã 51:01 [38]. In addition to the protein-protein interactions at the two peptide ends, a strong H-bond was formed between the backbone of the P5 Ala residue and the side chain of the T73 residue on α1 helix in the P4-P6 pocket in both HLA-B alleles (> 80%, Fig 5A and 5B). This is congruent with the previous observation that the P5 residue of HIV epitopes (Ile, Val and Pro) was found as an unexpected anchor deeply pointing into the sub-pocket of the binding groove of HLA-B Ã 51:01 [38]. In addition, these two HLA-B alleles had a strong H-bond between the P8 backbone of the MICA-TM nonapeptide and the W147 indole ring on the α2 helix.
In the case of the two HLA-A alleles (Fig 5C and 5D), the N-and C-terminal regions of the MICA-TM peptide showed a firmly established interaction, P1 with the Y159 phenyl ring (~100%) and P9 formed two H-bonds or salt bridge interactions with the K146 ammonium group (~70-80% in the BD-associated HLA-A Ã 26:01, but < 60% in the non-associated one), respectively, at the α2 helix. Only HLA-A Ã 26:01 ( Fig 5C) contained a strong interaction at the nearly C-terminal end (P8 residue) with the α2 helix residue W178 (~70%) and the anchor P5 residue weakly interacted with the α1 helix residue T73 (~40%). The Y99 residue on the βsheet of HLA-A Ã 26:01 somewhat stabilized the P3 residue of MICA-TM (~40%), while the α1 helix residue N66 in HLA-A Ã 11:01 supported the P2 and P3 residues (~40-50%). It is worthy to note that the high protein flexibility of HLA-A Ã 11:01 (Fig 2D) had consequently led to a low binding of the incoming peptide, and in particular at the P5-P9 residues.
Based on the formed H-bonds, the peptide binding recognition was better distinguished in the HLA-A alleles. The lowered binding strength at the C-terminal P9 residue observed in HLA-B Ã 35:01 was suspected to be the most important reason for the low selective binding affinity of this non-BD-associated HLA-B Ã 35:01 allele. Thus, the relatively high protein flexibility in the non-associated HLA-A Ã 11:01 allele led to a decreased H-bond strength with the P9 residue and no stabilization for the P5 and P8 residues. All simulations fairly agreed with the crystal structures of the HLA/peptide complexes in which this P6 residue outwardly located off the binding groove [31].

HLA/peptide binding affinity
To determine the MICA-TM peptide binding strength towards the four studied HLA alleles, the molecular mechanic/Poisson-Boltzmann surface area (MM/PBSA) approach was applied on the same set of 100 snapshots taken from the production phase. The MM/PBSA approach has been extensively used to predict the overall Gibbs free energy (ΔG bind ) in biomolecular systems [52]. The ΔG bind of the HLA alleles with the peptide bound at the binding groove was calculated from a summation of the total MM energy (ΔE MM ), the solvation free energy (ΔG sol ) and the entropic term (TΔS). The last term was estimated from the normal mode analysis [53]. Although the ΔE MM components suggested that the short peptide was better stabilized by the non-BD-associated HLAs in the gas phase, the large destabilization from the solvation effect on the HLA/peptide complex led to a lower total ΔG bind . Rather, within each paired HLA-I gene (HLA-A or HLA-B), the MICA-TM peptide interaction was stronger with the BD-associated allele, where B Ã 51:01 (-56.10 kcal/mol) > B Ã 35:01 (-48.88 kcal/mol) and A Ã 26:01 (-42.97 kcal/mol) > A Ã 11:01 (-31.14 kcal/mol) in Table 1. The obtained results from the present study help to differentiate the HLA alleles and explain a source of BD.

Conserved interaction of the BD associated HLAs
Based on MM/PBSA method, the alanine scanning mutagenesis commonly used to signify the important function of residues was carried out on the two BD associated HLA residues within the 5 Å sphere of the MICA-TM peptide. The results of the relative binding free energy (ΔΔG binding = ΔG wild-type − ΔG mutant ) are given in Table 2. Note that the entropy term of the complex is not significantly changed by only one residue substitution with alanine; therefore the entropic calculation was neglected. By a mutation to alanine on the 40 residues, the HLA residue numbers 70, 73, 99, 146, 147 and 159 contributing the ΔΔG value < -2 kcal/mol in either HLA-B Ã 51:01 or HLA-A Ã 26:01 (Table 2) were considered as the important residues. These six residues could provide the conserved interaction of both BD associated HLAs for binding of the incoming MICA-TM.

Conclusions
BD is a multi-organ inflammatory disorder with vasculitis, in which the main cause and mechanism of action are still not well understood. From clinical investigations, the HLA-B Ã 51:01 and HLA-A Ã 26:01 with MICA Ã 009 containing the nonamer peptide (AAAAAIFVI or MICA-TM) are most frequently detected in BD patients. Herein, we aimed to search for the selective correlation between the two BD-associated HLA alleles (HLA-A Ã 26:01 and HLA-B Ã 51:01) and the MICA-TM peptide in comparison with two class matched non-BD-associated HLA alleles (HLA-A Ã 11:01 and HLA-B Ã 35:01) by MD simulations. From the simulations, more contact residues at the binding groove of the BD-associated HLA alleles stabilized the incoming MICA-TM peptide than at the non-associated alleles of the same class (22 and 16 residues for B Ã 51:01 and A Ã 26:01; 18 and 14 residues for B Ã 35:01 and A Ã 11:01). The vdW force was found to be the main protein-protein interaction, but in addition strong H-bonds (>70% occupation) were likely formed with the backbone of the nonpolar peptide, and these were stronger in the BD-associated HLA alleles. The P2/P3 and P9 residues (close to and at the peptide ends, respectively) acted as the anchor for the peptide accommodation at the binding groove of the BD-associated HLAs. The total binding free energy of the HLA/peptide complex suggested a significantly stronger binding strength of the MICA-TM peptide with the Data are shown as the mean ± SD, derived from independent simulations. Means within a paired row (HLA-A or HLA-B alleles that are associated with BD versus that are not) followed by a different letter are significantly different. ΔG bind is the binding energy with inclusion of entropic term.   Table. HLA residues with a total decomposition (DC) free energy (ΔG residue ) of less than -0.5 kcal/mol.