Research Article  Open Access
Geometry Dynamics of αHelices in Different Class I Major Histocompatibility Complexes
Abstract
MHC αhelices form the antigenbinding cleft and are of particular interest for immunological reactions. To monitor these helices in molecular dynamics simulations, we applied a parsimonious fragmentfitting method to trace the axes of the αhelices. Each resulting axis was fitted by polynomials in a leastsquares 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 secondorder 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 surfacebound 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/peptideMHC (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 surfacebound 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 MHCbound 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 antigenpresenting 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 [5–8]. 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 tyrosinebased 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 wellestablished, 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 selfpeptide (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 Xray structure of the ternary complex 44:05, ABCD3 nonapeptide, and LC13 TCR that is accessible on the protein database (PDB, http://www.pdb.org/, PDBID: 3KPS). Due to extensive sequence identity we were able to use in silico sitedirected mutagenesis to obtain 3D structures of the missing TCR/pMHC complexes. MD simulations in the nanosecond range could possibly show shortlived 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 peptidebinding 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 proteinprotein interface. Then we extracted the αhelical axis and finally determined a polynomial that approximates this axis in a leastsquares sense.
Various methods to extract a helix axis have been developed [16], including rotational fitting, using atoms as control points of Bsplines [17] or fitting to a helix. We used a fragmentfitting 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 leastsquares 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 peptidebinding cleft that is ligated to the TCR. Dynamics in the shape of the proteinprotein 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 Xray structure of TCR/pMHC B4405 (number 3 in Table 1) was taken from the PDB (PDBID: 3KPS). Structures B4402 and B4403 were engineered by means of in silico sitedirected 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.

(a)
(b)
(c)
(d)
(e)
The full amino acid sequence of HLAB 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 Xray structure as plasma membrane structures and flexible protein parts are hard to determine using Xray 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 cutoff of 1.4 nm. Longrange 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 ProteinProtein Interface
As outlined in the introduction, the MHC protein comprises αhelices and betasheets. 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 cutoff value of 0.5 is justified by the distinctly bimodal distribution (see Figure 11(a)).(ii)The relative presence of close contacts at the proteinprotein interface: atoms that are less than 1.4 nm apart (i.e., the cutoff for electrostatic interactions in our MD simulations) are defined as being in close contact. Atomatom 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 cutoff value of 0.5 is justified by the distinctly bimodal distribution (see Figure 11(b)). To get the residuewise relative contact time, we averaged the atomwise 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 cutoff set to 1.4 nm. Note that at the end of the MHC’s peptidebinding pocket both GALPHA1 and GALPHA2 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.
(a)
(b)
(c)
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 leastsquares 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 leastsquares 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 GALPHA1 and GALPHA2 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 leastsquare 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 welldefined 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.
(a)
(b)
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 righthanded surface torsion, indicates lefthanded 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 fragmentfitting 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 fragmentfitting method and the resulting helix axis was fitted by a polynomial in a leastsquares 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 fragmentfitting method. A polynomial of th order was fitted to the traced αhelical axis in a leastsquares 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 fragmentfitting method (Section 2.4), fitting the resulting axis by a polynomial in leastsquares 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 eighthdegree 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 secondorder 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 secondorder 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, GALPHA1 and GALPHA2, 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 seconddegree polynomials fitted to helical axes.
 
We compared helix conformations between different phases of the simulation, that is, before and after an inflection point or continous decrease/increase in the curvature integral (see Figure 7, upper row). Median values of the RMSD distribution are shown in this table. 
(a)
(b)
(c)
3.2. Geometric Quantities Characterizing the Shape of MHC Helices
The MHC peptidebinding groove comprises two αhelices that interact with the TCR. Using secondorder 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 peptidebinding 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. GALPHA1 of B4403 and B4405 each undergoes abrupt fluctuations in helical conformation, reflected in the curvature integral, but is stable in the phases inbetween. The curvature integral for GALPHA1 of B4402 shows a gradually increasing trend. We notice that the curvature integral for GALPHA2 is generally higher than that for GALPHA1, reflecting the kink near its Nterminal end. The two helical parts of GALPHA2 form a hinge. GALPHA2 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.
(a)
(b)
(a)
(b)
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, secondorder 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 peptidebinding 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 lefthanded 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 peptidebinding cleft at the Nterminus of GALPHA1 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.
(a)
(b)
(a)
(b)
(c)
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 proteinprotein 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 [32–38]. 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 finetuning 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 selfpeptide. 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 fragmentfitting 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 fragmentfitting method to reproduce the curvature integral of helical bending and hinge motions. We found that secondorder 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 peptidebinding 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 singleresidue 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 righthanded 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 fragmentfitting 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
(a)
(b)
(c)
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. DNTSA 012/2013 and WTABG 06/2013.
References
 O. Acuto, V. Di Bartolo, and F. Michel, “Tailoring Tcell receptor signals by proximal negative feedback mechanisms,” Nature Reviews Immunology, vol. 8, no. 9, pp. 699–712, 2008. View at: Publisher Site  Google Scholar
 D. N. Garboczi, P. Ghosh, U. Utz, Q. R. Fan, W. E. Biddison, and D. C. Wiley, “Structure of the complex between human Tcell receptor, viral peptide and HLAA2,” Nature, vol. 384, no. 6605, pp. 134–141, 1996. View at: Publisher Site  Google Scholar
 Y.H. Ding, K. J. Smith, D. N. Garboczi, U. Utz, W. E. Biddison, and D. C. Wiley, “Two human T cell receptors bind in a similar diagonal mode to the HLA A2/Tax peptide complex using different TCR amino acids,” Immunity, vol. 8, no. 4, pp. 403–411, 1998. View at: Publisher Site  Google Scholar
 J. K. Burkhardt, “Seeing is believing: sorting out signaling events at the immunological synapse,” The Journal of Immunology, vol. 194, no. 9, pp. 4059–4060, 2015. View at: Publisher Site  Google Scholar
 L. E. Samelson, M. D. Patel, A. M. Weissman, J. B. Harford, and R. D. Klausner, “Antigen activation of murine T cells induces tyrosine phosphorylation of a polypeptide associated with the T cell antigen receptor,” Cell, vol. 46, no. 7, pp. 1083–1090, 1986. View at: Publisher Site  Google Scholar
 G. R. Crabtree, “Contingent genetic regulatory events in T lymphocyte activation,” Science, vol. 243, no. 4889, pp. 355–361, 1989. View at: Publisher Site  Google Scholar
 L. A. Timmerman, N. A. Clipstone, S. N. Ho, J. P. Northrop, and G. R. Crabtree, “Rapid shuttling of NFAT in discrimination of Ca^{2+} signals and immunosuppression,” Nature, vol. 383, no. 6603, pp. 837–840, 1996. View at: Publisher Site  Google Scholar
 A. Weiss, R. Shields, M. Newton, B. Manger, and J. Imboden, “Ligandreceptor interactions required for commitment to the activation of the interleukin 2 gene,” Journal of Immunology, vol. 138, no. 7, pp. 2169–2176, 1987. View at: Google Scholar
 F. Pappalardo, V. Brusic, F. Castiglione, and C. Schönbach, “Computational and bioinformatics techniques for immunology,” BioMed Research International, vol. 2014, Article ID 263189, 2 pages, 2014. View at: Publisher Site  Google Scholar
 R. Ribarics, R. Karch, N. Ilieva, and W. Schreiner, “Geometric analysis of alloreactive HLA αhelices,” BioMed Research International, vol. 2014, Article ID 943186, 8 pages, 2014. View at: Publisher Site  Google Scholar
 M. Kenn, R. Ribarics, N. Ilieva, and W. Schreiner, “Finding semirigid domains in biomolecules by clustering pairdistance variations,” BioMed Research International, vol. 2014, Article ID 731325, 13 pages, 2014. View at: Publisher Site  Google Scholar
 W. A. Macdonald, Z. Chen, S. Gras et al., “T cell allorecognition via molecular mimicry,” Immunity, vol. 31, no. 6, pp. 897–908, 2009. View at: Publisher Site  Google Scholar
 W. Stacklies, M. C. Vega, M. Wilmanns, and F. Gräter, “Mechanical network in titin immunoglobulin from force distribution analysis,” PLoS Computational Biology, vol. 5, no. 3, Article ID e1000306, 2009. View at: Publisher Site  Google Scholar
 M. Karplus and J. Kuriyan, “Molecular dynamics and protein function,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 19, pp. 6679–6685, 2005. View at: Publisher Site  Google Scholar
 B. Hischenhuber, F. Frommlet, W. Schreiner, and B. Knapp, “MH2c: characterization of major histocompatibility αhelices—an information criterion approach,” Computer Physics Communications, vol. 183, no. 7, pp. 1481–1490, 2012. View at: Publisher Site  Google Scholar
 J. A. Christopher, R. Swanson, and T. O. Baldwin, “Algorithms for finding the axis of a helix: fast rotational and parametric leastsquares methods,” Computers and Chemistry, vol. 20, no. 3, pp. 339–345, 1996. View at: Publisher Site  Google Scholar
 J. A. Lopera, J. N. Sturgis, and J.P. Duneau, “Ptuba: a tool for the visualization of helix surfaces in proteins,” Journal of Molecular Graphics and Modelling, vol. 23, no. 4, pp. 305–315, 2005. View at: Publisher Site  Google Scholar
 B. Hischenhuber, H. Havlicek, J. Todoric, S. HöllriglBinder, W. Schreiner, and B. Knapp, “Differential geometric analysis of alterations in MH αhelices,” Journal of Computational Chemistry, vol. 34, no. 21, pp. 1862–1879, 2013. View at: Publisher Site  Google Scholar
 B. Hischenhuber, H. Havlicek, J. Todoric, S. HöllriglBinder, W. Schreiner, and B. Knapp, “Corrigendum: differential geometric analysis of alterations in MH alphahelices,” Journal of Computational Chemistry, vol. 34, no. 32, p. 2834, 2013. View at: Google Scholar
 P. K. Warme, F. A. Momany, S. V. Rumball, R. W. Tuttle, and H. A. Scheraga, “Computation of structures of homologous proteins. αlactalbumin from lysozyme,” Biochemistry, vol. 13, no. 4, pp. 768–782, 1974. View at: Publisher Site  Google Scholar
 C. Oostenbrink, A. Villa, A. E. Mark, and W. F. Van Gunsteren, “A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS forcefield parameter sets 53A5 and 53A6,” Journal of Computational Chemistry, vol. 25, no. 13, pp. 1656–1676, 2004. View at: Publisher Site  Google Scholar
 U. Omasits, B. Knapp, M. Neumann et al., “Analysis of key parameters for molecular dynamics of pMHC molecules,” Molecular Simulation, vol. 34, no. 8, pp. 781–793, 2008. View at: Publisher Site  Google Scholar
 W. Kabsch and C. Sander, “Dictionary of protein secondary structure: pattern recognition of hydrogenbonded and geometrical features,” Biopolymers, vol. 22, no. 12, pp. 2577–2637, 1983. View at: Publisher Site  Google Scholar
 L. Pauling, R. B. Corey, and H. R. Branson, “The structure of proteins: two hydrogenbonded helical configurations of the polypeptide chain,” Proceedings of the National Academy of Sciences of the United States of America, vol. 37, no. 4, pp. 205–211, 1951. View at: Publisher Site  Google Scholar
 W. Kabsch, “A solution for the best rotation to relate two sets of vectors,” Acta Crystallographica A, vol. 32, no. 5, pp. 922–923, 1976. View at: Publisher Site  Google Scholar
 C. Runge, “Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten,” Zeitschrift für Mathematik und Physik, vol. 46, no. 1, pp. 224–243, 1976. View at: Google Scholar
 J. P. Boyd and F. Xu, “Divergence (Runge Phenomenon) for leastsquares polynomial approximation on an equispaced grid and MockChebyshev subset interpolation,” Applied Mathematics and Computation, vol. 210, no. 1, pp. 158–168, 2009. View at: Publisher Site  Google Scholar
 M. Peternell, H. Pottmann, and B. Ravani, “On the computational geometry of ruled surfaces,” CAD Computer Aided Design, vol. 31, no. 1, pp. 17–32, 1999. View at: Publisher Site  Google Scholar
 W. Kühnel, Differentialgeometrie Kurven—Flächen—Mannigfaltigkeiten, Vieweg+Teubner, Wiesbaden, Germany, 2008.
 P. W. Verbeek and L. J. van Vliet, “Curvature and bending energy in digitized 2D and 3D images,” in Proceedings of the 8th Scandinavian Conference on Image Analysis, pp. 1403–1410, Tromsø, Norway, May 1993. View at: Google Scholar
 M. C. Lawrence and P. M. Colman, “Shape complementarity at protein/protein interfaces,” Journal of Molecular Biology, vol. 234, no. 4, pp. 946–950, 1993. View at: Publisher Site  Google Scholar
 D. Aivazian and L. J. Stern, “Phosphorylation of T cell receptor zeta is regulated by a lipid dependent folding transition,” Nature Structural Biology, vol. 7, no. 11, pp. 1023–1026, 2000. View at: Publisher Site  Google Scholar
 K. Choudhuri and P. A. van der Merwe, “Molecular mechanisms involved in T cell receptor triggering,” Seminars in Immunology, vol. 19, no. 4, pp. 255–261, 2007. View at: Publisher Site  Google Scholar
 M. L. Dustin and D. Depoil, “New insights into the T cell synapse from single molecule techniques,” Nature Reviews Immunology, vol. 11, no. 10, pp. 672–684, 2011. View at: Publisher Site  Google Scholar
 S. T. Kim, K. Takeuchi, Z.Y. J. Sun et al., “The αβ T cell receptor is an anisotropic mechanosensor,” Journal of Biological Chemistry, vol. 284, no. 45, pp. 31028–31037, 2009. View at: Publisher Site  Google Scholar
 C. Xu, E. Gagnon, M. E. Call et al., “Regulation of T cell receptor activation by dynamic membrane binding of the CD3epsilon cytoplasmic tyrosinebased motif,” Cell, vol. 135, no. 4, pp. 702–713, 2008. View at: Publisher Site  Google Scholar
 M. S. Kuhns, A. T. Girvin, L. O. Klein et al., “Evidence for a functional sidedness to the αβTCR,” Proceedings of the National Academy of Sciences of the United States of America, vol. 107, no. 11, pp. 5094–5099, 2010. View at: Publisher Site  Google Scholar
 H. Zhang, S.P. Cordoba, O. Dushek, and P. A. van der Merwe, “Basic residues in the Tcell receptor zeta cytoplasmic domain mediate membrane association and modulate signaling,” Proceedings of the National Academy of Sciences of the United States of America, vol. 108, no. 48, pp. 19323–19328, 2011. View at: Publisher Site  Google Scholar
 N. MartinezMartin, R. M. Risueno, A. Morreale et al., “Cooperativity between T cell receptor complexes revealed by conformational mutants of CD3ε,” Science Signaling, vol. 2, no. 83, article ra43, 2009. View at: Publisher Site  Google Scholar
 G. Ryan, “T cell signalling: CD3 conformation is crucial for signalling,” Nature Reviews Immunology, vol. 10, no. 1, p. 7, 2010. View at: Publisher Site  Google Scholar
 R. Blanco, A. Borroto, W. Schamel, P. Pereira, and B. Alarcon, “Conformational changes in the T cell receptor differentially determine T cell subset development in mice,” Science Signaling, vol. 7, no. 354, Article ID ra115, 2014. View at: Publisher Site  Google Scholar
 A. A. Petruk, S. Varriale, M. R. Coscia, L. Mazzarella, A. Merlino, and U. Oreste, “The structure of the CD3 ζζ transmembrane dimer in POPC and raftlike lipid bilayer: a molecular dynamics study,” Biochimica et Biophysica Acta (BBA)—Biomembranes, vol. 1828, no. 11, pp. 2637–2645, 2013. View at: Publisher Site  Google Scholar
 W. Humphrey, A. Dalke, and K. Schulten, “VMD: visual molecular dynamics,” Journal of Molecular Graphics, vol. 14, no. 1, pp. 33–38, 1996. View at: Publisher Site  Google Scholar
 M. L. Connolly, “Analytical molecular surface calculation,” Journal of Applied Crystallography, vol. 16, no. 5, pp. 548–558, 1983. View at: Publisher Site  Google Scholar
 W. Schreiner, R. Karch, B. Knapp, and N. Ilieva, “Relaxation estimation of RMSD in molecular dynamics immunosimulations,” Computational and Mathematical Methods in Medicine, vol. 2012, Article ID 173521, 9 pages, 2012. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2015 Reiner Ribarics et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.