Relative Movements of Domains in Large Molecules of the Immune System

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.


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 (T C ), but is also found on natural killer cells, cortical thymocytes, and dendritic cells [1]. The extracellular domain of CD8 binds to the 3 -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 [3][4][5]. 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 antigenantibody binding [6][7][8].
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 3domain relative to the whole CD8 (since the 3 -domain is the binding site for the CD8-coreceptor).

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 nm 3 , 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 Berendsenthermostat 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 v-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.

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 = {59, 60, . . . , 83}. 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 x ( ) 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 then Both domains are shifted with their mean C coordinates into the origin before defining their relative orientation. Figure 2: Relative orientation of two submolecular domains.
Orthonormal eigenvectors for domains V and and standard basis for the laboratory system L. Rotation matrix R transforms eigenvectors from domain into domain . Note that eigenvectors share the same coordinate system origin. For better visualization, eigenvectors are displayed as if they had different origins.

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 [17][18][19]. 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 T ( ) = [k 1 k 2 k 3 ] of three orthonormal eigenvectors k 1 , k 2 , k 3 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 T ( ) = [w 1 w 2 w 3 ] and a local coordinate system for . (a) For each frame we compute the transformation of all coordinates x ( ) 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 R ( ) 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 : which also characterizes the relative orientation of both domains. From (3) we obtain since the inverse of an orthogonal matrix equals its transpose. Rewriting the matrix R ( ) in terms of its column vectors the Euler angles in -convention [21] between the two sets of the orthonormal eigenvectors [k 1 k 2 k 3 ] and [w 1 w 2 w 3 ] can be read as with 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 . Domain CD8 is shown in cyan, 3 in purple, and the remaining parts of the MHC as well as the TCR (labelled "protein") in black. Eigenvectors (k 1 , k 2 , k 3 and w 1 , w 2 , w 3 , resp.) are shown for the first frame of the trajectory and colored (k 1 , w 1 : blue, k 2 , w 2 : red, k 3 , w 3 : yellow) for each of the domains.

Relative Movements CD8-MHC.
As a first pair of domains we considered the CD8 homodimer (CD8 1 , CD8 2 ) as domain and domain 3 of the MHC (see Table 1) as domain . Domain 3 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.

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, k 1 and w 1 , has only a short memory (the ACF passes through zero around 20 ns), a massively prolonged memory is seen for both smaller eigenvectors k 2 versus w 2 and k 3 versus w 3 : 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).

Relative Movements of Two MHC -Helices.
The twohelices 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 k 1 ( ) and w 1 ( ) 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 k 1 (1) and w 1 (1) we actually use for frame 1 (and for all other frames) have orientations similar to k 1 ( ) and w 1 ( ), 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 < 0.01 and frames = 503); see Figure 10(a).

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 < 0.01 and f rames = 406); see Figure 10(b). To check if this difference in correlation coefficients is statistically significant, we computed the 95%-confidence intervals [23]   ( frames = 503) and [0.553, 0.659] for = 0.6 ( f rames = 406). These intervals do not overlap and the difference in correlation coefficients may thus be considered statistically significant.

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.

Methods to Characterize Relative Orientations.
There is a plethora of ways to characterize relative movements of intramolecular domains (e.g., [17][18][19][24][25][26][27]). 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 3 ; 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.

Unifying
Orientations. PCA as such yields eigenvectors unique in directions but ambiguous in orientation.
For example, the first eigenvector k 1 may also result as k * 1 = −k 1 , merely as a consequence of minute numerical noise in data substantially equal. The same is true for k 2 . The third eigenvector always has to satisfy k 3 = k 1 × k 2 (7) and is hence determined by k 1 and k 2 . 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 k 1 given a well-defined orientation within domain . This can be accomplished by selecting two C atoms and setting the orientation of k 1 such that a positive cosine of the angle between k 1 and the vector joining the two C atoms is obtained.
(ii) The same is done for k 2 , with a second pair of appropriate C atoms selected. k 3 is computed via the cross product of k 1 and k 2 ; 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