Abstract

MHC α-helices form the antigen-binding cleft and are of particular interest for immunological reactions. To monitor these helices in molecular dynamics simulations, we applied a parsimonious fragment-fitting method to trace the axes of the α-helices. Each resulting axis was fitted by polynomials in a least-squares sense and the curvature integral was computed. To find the appropriate polynomial degree, the method was tested on two artificially modelled helices, one performing a bending movement and another a hinge movement. We found that second-order polynomials retrieve predefined parameters of helical motion with minimal relative error. From MD simulations we selected those parts of α-helices that were stable and also close to the TCR/MHC interface. We monitored the curvature integral, generated a ruled surface between the two MHC α-helices, and computed interhelical area and surface torsion, as they changed over time. We found that MHC α-helices undergo rapid but small changes in conformation. The curvature integral of helices proved to be a sensitive measure, which was closely related to changes in shape over time as confirmed by RMSD analysis. We speculate that small changes in the conformation of individual MHC α-helices are part of the intrinsic dynamics induced by engagement with the TCR.

1. Introduction

T cells play a major role in both innate immunity and adaptive immunity. Their surface-bound T cell receptors (TCRs) recognise antigens and thereby detect hazardous organisms or changes inside cells. TCRs recognize short peptide fragments (p) that are bound to major histocompatibility complexes (MHC). Different interactions of TCR/peptide-MHC (TCR/pMHC) are believed to be the basis for distinctive stimuli that lead to and trigger different fates of the T cell, for example, T cell development, thymic selection, lineage commitment, differentiation into effector cells, or memory T cell responses to foreign antigens [1].

MHCs are surface-bound proteins and their role is to present peptide fragments to TCRs, be they self- or alloantigens. To achieve this, MHCs have a cleft that is able to bind peptide fragments. This cleft comprises two α-helices and five subjacent lateral β-strands forming a sheet or floor of the cleft. TCRs engage peptide MHCs in a diagonal arrangement that seems to be a common mode of interaction across TCRs [2, 3]. The two MHC α-helices interact with the TCR complementarity determining regions (CDR) 1 and 2, while the MHC-bound peptide interacts with CDR3, although CDR1 has also been shown to interact with the terminal parts of the peptide. Most of the sequence variability of TCRs is found within CDRs; these regions are also referred to as hypervariable regions. The two MHC α-helices are of particular interest as they represent stable secondary structural domains interacting with TCRs.

Adhesion and signaling proteins together with a set of TCR/pMHC complexes form a junction between T cell and an antigen-presenting cell that is referred to as an immunological synapse. It serves as a platform for assembly and segregation of signaling complexes, which cooperatively decide the outcome of T cell activation and effector function. New findings show that this supramolecular complex forms too late to be relevant for initial TCR signaling that happens within seconds after pMHC engagement, with the TCR initiating a tyrosine phosphorylation cascade [4]. Such an early signal may be sufficient to trigger effector function like killer T cell cytolysis of target cells. Stimulation of proliferation, however, requires engagement and signaling for many minutes or even hours [58]. It is important to note that αβ-TCR itself has no signaling motif but associates with homo- and heterodimeric cluster of differentiation 3 (CD3) subunits, , , and in a noncovalent way. These subunits contain immunoreceptor tyrosine-based activation motifs (ITAMs) that can be phosphorylated and initiate downstream signaling of TCR activation.

Protein structures as found in protein databases (e.g., Protein Data Bank (PDB)) show a static image. In contrast to that, in molecular dynamics (MD) simulations, the protein explores many conformations and allows one to capture its dynamics. Computational immunoinformatics is a well-established, rapidly evolving field [9]. In previous papers [10, 11] we presented a first evaluation of three TCR/pMHC systems that differ only slightly in the MHC amino acid sequence. Macdonald et al. [12] determined binding characteristics and immunogenicity of MHC alleles 44:02, 44:03, and 44:05 in complex with the ABCD3 self-peptide (EEYLQAFTY) and LC13 TCR. 44:02 and 44:05 trigger an immune response when bound to the LC13 TCR, while 44:02 does not despite extensive amino acid sequence homology. This renders these HLA alleles an interesting set to study TCR allorecognition. Macdonald et al. [12] determined the X-ray structure of the ternary complex 44:05, ABCD3 nonapeptide, and LC13 TCR that is accessible on the protein database (PDB, http://www.pdb.org/, PDB-ID: 3KPS). Due to extensive sequence identity we were able to use in silico site-directed mutagenesis to obtain 3D structures of the missing TCR/pMHC complexes. MD simulations in the nanosecond range could possibly show short-lived changes in dynamic behaviour, conformation, propagation of forces, or early activation signals [13, 14].

The aim of the present paper is to put forward adequate tools for identifying and monitoring of conformational changes with possible functional relevance in MHC α-helices, and in particular to monitor geometric characteristics of the MHC peptide-binding groove. Here, we present (i) a robust and parsimonious method to find an approximation of a protein’s α-helical axis, (ii) an evaluation of polynomial degree adequacy in describing bending or hinge movements of particular α-helical axes, and (iii) application of polynomial fitting of α-helical axis to monitor α-helical conformations along MD simulations in different TCR/pMHC immunological complexes.

Our previous studies established the use of polynomials to model α-helices in MHC molecules and monitor their dynamic behaviour [15]. To mathematically describe the structural dynamics of MHC α-helices at the TCR/pMHC binding interface we first identified those helical regions which were stable and therein those which are close to the protein-protein interface. Then we extracted the α-helical axis and finally determined a polynomial that approximates this axis in a least-squares sense.

Various methods to extract a helix axis have been developed [16], including rotational fitting, using atoms as control points of B-splines [17] or fitting to a helix. We used a fragment-fitting method, based on previous work [16], to locate the axis of α-helices as follows: an ideal, linear α-helix fragment comprising four atoms is superimposed on successive pieces of MHC α-helices in a least-squares sense. Along the axis of the fitted helical fragment, we adopt points as estimates of the MHC α-helix axis and fit a polynomial through these points. From the polynomial, geometric parameters can be derived to monitor conformational changes. Polynomials fitted to the α-helical axis can, in principle, be of any degree. However, polynomials of higher order tend to oscillate, adding noise to geometrical quantities computed thereof. We therefore evaluated polynomials of different degrees for their ability to reproduce bending and hinge motions of an α-helix with minimum relative error. Between two adjacent α-helices, as found in MHC proteins, the polynomials serve to span a ruled surface. This interhelical surface lends itself to derive several quantitative characteristics of shape: (a) total area [18, 19], (b) a profile of interhelical distances along the binding cleft, and (c) heuristic “centre line of the cleft” which may be constructed, along which the surface torsion, that is, a twist or screw of the interhelical surface, can be computed. The latter characterizes the positions and bending of helices relative to each other and defines the geometrical shape of the peptide-binding cleft that is ligated to the TCR. Dynamics in the shape of the protein-protein interface might modulate the TCR/pMHC binding affinity.

2. Methods

2.1. Homology Modelling of TCR/pMHC Complexes

Three TCR/pMHC systems listed in Table 1 were simulated. The X-ray structure of TCR/pMHC B4405 (number 3 in Table 1) was taken from the PDB (PDB-ID: 3KPS). Structures B4402 and B4403 were engineered by means of in silico site-directed mutagenesis [20] using B4405 as a structural template. We introduced(i)mutation Y116D to the MHC molecule to get LC13/ABCD3/44:02 complex (B4402),(ii)mutations Y116D and D156L to the MHC molecule to get LC13/ABCD3/44:03 complex (B4403).See Figure 1(b) for a 3D representation of amino acid positions Y116 and D156. 3D structures were edited and mutations introduced with the Swiss PDB Viewer. The program automatically browses a rotamer library and selects an amino acid rotamer minimizing a scoring function in order to fit the new amino acid in its surrounding and avoid steric clashes with other atoms. All MHC molecules simulated have an amino acid sequence identity of more than 99% and stay in very similar 3D fold. MHC molecules appeared stable during all our MD simulations as seen in RMSD plots in Figure 15.

The full amino acid sequence of HLA-B alleles is accessible from the HLA library (https://www.ebi.ac.uk/ipd/imgt/hla/) and a description of its topology is accessible from UniProt (http://www.uniprot.org/uniprot/Q95365). Not surprisingly, a sequence comparison showed that a transmembrane helix (24 amino acids) and cytoplasmic tail (30 amino acids) are missing from the MHC X-ray structure as plasma membrane structures and flexible protein parts are hard to determine using X-ray crystallography. The LC13 TCR is also missing its transmembrane helix.

2.2. Molecular Dynamics Simulations

GROMACS 4.0.7 was used for molecular dynamics simulation. First, water molecules were added to the protein structure, immersing it in an artificial water bath of rectangular form and allowing a minimum distance of 2 nm between the protein and the box boundaries. Second, water molecules were replaced by sodium and chloride ions to yield a salt concentration of 0.15 mol/L and neutralize the protein net charge. Third, the energy of the solvated system was minimized using a steepest descent method and then the system was gradually heated up to 310 K during 100 ps position restraint MD simulation. Finally, MD production runs were done with the LINCS constraint algorithm acting on all bonds using the Gromos 53A6 force field [21]. Hydrogen and fast improper dihedral motions were removed, allowing for an integration step of 5 fs. Van der Waals and Coulomb interactions were computed using a cut-off of 1.4 nm. Long-range Coulomb interactions were computed by Ewald summation. Velocity rescale temperature coupling was set to 310 K and Berendsen isotropic pressure coupling was set to 1 bar. The selection of the force field and MD parameters for pMHCs were evaluated by Omasits et al. [22] and set accordingly.

2.3. Finding Dynamically Stable -Helices at the Protein-Protein Interface

As outlined in the introduction, the MHC protein comprises α-helices and beta-sheets. We are interested in monitoring atoms that form stable α-helices and are in close contact with the TCR. From MD simulation data we calculated the following:(i)The relative presence of α-helical structure in the protein complexes over the simulation time: we used the DSSP algorithm [23] as implemented in GROMACS to identify secondary structural elements of the protein complex over simulation time. α-helices were considered stable if the ratiowith being the time the main chain (backbone atoms plus carbonyl oxygen atoms) meets the DSSP criterion for an -helix and being the total simulation time. The cut-off value of 0.5 is justified by the distinctly bimodal distribution (see Figure 11(a)).(ii)The relative presence of close contacts at the protein-protein interface: atoms that are less than 1.4 nm apart (i.e., the cut-off for electrostatic interactions in our MD simulations) are defined as being in close contact. Atom-atom contacts between TCR and MHC are defined stable if the ratio is the time during which two atoms are in close contact, and is the simulation time. Again, the cut-off value of 0.5 is justified by the distinctly bimodal distribution (see Figure 11(b)). To get the residue-wise relative contact time, we averaged the atom-wise relative contact time (defined in (2)) over all atoms per residue.The resulting sets of amino acid residues defined by procedures (i) and (ii) were intersected (workflow, see Figure 2; results, see Figure 3) before applying further methods described in Sections 2.4 and 2.5. Dynamically stable helices (green atomic surface in Figure 3) as defined by procedure (i) overlap well with the atoms in close contact of the TCR (red atomic surface in Figure 3) especially with the cut-off set to 1.4 nm. Note that at the end of the MHC’s peptide-binding pocket both G-ALPHA1 and G-ALPHA2 exhibit a kink followed by a short α-helix. These short α-helices are present in the crystal structure, but we found them being unstable during MD simulations as they fold and unfold. We do not consider these helices in our analysis, as they are not in close contact with the TCR.

2.4. Approximating the Axis of an -Helix

In order to mathematically describe and quantify α-helical geometry and movements, polynomials of degree were fitted to the α-helices, where represents a vector of 3D coordinates and is the curve parameter. Prior to fitting we extracted atom coordinates of those amino acids, which fulfil the criterion of stable α-helices and close contacts as described in Section 2.3.

According to the structural definition of -helices [24] a model of one ideal α-helical turn, that is, a fragment consisting of four atoms, is constructed, with its axis coinciding with -axis. The coordinates of its th were assigned aswith pitch (advance from one to the next), helix radius , and . Along the axis of this helical fragment we consider three axis points, the initial , the intermediate , and the final .

Out of an α-helix with    atoms we pick moving groups of four successive atoms each, to which we fit the fragment model defined above in a least-squares sense [25]. Along -axis of the fitted fragment model we adopt points as estimates of the axis of the α-helix, see Figure 4. From the very first fitted fragment ( atoms ) we adopt two points as points of the helix: the transformed initial axis point as and the transformed intermediate point as . From fragments fitted subsequently we adopt only the respective intermediate points, and from the last fragment (fitted to atoms ), we again adopt 2 points, its “intermediate” as and its final point as . Thus, points represent the axis of the α-helix.

2.5. Fitting the Axis of an -Helix

The points representing the helical axis are fitted by polynomials of degree in a least-squares sense. Separate regression functions , , and are computed for -, -, and -coordinates:for example, for -coordinate,The total number of parameters for this model is . When fitting the curve, parameter is evaluated at discrete values of pitch, corresponding to the positions of estimated points along the helical axis. After the regression has been performed, in (5) may be evaluated for arbitrary, continuous values , yielding a continuous representation of the helical axis.

Both G-ALPHA1 and G-ALPHA2 helices were modelled in the same way, yielding models and with equal polynomial degrees . It is well known that (with equidistant data points) interpolating polynomials of too high a degree may exhibit severe oscillations near the ends of the interpolation interval [26]. (Interpolating polynomials actually pass through all data points.) This is also true for approximating polynomials [27]. (Approximating polynomials need not pass through data points but rather approximate the shape of their functional dependence in a least-square sense.) We therefore kept the polynomial degrees as low as possible.

2.6. Geometric Quantities
2.6.1. Interhelical Distance and Area of Interhelical Surface

For each polynomial, the curve parameter ranges within , with possibly different values for each helix. We consider equidistant values of , yielding reference points on each helix model. Connecting corresponding reference points by straight lines yields rulings, , which span a ruled surface (see Figure 5). (The rule for defining “corresponding” points has to be adopted in a reasonable way but finally remains to some degree arbitrary.) From rulings we estimate distances between polynomials across the cleft. The resulting polygon mesh was triangulated and interhelical area was calculated, as previously outlined [18, 28]. Each of these quantities may be monitored over time, for example, ; see Figure 9. Respective graphs provide well-defined estimates of changes in width of the intrahelical gap (binding cleft) as a function of both helical position and time. Likewise, median, quartiles, and extreme values of interhelical distances for each and of can be obtained.

2.6.2. Torsion of Interhelical Surface

The ruled surface between both polynomial helix models may bend and wind in various ways. Describing all aspects would call for a comprehensive mathematical treatment in terms of differential geometry, from which we refrain. Instead, we restrict ourselves to describe something like the “torsion of the interhelical surface” in a simple, intuitive way; see Figure 5. (In everyday terminology “torsion” is often called “winding.” It may apply to (curved) lines as well as to surfaces in 3D space. We use both terms in parallel with equal meaning.)

To this end we note that the observer may slide across a surface along various paths and may, in general, observe different values of surface torsion along each of the paths. We notice that torsion in strict mathematical terms is a local characteristic of the surface, and even more, it depends upon the path one takes to inspect the surface. For the sake of simplicity we deliberately adopt the centre line between both helix models:and let be the unit tangent vector of . Of note, the centre line intersects with each of the rulings. While proceeding along the centre line we inspect the directional change of rulings. This provides us with a particular characterization of the deformation of the interhelical surface. To quantify this intuitive description, we introduce a parameter “surface torsion.” It is based on the change in direction between successive rulings [29] and is determined by the derivative . (In strict mathematical terms (of differential geometry) this quantity is called “parameter of distribution,” which we think is a very nonintuitive label, prone to be mixed up with probability distributions. Hence, we prefer the more intuitive term “surface torsion” to quantify the winding of a surface while proceeding along a given path (in our case the centre line), not that surface torsion is generally different in different directions, even in the very same point of a surface.) What we call surface torsion is given by [29]The integral valuequantifies the relative tilt of the surface between start and endpoint of the path. Note that indicates right-handed surface torsion, indicates left-handed surface torsion, and indicates that the surface at this point is developable and rulings are parallel. (“Developable” means that a sheet of paper could be bended to exactly match the surface; the surface could be “flattened.”) In summary, the definition of the intrahelical surface depends on a reasonable selection of rulings, , and a path, , which are both arbitrary to some extent. Despite these shortcomings in a strict mathematical sense, the concept of torsion, as introduced here, mirrors some essential features in describing the interplay between the shape of -helices and their relative orientation in forming the MHC binding cleft.

2.6.3. Curvature of Helices

Derivatives of each polynomial helix model can be obtained analytically for each value of the parameter , yielding the vectorsCurvature is a scalar quantity being per definition positive in Euclidian 3D space and is obtained via [29]:

2.6.4. Construction of a Curved Helix Model and Helix Hinge Model

So far we have described how an α-helical axis is reconstructed by our newly proposed fragment-fitting method and fitted by a polynomial of a certain degree. From this polynomial, we calculate the curvature integral as described in Section 2.6.3. We wanted to find an optimal degree for these polynomials in order to retrieve conformational deformations with minimal errors. To do so we modelled two different motions such that MHC α-helices could perform: a bending motion and a hinge motion (see Figure 6).

To model a bending motion we created a curved helix backbone model with known curvature, as described in detail in Appendix A. Following polynomial representation of the backbone, the curvature integral can be obtained analytically, yielding the reference, (see (A.2) in the appendix).

To test our method, the curved helix model was subjected to the fragment-fitting method and the resulting helix axis was fitted by a polynomial in a least-squares sense as described in Sections 2.4 and 2.5. The polynomial was evaluated at 100 equidistant points and curvature and the curvature integral were calculated numerically, yielding (see (10)). As a quality criterion for regaining the correct curvature integral we used the relative error and evaluated it for polynomial degrees 1 to 8. The results are depicted in Figure 6(b).

To model a hinge motion, we constructed an ideal, linear -helix comprising 31 atoms. Subsequently, the helix was split into two parts: one atom was selected and the remaining part of the helix rotated around the selected atom to model a hinge motion (see Figure 6(c)). The position to split the helix (number of selected) was varied from atom 5 to 15. The aim was to examine α-helix hinges with two legs unequal in length. From this series of α-helical models (different positions of the hinge and varying hinge angles, ) we then reconstructed the α-helical axis using the fragment-fitting method. A polynomial of th order was fitted to the traced α-helical axis in a least-squares sense and the curvature integral calculated numerically, yielding . The relative error was obtained for polynomials of degrees 1 to 8, with being the angle between the two parts of the helix.

3. Results

3.1. Finding the Optimal Polynomial Degree for Monitoring Helix Motions

To monitor deformations of α-helices in MD simulations we propose an approximation of the α-helical axis using a fragment-fitting method (Section 2.4), fitting the resulting axis by a polynomial in least-squares sense (Section 2.5), and derive geometric quantities from it (Section 2.6). To find an optimal degree of polynomials we applied our method by applying it to modelled α-helical bending and hinge motions (Figures 6(a) and 6(b)) and evaluated the polynomial degree that retrieves a certain α-helical motion with minimal relative error. For the hinge motion we noticed that relative errors increase with polynomials of higher order. For the bending motion, choosing sixth- or eighth-degree polynomials increases the degree of freedom while at the same time failing to add adequate improvement in relative error. Generally, relative errors in the detection of hinge motions are larger than those for bending motions. The second-order polynomials reproduced the measured quantities of bending and hinge motions with low relative error and hence were adopted for monitoring MHC α-helices in all subsequent analyses. We admit that second-order polynomials may fail to model α-helical axes in full detail, but it seems appropriate for capturing trends in helical motions and shape during an MD simulation with minimal relative error.

We applied the geometric analysis to model MHC α-helices of TCR/pMHC complexes B4402, B4403, and B4405. By inspecting the curvature integral for MHC α-helices, G-ALPHA1 and G-ALPHA2, three interesting cases could be identified showing changes along the 250 ns MD simulation (Figure 7, panels in upper row). Phases during which the curvature integral changes either abruptly ((b) and (c)) or gradually (a) are highlighted in orange and red. Helix conformations (“bundles”) corresponding to these phases are shown in panels in the middle row. To check if changes in the curvature integral actually indicate a conformational change we calculated RMSD matrices of conformations within and between orange and red helix bundles (panels in lower row). In all three cases, RMSD between phases is distinctly larger (shifted to the right) as compared to RMSD within phases. A shift towards higher RMSD values indicates a conformational change and is seen in all three cases (see Table 2 for median RMSD values). These case studies demonstrate that changes in conformation of α-helices during MD simulations can be monitored using second-degree polynomials fitted to helical axes.

3.2. Geometric Quantities Characterizing the Shape of MHC -Helices

The MHC peptide-binding groove comprises two α-helices that interact with the TCR. Using second-order polynomials, we computed (i) the integral of the curvature of individual MHC α-helices, (ii) the area of interhelical surface, and (iii) the surface torsion along the imaginary centre line derived from both polynomials that model MHC α-helices for single time steps in MD simulations. Items (ii) and (iii) are geometric properties of the ruled surface. They are used to quantify the geometric relation between the two helices, for example, their relative orientation, and thus capture important aspects of the geometry of the peptide-binding groove, as illustrated in Figure 5.

Curvature is a local feature of a curve. Considering the curvature integral we obtain a measure of the overall bending of the whole curve [30]. Curvatures of polynomials modelling single helices were integrated and monitored over time as seen in Figure 8. G-ALPHA1 of B4403 and B4405 each undergoes abrupt fluctuations in helical conformation, reflected in the curvature integral, but is stable in the phases in-between. The curvature integral for G-ALPHA1 of B4402 shows a gradually increasing trend. We notice that the curvature integral for G-ALPHA2 is generally higher than that for G-ALPHA1, reflecting the kink near its N-terminal end. The two helical parts of G-ALPHA2 form a hinge. G-ALPHA2 of B4403 and B4405 shows only minor fluctuations and no abrupt changes in the curvature integral indicating that the hinge angle stays stable. For B4402, the hinge angle decreases in the first half of the MD simulation and remains stable thereafter.

The area of interhelical surface depends on the relative location of both helices. It changes, as the helices drift apart or elongate. It also changes when both helices bend in opposite directions as these amount to a distension of the surface. Whenever helices bend in similar ways in the same direction, interhelical area will be relatively unaffected. Depending on the complexity in shape of the helical axis, second-order polynomials might not adapt to the axis’ path accurately. Figure 9 shows an increasing trend of interhelical area for complexes B4403 and B4405, while a declining trend of interhelical area is seen for B4402. These changes occur as the helices of the peptide-binding groove move closer together or further apart, respectively. This is also reflected in the distance of α-helical centres of masses (data not shown). However, inspecting Figure 5 clearly demonstrates that the actual shape of an interhelical surface cannot be characterized by a single quantity such as the area of interhelical surface.

A more elaborate descriptor of the shape of a surface is the surface torsion. This parameter describes the change in angle between subsequent tangent planes along a given path, in our case the centre line. High values of surface torsion, regardless of being positive or negative, describe a rapid change in the angle. A pair of helices, each being somehow deformed and both being in varying positions to each other, gives rise to an interhelical surface with a plethora of possible shapes. For describing these shapes geometrically, surface torsion is an important concept, lending itself to quantify the twist in the surface along a prescribed path.

For B4402, surface torsion along the centre line (see Figure 5(a)) is left-handed most of the time, reaching a minimum (ca. −2, data not shown) at half of its length. Positive values of surface torsion (ca. 0.2, data not shown) are seen only in rare cases and near both surface termini. The ruled surface of the peptide-binding cleft at the N-terminus of G-ALPHA1 appears to be more stable for B4405 and B4403 than for B4402. Similar trends are seen in the surface torsion integral, see Figure 10. It is stable for B4405, but trends are observed for the other two complexes: B4402 shows a drop in the surface torsion integral during the first 60 ns followed by a gradual rise. B4403 is rather stable during the first 60–70 ns and then shows a gradual decrease in surface torsion.

We also looked at the shape complementarity () (Lawrence and Colman introduced a method to measure how well the surfaces of two proteins at their interface match [31]. The parameter output by this method is called shape complementarity ().) of the protein-protein interface surface and backbone RMSD. We found to be stable for all three TCR/pMHC systems (see Figure 14). As is unaffected by deformations of the binding cleft (as reflected by surface torsion, see above) we conclude that TCRs follow the conformational changes of the MHC surface. RMSDs of B4402 and B4403 reach a low plateau after a few nanoseconds and remain stable thereafter. RMSD of B4405 shows an increasing trend over 250 ns simulation time, rising from approximately 0.3 nm to 0.6 nm; see Figure 15.

4. Discussion

The mechanism of TCR activation is still controversially debated and several models have been proposed [3238]. A recent work of Dustin and Depoil [34] summarizes new insights into the function of the T cell synapse. The authors grouped T cell synapse into three interactive layers including interactions of receptors, a signaling layer, and a cytoskeleton layer, all contributing to TCR activation, regulation, and fine-tuning of signaling and responses. Conformational changes of the TCR complex have been demonstrated to be relevant for signaling of the TCR [39]. Association of CD3 proteins to the TCR/pMHC complex is necessary to transmit the activation signal to intracellular signaling molecules [40]. Evidence suggests that TCR conformational changes are required for full activation, but there are certain signaling pathways that can also be activated in the absence of conformational changes [41]. The three TCR/pMHC complexes analysed in this work differ only by one or two amino acids in the MHC molecule. 44:05 (B4405) and 44:02 (B4402) MHC types trigger LC13 TCR activation in the presence of the ABCD3 self-peptide. Surprisingly 44:03 (B4403) does not trigger TCR activation.

To characterise the dynamics of the MHC antigen binding cleft, we applied (similar to the methods reported by Christopher et al. [16]) a fragment-fitting method to model stable α-helical regions that are in close proximity (≤1.4 nm) to the TCR and monitored their geometric parameters. However, it is known that geometric quantities derived from polynomial approximations may vary substantially depending on the polynomial degree chosen [10]. To select an appropriate polynomial degree, we tested the ability of polynomials with different degrees to retrieve predefined parameters of helical motions and deformation. In a simplified model, we tested the ability of the fragment-fitting method to reproduce the curvature integral of helical bending and hinge motions. We found that second-order polynomials are best suited to model these α-helical motions with low relative error. The curvature integral derived from polynomials of individual α-helices can be related to conformational changes of α-helices. Between the two MHC α-helices a ruled surface can be spanned, of which we computed the area as an estimate for the size of the peptide-binding cleft. We also calculated the surface torsion along an imaginary centre line characterizing the orientation of both α-helices relative to each other. We applied this method to MD simulations of three TCR/pMHC complexes. However, we were not able to find correlations between immunogenicity and certain patterns in the α-helical movements.

4.1. Limitations

A limitation of the current analysis is that the ruled surface between the MHC α-helices that we use to model the MHC surface presented to the TCR does not consider the shape and dynamics of the peptide that lies between the two helices. We cannot assume that phase space has been sampled comprehensibly for these large molecules. Stepwise fluctuations in measured variables are visible (see Figure 7(b), upper row, and Figure 15). Also, even with highly optimized simulation performance of 15 ns/day on 1024 computing cores of the IBM BlueGene, statistics to discriminate between different simulated systems is not feasible. Enhanced sampling techniques and adequate collective variables might be useful to identify adequate collective variables for such systems. Interpretation of single simulations should therefore be done with caution. Signaling may involve a series of other proteins of the immunological synapse and interactions that are not considered in these simplified TCR/pMHC models. Interactions between TCR and pMHC take place between two cells that are in close contact to each other. It has been shown that plasma membrane lipids affect the activity of signaling networks [42] and some models of TCR/pMHC interaction propose involvement of the plasma membrane [33].

4.2. Conclusion

In this work, changes of MHC shape and their dynamics were quantified. We applied the quantification method to three large TCR/pMHC complexes, due to their size being accessible by MD simulation studies only since recently. We saw that MHC α-helices undergo rapid changes in conformation by either bending motion or hinge motion. Surface torsion used for characterizing the MHC surface presented to the TCR is stable in B4405, which is the most immunogenic complex. We speculate that rapid changes in helical conformation are part of the intrinsic dynamics of MHCs when engaging with TCRs.

Though we were not able to find a clear correlation between immunogenicity and certain patterns in the α-helical movements, we could demonstrate that single amino acid polymorphisms in the MHC seem to have a subtle influence on the helixes’ shape dynamics and that it would be interesting to apply the same method in the case of peptide polymorphism.

In summary, the present work demonstrates the feasibility and reliability of deriving shape parameters from simulation data. Next, the influence of the detected small conformational changes on the microscopic dynamics will be investigated to clarify their relation to the biological functions of the complexes of interest. Conclusions regarding functional differences between TCR/pMHC systems characterized by a single-residue polymorphism certainly require advanced sampling techniques in order to sample the conformational phase space appropriately for molecules of this size. Future studies might investigate if the small conformational changes in MHC α-helices transmit forces to the TCR.

Appendices

A. Construction of a Curved Helix Model

The backbone is modelled along a cosine curve. We compute equidistant points along the curve and create atoms with constant normal distances to the backbone. Angles between successive atoms were set to . Formally, coordinates along the axis of the helix are given bywhere represents the maximum elongation (amplitude) of the curved helix, compared to , corresponding to a straight model. Increasing in a stepwise fashion generates models of increasing curvature.

For the curvature we obtain In order to keep the length of backbones constant, the limits for parameter have to be adjusted appropriately; . In fact is chosen to make the arc length equal to 1: By varying within different models were created. Curvature decreases with decreasing ; in particular corresponds to a straight line without curvature. In order to find the appropriate values for with , an elliptical integral (see (A.3)) has to be solved. To find equidistant points along the backbone , for , we proceed, similar to the normalization of the arc length, via numerically solving the elliptical integral:In the next step, atoms are positioned in distance to the backbone and successively rotated by . The tangent vector of unit length at position is given byNext, the radius of α-helical turns around the axis is set in proper relation to the total length of the helix:and the vector to be rotated is given byThe point is rotated around the axis by angles to create coordinates of successive atoms . Finally, all coordinates are scaled via to arrive at proper dimensions. Thus, we obtain a chain of atoms rotated in a clockwise manner to represent a curved right-handed helix along a predefined, cosinusodial backbone.

The curved helix model has a few limitations. (i) Due to the curvature of the cosine arc and the construction of atoms with normal distances, the typical distance between successive atoms is not maintained. (ii) We refrained from constructing the full backbone also including nitrogen, oxygen, and hydrogen atoms and restricted the model to carbon atoms only. This seems justifiable, however, since the fragment-fitting algorithm also takes into account atoms only as does the DSSP algorithm and the calculation of close contacts. Therefore, modelling the full backbone would not add more information to the model. (iii) The cosine arc is planar and hence cannot incorporate any torsion within the curved helix model.

B. Additional Data

See Figures 1116.

Conflict of Interests

Reiner Ribarics is an employee and stockholder of Gilead Sciences. The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The MD trajectories used in the present work were generated on the IBM BlueGene/P computer facility at Bulgarian National Centre for Supercomputing Applications (NCSA, http://www.scc.acad.bg/). The work was supported in part by BSF and OeAD under Grants nos. DNTS-A 01-2/2013 and WTA-BG 06/2013.