Abstract

Molecular dynamics was used to simulate large molecules of the immune system (major histocompatibility complex class I, presented epitope, T-cell receptor, and a CD8 coreceptor.) To characterize the relative orientation and movements of domains local coordinate systems (based on principal component analysis) were generated and directional cosines and Euler angles were computed. As a most interesting result, we found that the presence of the coreceptor seems to influence the dynamics within the protein complex, in particular the relative movements of the two α-helices, Gα1 and Gα2.

1. Introduction

The interaction between major histocompatibility complexes (MHCs) and T-cell receptors (TCRs) plays a key role in triggering adaptive immune responses. TCRs bind the highly polymorphic MHC proteins that present peptide fragments (p) derived from the host proteome, pathogens, or tumour antigens on the cell surface. TCR and pMHC represent the core of the immunological synapse that in turn comprises many proteins, both membrane-bound and in the cytosol that could relay signals and/or act as an adjustable screw to fine-tune TCR sensitivity. Cluster of differentiation 8 (CD8) is a transmembrane, mostly heterodimeric glycoprotein that functions as a coreceptor for the TCR. It is mainly expressed by cytotoxic T-cells (), but is also found on natural killer cells, cortical thymocytes, and dendritic cells [1]. The extracellular domain of CD8 binds to the -domain of the MHC class I heavy chain [2]. It is well-known that CD8 and CD4 coreceptors are able to enhance T-cell responses to antigen stimulation [35]. Also, when subjected to an immune response, CD8+ T-cells can substantially increase in sensitivity by the mechanism of functional avidity maturation, that is, maturation of strength of multivalent antigen-antibody binding [68].

According to the literature, the major mechanism for stabilizing TCR-pMHC interaction by CD8 is the CD8-MHC interaction that increases the TCR-pMHC rebinding probability [9, 10]. A less obvious mechanism is stated by Borger et al. [11] and includes affected binding rates of TCR-pMHC. They propose a two-stage reversible reaction mechanism of pMHC with either TCR or CD8, similar to the mechanism found by Liu et al. [12].

Molecular dynamics (MD) simulations of TCR/pMHC (PDB ID: 3KPS) and TCR/pMHC (3KPS) plus CD8αα homodimer were performed. The topology of the TCR/pMHC complex with and without a CD8 coreceptor is shown in Table 1 and Figure 1. We chose to monitor the relative movements of MHC α-helices, Gα1 and Gα2, with and without the presence of CD8 (as these helices constitute part of the binding cleft for the peptide) and of the MHC -domain relative to the whole CD8 (since the -domain is the binding site for the CD8-coreceptor).

2. Methods

2.1. Molecular Modelling

The structure of LC13 TCR, ABCD3 peptide, and MHC (TCR/pMHC) of HLA-B44:05 type has been resolved (PDB-ID: 3KPS, www.pdb.org). Also, molecular structures of CD8 coreceptors bound to MHCs are available (PDB-ID: 1AKJ).

To our knowledge, a TCR/pMHC/CD8 complex has not yet been cocrystallized. To model the TCR/pMHC/CD8 complex we localized the CD8/MHC binding site in the 1AKJ crystal structure by finding all atoms of the MHC within the range of 0.8 nm to CD8. Structures of TCR/pMHC (3KPS) and MHC/CD8 (1AKJ) were merged into one file, both MHC binding sites superimposed so as to minimize RMSD in a least-squares sense and the MHC molecule from the MHC/CD8 complex deleted to get the TCR/pMHC/CD8 complex; see Figure 1.

2.2. Molecular Dynamics Simulations

Molecular dynamics simulations were performed with GROMACS 4.0.7 [13] using the gromos53a6 force field. The whole system counts about 274000 atoms, including the solvent (protein atoms only: about 8500), within a simulation box sized 13 × 13.5 × 16.5 nm3, to ensure 2 nm minimal distance between the protein atoms and the box walls, and with periodic boundary conditions imposed. The solvent was described with the SPC water model [14], the system neutralized at a salt concentration of 0.15 mol/L, and its energy was minimized by the steepest-descent method. The temperature was then gradually increased to 310 K within a 100 ps position-restraint simulation. Temperature was controlled by a Berendsen-thermostat with a time constant of 0.1 ps and the pressure controlled by a Berendsen-barostat set to 1 bar with a time constant of 0.5 ps, both chosen for being the most efficient in the beginning of the simulation. Constraints on all bonds were imposed with the LINCS algorithm [15], and the particle mesh Ewald (PME) method [16] was used to compute the long-range electrostatic interactions, with van der Waals and Coulomb cutoff radii of 1.4 nm. For the MD simulation runs of 200 ns with a time step of 5 fs, enabled by using virtual sites for hydrogen atoms, the thermostat was set to -rescale, with the same time constants in order to guarantee the generation of a proper canonical ensemble [13]. Coordinates were written every 50 ps, giving thus rise to 4000 frames. Prior to the analysis of domain movements, translational and rotational motions relative to the energy-minimized protein structure were removed.

2.3. Relative Location of Domains

Within the biomolecules we consider the Cα atoms in the backbones of protein chains. These Cα atoms are addressed via their indices to define domains or even subdomains; see Table 1. We define a first domain, , by enumerating the Cα atoms contained; for example, to select the α-helix Gα1 we have . Similarly, we define a second domain, . Note that domains may consist of several parts such as the β-sheet; see Table 1.

Considering a domain containing Cα atoms with Cartesian coordinates in MD-frame , the coordinates of its geometrical center are given by The distance between the centers of two domains, and , in MD-frame is thenBoth domains are shifted with their mean Cα coordinates into the origin before defining their relative orientation.

2.3.1. Rigid Axes within Deformable Domains

To quantify the relative orientation of two domains one has to bear in mind that each domain is deformed from time step to time step of an MD simulation and hence no unique frame of reference can easily be assigned as is possible for a rigid body. Often a principal component analysis (PCA) is performed to obtain three main characteristic axes of a given domain [1719]. However, principal components are not fully rigorously defined: for example, the orientation of the first PC-eigenvector may swap by nearly 180° for virtually the identical coordinates, just because of numerical noise. Similarly, the second and third PC-eigenvectors may interchange roles from time to time for an atomic domain almost cylindrical in shape.

In order to avoid such artifacts we refrained from computing principal components repeatedly for each MD time step and adopted the following procedure, see Figure 2:(1)For domain we select a reference frame, , from the whole trajectory. This is done by computing the sum of RMSD-displacements of relative to itself [20] over all frames of the trajectory and adopting the frame with minimum sum of RMSD. This frame is in a geometrical sense considered most “central” for domain within the whole trajectory.(2)The Cα coordinates of domain at frame are subjected to a PCA, yielding the orthogonal matrix of three orthonormal eigenvectors , , of the covariance matrix of the coordinates of Cα atoms in domain at frame . These vectors define a reference frame (i.e., a local coordinate system) for domain .(3)Steps (1) and (2) are performed also for domain , yielding the respective central frame with an eigenvector matrix and a local coordinate system for .

2.3.2. Computing Robust Relative Orientations

Given all MD-frames of a trajectory, the relative orientation of domains and is computed as follows:(a)For each frame we compute the transformation of all coordinates of domain from its position within the central frame into its position at frame according to minimum RMSD using Kabsch’s method [20].(b)The rotational part of the above transformation is applied to the eigenvectors of domain at frame (the reference frame) to obtain the position of the eigenvectors of at frame : .(c)Steps (a) and (b) are performed also for domain , yielding and transformed eigenvectors for each frame .(d)For each frame we note that the directional relation between the two sets of eigenvectors and can be represented via a rotation matrix aswhich also characterizes the relative orientation of both domains. From (3) we obtainsince the inverse of an orthogonal matrix equals its transpose. Rewriting the matrix in terms of its column vectors the Euler angles in -convention [21] between the two sets of the orthonormal eigenvectors and can be read aswith all quantities depending on frame f (dependencies suppressed in the notation).As a result, , , , and define the relative spatial relation between domains and for MD frame .

3. Results

3.1. Relative Movements CD8-MHC

As a first pair of domains we considered the CD8 homodimer (CD8 , CD8 ) as domain and domain of the MHC (see Table 1) as domain . Domain is the binding site for CD8 (see Figure 3) and therefore most interesting regarding relative movements.

3.1.1. Relative Distance

The distance between domains, (2), was computed over the time; see Figure 4. This distance, initially around 2.75 nm (as modelled, based on the crystallographic structure), gradually decreases to about 2.55 nm during the first 50 ns of the simulation. Apparently, dynamics lets CD8 get somewhat closer to the TCR/pMHC complex.

3.1.2. Relative Orientation

For each domain, the sum of RMSD of each frame to all other frames of the trajectory was computed to determine the central (reference) frames, and ; see Figure 5.

Relative positions of both domains were then computed for each frame (4000 all in all) as outlined above leading to the following results; see Section 2.3.2.

A first and direct measure is provided by the orientation cosines between corresponding eigenvectors; see Figure 6(a). They clearly reflect an initial phase different from the remaining trajectory, corresponding to the phase of decreasing interdomain distance, already seen in Figure 4. Interestingly, this effect is not at all revealed by the Euler angles; see Figure 6(b).

To further characterize the evolution of the geometry of the complex with time we computed the autocorrelation functions (ACFs) for cosines between eigenvectors and for Euler angles; see Figure 7(b). Euler angles exhibited relatively short autocorrelations, passing through zero already at approximately 20 ns, and oscillating around zero for larger time lags.

A completely different picture results from inspecting the ACFs of the cosines between corresponding eigenvectors; see Figure 7(a). While the relative orientation of the eigenvectors corresponding to the largest eigenvalues, and , has only a short memory (the ACF passes through zero around 20 ns), a massively prolonged memory is seen for both smaller eigenvectors versus and versus : they exhibit a long negative tail and do not become stochastic throughout the whole simulation time. This reflects the fact already seen in the time course itself (Figure 6(a)): a long equilibration phase extends up to about 100 ns (which amounts to half the simulation time).

3.2. Relative Movements of Two MHC α-Helices

The two α-helices of the MHC, Gα1 and Gα2, together with a β-floor form the binding cleft for the peptide. To evaluate their relative motion we chose the helices as domains and , located the central (reference) frames ( = 279 ≡ 138.4 ns and = 345 ≡ 171.2 ns) of each, and computed corresponding eigenvectors; see Figure 8(b). It is interesting to see that the first eigenvectors and have similar directions and orientations, when computed for the central frames. As opposed to this, they show almost opposite directions when computed from the first frame of the trajectory; see Figure 8(a). These opposite directions of eigenvectors are formally obtained from the very same matrix algebra, although the domains themselves have by no means turned upside down, as one can see by simple inspection. We have addressed this issue in Section 2 and display an example here. Our method of fitting reference domains via the Kabsch method has been designed to avoid these effects. In fact, the eigenvectors and we actually use for frame 1 (and for all other frames) have orientations similar to and , except for the actual, relative movements of both domains.

Relative motions of Gα1 and Gα2 are characterized by cosines (see Figure 9(a)) and Euler angles (Figure 9(b)) between eigenvectors attached to each domain. Cosines as well as Euler angles reveal a correlation in movements of eigenvectors 2 and 3. These eigenvectors point away from the axis of the helices at right angles, and their correlated changes indicate a synchronized oscillating “rolling” of both helices. In fact this kind of movement is also evident when visually inspecting the trajectories in VMD [22].

The above finding is nicely quantified via the correlation coefficient ρ = 0.79 (with and = 503); see Figure 10(a).

3.3. Impact of CD8 Presence

We have also analyzed the relative motions of Gα1 and Gα2 from our MD simulation of TCR/pMHC/CD8 with CD8 attached to the MHC and compared it with the above results without CD8. Interestingly, the relative motions do not differ significantly (figures not shown); however, the correlation of “rolling” oscillations is almost lost in the presence of CD8: correlation coefficient ρ = 0.6 (with and = 406); see Figure 10(b). To check if this difference in correlation coefficients is statistically significant, we computed the 95%-confidence intervals [23], resulting in for ρ = 0.79 ( = 503) and for ρ = 0.6 ( = 406). These intervals do not overlap and the difference in correlation coefficients may thus be considered statistically significant.

4. Discussion

The methodological parts of this work describe a new computational technique to obtain relative orientations of intramolecular domains. In the application parts this method is used to analyze the molecular dynamics of two systems, TCR/pMHC and TDC/pMHC/CD8, respectively.

4.1. Methods to Characterize Relative Orientations

There is a plethora of ways to characterize relative movements of intramolecular domains (e.g., [1719, 2427]). The most direct one is to compute average distances between groups of atoms, as they change over time. We did this by computing the distance between CD8 and MHC ; see Figure 4. By appropriate selection of target groups, some basic information can be obtained also on relative orientations.

In this work we present a more sophisticated approach by attaching local coordinate systems (of eigenvectors) to each domain and calculating the rotational relations between them. It is a well-known drawback of eigenvector-based techniques that the eigenvector orientation is not well defined and may suddenly switch into almost opposite directions. Up to now, this was in most cases mended by some logical condition in the code selecting the appropriate orientation with reference to some atoms that have to be individually specified. These drawbacks are even more severe when eigenvectors are computed from internally deformable sets of data points, such as MD-frames.

We have solved this problem by computing eigenvectors only once (for each domain) from a very specific frame (the central frame), see Section 2.3.1. In this first step, the orientation of both systems of eigenvectors may be corrected if desired. Thereafter, relations will remain stable without any intervention on the side of the researcher. Stability of local coordinate systems is achieved by fitting the atoms within each domain and carrying along the eigenvectors accordingly. This results in robust relative orientations.

Note that although the method seems computationally demanding, it is in fact fast. The reason is that the Kabsch algorithm is a direct matrix operation, not an iterative procedure, as one might expect. It requires no larger computational effort than singular value decomposition.

Once mutual orientations have been computed (rotation matrix), some attention was given to select suitable quantities to characterize relative orientations. We presented cosines between corresponding eigenvectors, since they can easily be interpreted. In addition, we used Euler angles as a well-known concept for relative orientations of rigid bodies.

4.2. Unifying Orientations

PCA as such yields eigenvectors unique in directions but ambiguous in orientation. For example, the first eigenvector may also result as , merely as a consequence of minute numerical noise in data substantially equal. The same is true for . The third eigenvector always has to satisfyand is hence determined by and . When comparing the relative orientation of two domains via their eigenvectors, orientation ambiguity has to be coped with, which could be done as follows:(i)After performing a PCA for the reference frame of domain , results are manually inspected and given a well-defined orientation within domain . This can be accomplished by selecting two Cα atoms and setting the orientation of such that a positive cosine of the angle between and the vector joining the two Cα atoms is obtained.(ii)The same is done for , with a second pair of appropriate Cα atoms selected. is computed via the cross product of and ; see above.In order to arrive at standardized eigenvectors allowing for comparison between different MD-runs the resulting eigenvectors should be reoriented (if necessary) after performing PCA for the reference frame of domain : a criterion for reorientation could be positive cosines with the eigenvectors of domain :Note that flipping the orientation of eigenvectors between different frames of a trajectory is avoided intrinsically by our procedure (fitting of domains rather than repeatedly performing PCA to successive frames). However, initially selecting appropriate and definite orientations is only guaranteed by the above precautions.

4.3. Application of the New Methods to CD8 Coreceptor Dynamics

Inspecting the distance between CD8 and MHC over time (Figure 4(a)) we clearly observe an equilibration phase pertaining up to 50 ns, during which distance decreases. This observation suggests that more extensive MD simulations than those presented in this work are necessary to reliably investigate equilibrium properties of this system. This finding is also supported by the autocorrelation function, which passes through zero around 60 ns and shows a long negative tail afterwards; see Figure 4(b).

The existence of an equilibration phase raises the importance of a comparison of the two presentation methods employed: cosines between eigenvectors and Euler angles (see Figure 6). Clearly, cosines reflect the fact that there is a relaxation phase, while the Euler angles do not. Interestingly, the eigenvectors corresponding to the main extensions (, ) fail to indicate the relaxation but (, ) and (, ) clearly do. This indicates that the relaxation is made up of some rotation around the major axes, and , respectively.

This surprising finding is supported by comparing the corresponding autocorrelation functions; see Figure 7. Euler angles (Figure 7(b)) lose memory already after 20 ns and the ACF then oscillates around zero. As opposed to this, cosines between eigenvectors corresponding to the smaller components ((, ) and (, )) show excessively long autocorrelations, pertaining throughout the whole simulation time. This finding, even more than the 50 ns equilibration of distance, indicates the necessity of additional MD simulations to arrive at a more adequate sampling.

As a most interesting result, we found that the presence of CD8 seems to influence the dynamics within the MHC, in particular the relative movements of the two α-helices, Gα1 and Gα2. In this case, Euler angles proved the more sensitive tool. Negative correlation between cosines of (, ) and (, ) was in both cases beyond ρ < −0.9, and presence or absence of CD8 did not make much difference. For the Euler angles, however, we obtained the nice reduction of movement-correlation induced by the presence of CD8; see Figure 10.

All in all, both sets of orientation parameters presented (cosines and Euler angles) seem to have their merits and weaknesses that have to be explored in many more situations to arrive at a comprehensive judgment.

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.

Authors’ Contribution

Wolfgang Schreiner and Rudolf Karch contributed equally.

Acknowledgments

MD simulations were run on the IBM-BlueGene/P at Bulgarian NCSA. Partial support by Bulgarian Science Fund and Austrian Academic Exchange Programme under Grants DNTS-A 01-2/2013 and WTA-BG 06/2013 is acknowledged. Karin Schlangen helped with the statistical analysis. The software for the analysis is available on request from the authors.