Using molecular dynamics simulations to interrogate T cell receptor non-equilibrium kinetics.

An atomic-scale mechanism of T Cell Receptor (TCR) mechanosensing of peptides in the binding groove of the peptide-major histocompatibility complex (pMHC) may inform the design of novel TCRs for immunotherapies. Using steered molecular dynamics simulations, our study demonstrates that mutations to peptides in the binding groove of the pMHC - which are known to discretely alter the T cell response to an antigen - alter the MHC conformation at equilibrium. This subsequently impacts the overall strength (duration and length) of the TCR-pMHC bond under constant load. Moreover, physiochemical features of the TCR-pMHC dynamic bond strength, such as hydrogen bonds and Lennard-Jones contacts, correlate with the immunogenic response elicited by the specific peptide in the MHC groove. Thus, formation of transient TCR-pMHC bonds is characteristic of immunogenic peptides, and steered molecular dynamics simulations can be used in the overall design strategy of TCRs for immunotherapies.

Introduction

T cells discriminate peptide-major histocompatibility complexes (pMHCs) expressed on antigen presenting cells (APCs) via their T Cell Receptors (TCRs). Engineered T cells have elicited exciting clinical responses in cancer patients [1], [2], [3], [4], [5]; however, matching TCRs with specific pMHC targets to invoke an appropriate and targeted immune response remains a challenge, and has created a need for a deeper understanding of how the TCR engages its antigen to induce an immunogenic response.

Several features of the TCR-pMHC interaction have been evaluated in an effort to identify quantitative descriptors that predict in vitro T cell immunogenicity [6], [7], [8], [9], [10], [11]. For example, equilibrium TCR-pMHC kinetic parameters, measured by surface plasmon resonance, have generated a collection of theories to explain specificity including kinetic proofreading [12], [13], [14], [15], serial triggering [16], [17], [18], and the law of mass action [19], [20], [21], [22]. However, none of these theories can completely predict the immunogenic response across multiple TCR-pMHC systems [23], [24], [25], [26], indicating that equilibrium-based models are insufficient.

In vivo binding of a cytotoxic T cell to an APC is a dynamic process that introduces mechanical forces on the TCR-pMHC interaction [27], [28] indicating that these interactions may not occur in equilibrium. Using DNA-based tension probes, mechanical forces on the TCR-pMHC bond have been estimated in the 10–20 pN range [29], [30]. These tensile forces may be attributed to interstitial shear on cells [31], relative membrane sliding during cell motility [32], or actin polymerization during formation of the immunological synapse [33], [34]. To investigate the effects of these forces, several biophysical techniques — including the biomembrane force probe [35], [36], [37], optical tweezers [38], and the laminar flow chamber [39], [40] — have been used to characterize the force-dependent TCR-pMHC dissociation kinetics. These studies clearly demonstrate that dissociation kinetics at critical forces (∼10–20 pN) correlate with immunogenicity [35], [36], [37], [38], [39], [40], [41], highlighting the importance of TCR-pMHC bond strength, and suggest a force-dependent kinetic proofreading discrimination model. In such a model, the TCR would sustain and form transient bonds under load for sufficient time (i.e., stabilized) to initiate biochemical signalling. While intriguing, the underlying molecular level mechanisms for bond lifetime enhancement under load are unknown.

Emergence of mechanosensing by the T cell underscores the importance of the dynamic energy landscape and molecular mechanics of the TCR-pMHC interaction under tensile force. Steered molecular dynamics (SMD) simulations have been used to investigate these properties [36], [37]. These reports have argued that transiently formed Lennard-Jones contacts (LJ-Contacts – defined as a distance less than 0.35 nm between atoms) or hydrogen bonds (H-Bonds) are responsible for increased bond lifetime. In contrast, a recent study [42] argues that the constant regions of the TCR stabilize and preserve the interfacial pMHC H-Bonds and LJ-Contacts under physiological load (∼10–20 pN). However, these studies did not provide a comprehensive analysis of the TCR-pMHC bond under load at the atomic scale. This level of examination is necessary to decipher the relative energetic contributions to and locations of H-Bonds and LJ-Contacts in the overall TCR-pMHC bond. Finally, atomic fluctuations in receptor-ligand interactions result in a stochastic ensemble of equilibrium configurations and the consequential heterogeneity of SMD simulation results remains unaddressed [43].

To initiate a more detailed investigation into the atomic level interactions between the TCR and pMHC, we began with a TCR-pMHC complex with a known crystal structure – DMF5 TCR bound to the MART1(AAGIGILTV) restricted pMHC (PDB ID: 3QDJ). We performed SMD in triplicate to examine the stochastic ensemble of equilibrium configurations. We then investigated three specific peptides, L1 (LAGIGILTV) and GVA (GAGIGVLTA, b-endoxylanase), restricted to the pMHC, which have been reported as a super-agonist (increased immunogenic response compared to wildtype) [44] and antagonist (no immunogenic response) [45], respectively (Fig. 1, Supplementary Videos). We found that the peptide identity alters the conformation of the MHC binding groove at equilibrium, and these conformational changes result in characteristic sets of H-Bonds and LJ-Contacts – at equilibrium and during dynamic pulling – that significantly contribute to the overall dynamic bond strength under load. These sets of H-Bonds and LJ- Contacts are uniquely distributed across the pMHC interface, occur at characteristic time and distance scales, and correlate with the immunogenicity of the specific peptide. Furthermore, under load, all TCR-pMHCs maintain their secondary structure and essential atomic motion is predominantly in the direction of applied force. These results provide new atomic-level insight on TCR-pMHC bond strength that may inform TCR design strategies.

Fig. 1

Graphic representation of Steered Molecular Dynamics (SMD) simulations on a TCR-pMHC system. This includes highlighted interfacial substructures [MHC⍺ & MHCβ = green, TCR CDR⍺ = blue, TCR CDRβ = red, Peptide (L1, MART1, GVA) = purple, magenta, orange], direction of pulling (yellow arrows), and position of (yellow circles) and distance between COMs (scale bar at bottom). The non-interacting bodies of the TCR and pMHC are colored in gray. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Results & discussion

Equilibrated structures and interaction energy

We initiated our simulation with the crystal structure of the TCR-pMHC complex bound to the MART1 peptide (PDB ID: 3QDJ) [46]; biomolecular simulations of systems this size require tens of nanoseconds to reach equilibration [47], [48]. TCR-pMHC complexes bound to the L1, MART1, and GVA peptides were simulated for 100 ns and reached near equilibrated configurations around 75, 20, and 40 ns, respectively (Fig. S3 A-C). The equilibrated configurations are indicated by the flattening of the root mean square deviation plots with fluctuations<0.10 nm for all protein chains (Fig. S3 A-C). Additionally, the GROMOS clustering algorithm [49] was used to identify structural similarity. In the equilibration time window with a C⍺ RMSD cutoff of 0.20 nm, approximately 56%, 94%, and 94% of structures were in the top ten clusters for L1, MART1, and GVA, respectively (Fig. S3 D-F). Moreover, in the 90 to 100 ns time window, the root mean square fluctuations range from 0.05 − 0.10 nm for all TCR-pMHC interfacial substructures (Fig. 2A), consistent with equilibration. Although there are statistically significant differences in TCR substructure fluctuations (i.e., CDR1⍺, CDR2⍺, CDR3⍺, and CDR3β) (Fig. S4) in this time window, there is no observable trend, and the magnitude of the differences is very small (0.02 nm). Nonetheless, to account for the ensemble of protein configurations and capture the average behavior, we performed pull simulations in triplicate utilizing the equilibrated structures at 90, 95, and 100 ns for all three mutants (Fig. 2B-D). After structural equilibration, the sampled configurations (90, 95, 100 ns) for SMD were in the two most dominant clusters for L1, MART1, and GVA (Fig. S3 D-F).

Fig. 2

Equilibrated TCR-pMHC Structures from 100, 95, and 90 ns included with Interfacial Substructures. (A) Interfacial substructure root mean square fluctuations for each TCR-pMHC from 90 to 100 ns. The grouped bars are in order (L1, MART1, GVA) for each respective interfacial substructure. (B) L1 peptide (purple) and TCR-pMHC (C) MART1 peptide (magenta) and TCR-pMHC (D)GVA peptide (orange) and TCR-pMHC (E) Aligned pMHC structure at 100, 95, and 90 ns for each peptide (L1 = purple, MART1 = magenta, GVA = orange). The equilibrated MHC⍺ & MHCβ are indicated in the peptide color to distinguish structure (previously green in B-D). Substructures of mutants were statistically compared (n = 3): #p < 0.10, *p < 0.05, **p < 0.01, ***p < 0.001 by one-way ANOVA followed by Tukey-HSD post-hoc test. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Interestingly, mutations to the peptide determine the conformation of the surrounding MHC ⍺-helices (Fig. 2E) and TCR substructures (Fig. S4) and are independent of fluctuations at equilibrium. To broadly assess TCR-pMHC bond strength, we evaluated interaction energies, H-Bonds, and LJ-Contacts on center of mass (COM) distance and time scales. Error bars represent the standard error of measurement (SEM) for each peptide (L1, MART1, and GVA) over the three SMD simulations. In accordance with immunogenicity, longer reaction distance (Fig. 3A) and non-linear energetic decay (Fig. 3A-B) for L1 and MART1 is consistent with higher energetic bond strength than GVA. The interaction energy between the TCR-MHC (Fig. 3A) is an order of magnitude larger than the TCR-peptide (Fig. 3B). Thus, peptide-induced structural alterations to the MHC ⍺-helices (Fig. 2E) determine the interaction energy (and reaction distance) of the TCR-MHC (Fig. 3A) for L1, MART1, and GVA. As such, the unique sequence of the peptide alters the equilibrated structure of the MHC and the subsequent TCR-MHC interaction, thereby functioning as the key to unleash the immunogenic response. This concurs with experimental observations of peptide-altered MHC binding grooves [50] and peptide-altered immunogenic responses [44], [45].

Fig. 3

TCR-pMHC Interaction Energy. (A) TCR-MHC and (B) TCR-Peptide interaction energy is separated into Lennard-Jones and Coulombic potential for each peptide and plotted against the reaction coordinate defined as the TCR-pMHC COM distance. For TCR-pMHC SMD simulations (100, 95, 90 ns), interaction energies (with 1 ps resolution) were placed into ∼0.5 Å bins and error represents SEM.

For receptor-ligand interactions, Coulombic and Lennard-Jones potentials are dominated by H-Bonds and hydrophobic LJ-Contacts, respectively. The TCR-MHC Coulombic potential is higher in magnitude than the Lennard-Jones potential for all mutants along the reaction coordinate (Fig. 3A). Conversely, the TCR-Peptide Lennard-Jones potential is stronger than the Coulombic potential for all mutants along the reaction coordinate (Fig. 3B). Thus, these results indicate that there is localized energetic heterogeneity, where H-Bonds are more consequential for TCR-MHC interactions and LJ-Contacts are more consequential for TCR-Peptide interactions.

Hydrogen bonds, Lennard-Jones Contacts, and solvent accessible surface area

In order to quantify the magnitude and duration of transient H-Bonds or LJ-Contacts during the pull, we utilized probability densities and existence occupancies. Probability densities represent the likelihood that a variable has a specific value during the pull simulation and was applied to H-Bonds and peptide solvent accessible surface area (SASA). For existence occupancy, we used cutoff values that represent the percentage of time under load that unique H-Bonds or LJ-Contacts exist; thus, a higher occupancy represents H-Bonds or LJ-Contacts that were present during a longer portion of the pull simulation and therefore represents more stable transient interactions.

A trend of increasing total H-Bonds from the probability density distribution (Fig. 4A) for GVA, MART1, and L1 (4.47, 5.04, and 5.99, respectively) combined with less low occupancy (less stable) H-Bonds for GVA compared to MART1 (Fig. 4B) suggest that transiently formed H-Bonds – bonds that exist for very short periods of time – may impact TCR-pMHC bond strength. The trend of increasing H-Bonds with increased peptide immunogenicity continues for more high occupancy (more stable, up to 10%) H-Bonds suggesting that more stabilized interactions may contribute to immunogenicity (Fig. 4B); however, there is no statistical significance. No statistical significance of high occupancy H-Bonds between mutants is likely the result of one heterogenous sample for L1 and MART1 in the stochastic ensemble (n = 3) demonstrating the importance of evaluating the stochastic fluctuations at equilibrium by performing several independent simulations. Furthermore, there is qualitative agreement between the Coulombic potential reaction distance (Fig. 3A-B) and the number of H-Bonds for TCR-MHC (Fig. 4B) interactions.

Fig. 4

Hydrogen Bonds and Lennard-Jones Contacts. (A) For each peptide, probability density of H-Bonds between TCR and pMHC is plotted. For TCR-pMHC SMD simulations (100, 95, 90 ns), times (with 1 ps resolution) were distributed into number of H-Bonds bins and error represents SEM. (B) Number of unique H-Bonds are graphed with increasing existence occupancy between the TCR and pMHC. The interactions are separated into TCR-MHC H-Bonds (green) and TCR-Peptide H-Bonds (L1 = purple, MART1 = magenta, GVA = orange) and stacked. Error represents SEM over three SMD simulations. (C) For each peptide, probability density of solvent accessible surface area (SASA) is plotted. For TCR-pMHC SMD simulations (100, 95, 90 ns), times (with 1 ps resolution) were cut into bins and error represents SEM. (D) Number of unique LJ-Contacts are graphed with increasing existence occupancy between the TCR and pMHC. The interactions are separated into TCR-MHC LJ-Contacts (green) and TCR-Peptide LJ-Contacts (L1 = purple, MART1 = magenta, GVA = orange) and stacked. Error represents SEM over three SMD simulations. Mutants were statistically compared (n = 3): #p < 0.10, *p < 0.05, **p < 0.01, ***p < 0.001 by one-way ANOVA followed by Tukey-HSD post-hoc test.. Mutants at >0% existence occupancy are statistically compared (n = 3) with a student’s t-test and displayed in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In previous work with a different system, we have demonstrated that increased hydrophobic SASA yields greater aggregation propensity [51] and likely more hydrophobic interactions. Significantly higher mean probability density of peptide hydrophobic SASA for L1 and MART1 suggest more hydrophobic interactions/LJ-Contacts (Fig. 4C). More numerous low (>0%) and high occupancy (up to 25%) LJ-Contacts for L1 compared to GVA highlight the importance of transient and stabilized hydrophobic interactions (Fig. 4D). Moreover, there is qualitative agreement between the Lennard Jones potential reaction distance (Fig. 3A-B) and the number of LJ-Contacts for TCR-Peptide (Fig. 4D) interactions. These results suggest that TCR-pMHC bond strength is impacted by stabilized and transiently formed H-Bonds and LJ-Contacts. In contrast to previous reports that claim either transient or stabilized interactions are responsible for increased bond lifetime [36], [37], [42] our atomic-level simulations present evidence that a single interaction type cannot explain the observed behavior. Importantly, the generalizability of these results is limited to the TCR-pMHC interaction presented and future investigations may reveal alternate trends related to the energetic importance for disparate peptides and MHC alleles (e.g., peptides with charged residues). However, these results provide clarity that there is energetic heterogeneity at the TCR-pMHC interface, this energetic heterogeneity effects the overall strength of interaction, and the overall strength of interaction corresponds with biological function.

TCR-pMHC dissociation and essential atomic motion

To investigate whether conformational changes in the protein complexes occur during the pull and contribute to the interaction energy landscape, we investigated essential atomic motion. During the equilibrated state and under load, root mean square fluctuations do not exceed 0.10 nm for any TCR-pMHC interfacial substructure (Fig. 5). Although there are statistically significant differences in substructure fluctuations between equilibration and pull simulations (i.e., CDR3⍺, CDR1β, CDR3β, MHC⍺, MHCβ, and Peptide), there is no trend except for an increase in the RMSF for the peptide during the pull. This slight increase for the peptide is perhaps expected because, as the TCR is dissociated from the pMHC, there is more space for the peptide to fluctuate. Moreover, 0.02 nm differences in fluctuation are unlikely to have physical significance because these fluctuations are less than the global fluctuations (<0.10 nm) of the TCR-pMHC structure at equilibrium [52]. Additionally, secondary structure analysis of the TCR-pMHC interface under load reveals that all complexes maintain secondary substructures through the bond dissociation process (Fig. S6). Moreover, a principal component analysis on the TCR-pMHC backbone indicates that more than 96% of essential atomic motion is in the direction of the applied force (Fig. S7, Supplementary Videos) except for the L1 95 ns and L1 90 ns configurations. Principal component analysis on atomic motion indicates that for these L1 configurations, the TCR-pMHC interface is more energetically stabilized under load than the quaternary structure of the TCR or pMHC, respectively (Fig. S7, Supplementary Videos). This energetic stability results in longer bond lifetime (Fig. S5), longer reaction distance (Fig. 3), and more transient interactions (Fig. 4). Comparable interfacial fluctuations (Fig. 5) demonstrate that differences in bond strength are mostly attributed to differences in the equilibrated structures, and not conformational changes to the protein structure during the pull. This is supported by the evidence that the peptide alters the structure of the MHC ⍺-helices at equilibrium, which determines TCR-pMHC energetic bond strength (Fig. 3).

Fig. 5

Equilibrated (EQ) and SMD (Pull) interfacial substructure root mean square fluctuations (RMSF) for each TCR-pMHC. EQ error represents SEM on substructure atoms from 90 to 100 ns. The grouped bars are in order (L1, MART1, GVA) for each respective interfacial substructure. SMD error represents SEM on substructure atoms from 3 independent SMD simulations (100, 95, and 90 ns starting configurations). Equilibration and pull subsctructures of mutants were statisically compared: #p < 0.10, *p < 0.05, **p < 0.01, ***p < 0.001 by one-way ANOVA followed by Tukey-HSD post-hoc test.

TCR-pMHC signature interaction maps

To characterize the alterations more precisely, existence occupancy heat maps were created for H-Bonds and LJ-Contacts. To identify the most significant H-Bonds, the hydrogen bonding heat maps included all H-Bonds with a minimum of 5% existence occupancy and located in at least one of the nine (3 peptides and 3 simulations for each peptide) SMD simulations. For the more numerous LJ-Contacts, the heat maps include LJ-Contacts with a minimum of 80% occupancy and are found in at least one of the nine SMD simulations. Although the existence occupancy cutoff (5% or 80%, respectively) is somewhat arbitrary, this confines the heat maps to ∼40–50 of the most stable interactions which may provide TCR mutagenesis targets. Furthermore, TCR-MHC and TCR-Peptide interactions were separated to provide the ability to independently examine TCR-MHC and TCR-Peptide signatures.

To evaluate the precise H-Bond or LJ-Contact existence occupancy as a function of distance, the time axis was converted to COM distance by distributing results into 0.5 Å bins and calculating the fractional occupancy in each respective bin. These existence occupancy heat maps as a function of distance were averaged across the triplicate SMD simulations revealing an overall image of unique H-Bond and LJ-Contact signatures (Figs. 6 A-C and 7 A-C, respectively) for L1, MART1, and GVA. The identical TCR-MHC y-axis (Fig. S8 A-B) facilitates direct comparison of the existence heat maps between the three peptides and between H-Bonds and LJ-Contacts and provides insight into specific interactions that may be of particular importance in T cell activation. For example, the two H-Bonds donated from TCR⍺ Arginine 27 to MHC Glutamate 55 are absent for GVA, present through reaction coordinate ∼7.50 nm for MART1, and present through reaction coordinate ∼8.30 nm and at greater occupancy (more “heat”) for L1 (Fig. 6). Similarly, the LJ-Contacts between TCR⍺ Tyrosine 50 to MHC Glutamine 155 are absent for GVA, present through reaction coordinate ∼7.14 nm for MART1, and present at greater occupancy (more “heat”) through reaction coordinate ∼6.82 nm for L1 (Fig. 7). Because L1 and MART1 are immunogenic peptides, these interactions on the DMF5 TCR in the CDR1⍺ and CDR2⍺ regions, respectively, may indicate recognition hot spots. Moreover, the energetic importance (Fig. 3B), significantly larger peptide SASA (Fig. 4C), and number of LJ-Contacts for the L1 peptide (Fig. 4D) are consistent with the nine contacts with over 80% existence occupancy compared to zero for MART1 and GVA (Fig. 7).

Fig. 6

Average hydrogen bond heat maps as function of TCR-pMHC COM distance. The time axis of the H-Bond existence maps is converted to COM distance by distributing time points into ∼0.5 Å bins and calculating the fractional occupancy in each respective bin. The 100, 95, and 90 ns configuration maps of each peptide are arithmetically averaged. The fractional occupancy is represented by the heat scale on the y-axis (right). Peptides (A) L1, (B) MART1, and (C) GVA are comprised of three SMD simulations. The H-Bond acceptor and donor index is specified on the y-axis (left). The index corresponds to donor-acceptor pairs and is in Fig. S8. In each panel, H-Bonds are split into interactions between the TCR-MHC (top) and TCR-Peptide (bottom). For TCR-MHC interactions, donor-acceptor pairs with greater than 5% existence occupancy in at least 1/9 simulations are included. For TCR-Peptide interactions, donor-acceptor pairs with greater than 5% existence occupancy in at least 1/3 simulations (for each respective peptide) are included.

Fig. 7

Average Lennard-Jones contact heat maps as function of TCR-pMHC COM distance. The time axis of the LJ-Contact existence maps is converted to COM distance by distributing time points into ∼0.5 Å bins and calculating the fractional occupancy in each respective bin. The 100, 95, and 90 ns configuration maps of each peptide are arithmetically averaged. The fractional occupancy is represented by the heat scale on the y-axis (right). Peptides (A) L1, (B) MART1, and (C) GVA are comprised of three SMD simulations. The bond acceptor and donor index is specified on the y-axis (left). The index corresponds to donor-acceptor pairs and is in Fig. S8. In each panel, LJ-Contacts are split into interactions between the TCR-MHC (top) and TCR-Peptide (bottom). For TCR-MHC interactions, donor-acceptor pairs with greater than 80% existence occupancy in at least 1/9 simulations are included. For TCR-Peptide interactions, donor-acceptor pairs with greater than 80% existence occupancy in at least 1/3 simulations (for each respective peptide) are included.

Finally, intra-mutant bond lifetime heterogeneity may depend on the existence of interactions from the starting configuration. This can be uniquely interrogated by examining inter-simulation differences (i.e., differences between the 90, 95, and 100 ns equilibration times) in the existence maps for a single peptide (Figs. S9 and S10). For example, the donated H-Bond of TCR⍺ Lysine 66 in L1 is present in simulations with ∼40 ns bond lifetime yet absent in the simulation with ∼4 ns bond lifetime (Fig. S9). Likewise, for MART1, the donated hydrogen bond of TCRβ Lysine 57 is present in the simulations with ∼4 ns bond lifetime yet absent in the simulation with ∼1.5 ns bond lifetime (Fig. S9). This analysis provides insight that underscores the importance of performing independent simulations that capture the ensemble of equilibrium configurations.

Conclusion

Our results demonstrate that TCR-pMHC bond strength corresponds with immunogenicity. Bond strength is characterized by the formation and stabilization of transient bonds/contacts as the TCR and pMHC are separated. Interestingly, these bonds/contacts are determined by peptide-mediated conformational shifts in the MHC binding groove. These conformational changes at equilibrium result in unique TCR-MHC interaction signatures during separation. The localized energetics indicate that, for this TCR-pMHC system, H-Bonds are crucial for TCR-MHC interactions whereas LJ-Contacts are more consequential for TCR-Peptide interactions. Although the generalizability of these results is limited to the TCR-pMHC interaction presented, both H-Bonds and LJ-Contacts contribute to overall bond strength that corresponds with the immunogenic response. As a result, TCR mutagenesis strategies – aimed at targeting a predetermined pMHC – may make informed mutations to affect the localized energetics and analogously evaluate the relative strength of interaction in silico. The atomic level detail of the SMD simulations also provide the opportunity to identify recognition hot spots (e.g., CDR1⍺, CDR2⍺) and specific amino acids in the TCR and MHC (e.g., TCR⍺ Arginine 27, TCR⍺ Tyrosine 50, TCR⍺ Lysine 66) that significantly determine overall bond strength. This level of insight may inform computational strategies to enhance TCR immunogenicity towards novel peptide targets by making targeted mutations that affect the localized energetics at the TCR-pMHC interface. Moreover, SMD may be used as a potential pre-screening approach to overcome the daunting experimental task of matching TCRs to pMHC targets from sequenced repertoires [53], [54]. Continued advancements in TCR design principles will provide a discovery platform with broad application in immune-oncology and infectious disease.

Methods

Molecular dynamics setup

The crystal structure of the human DMF5 TCR complexed with agonist peptide MART1-HLA-A2 (PDB code: 3QDJ) was the initial structure for all simulations [46]. Amino acid substitutions were made to the MART1 peptide (AAGIGILTV) using the Mutagenesis plugin on Pymol Molecular Graphics System (Schrodinger, LLC; New York, NY) to generate starting configurations for the L1 (LAGIGILTV) and GVA (GAGIGVLTA) peptides. Interfacial substructures were defined by sequential residues from the corresponding chains: TCR⍺ (CDR1⍺: 24–32, CDR2⍺: 50–55, CDR2⍺: 89–99), TCRβ (CDR1β: 25–31, CDR2β: 51–58, CDR3β: 92–103), MHC⍺ (MHCβ: 50–85, MHC⍺: 138–179), and peptide (L1: 1–9, MART1: 1–9, GVA: 1–9). To determine protonation states, pKa values were calculated using propka3.1 [55], [56] and residues were considered deprotonated in Gromacs [57] if pKa values were below physiologic pH 7.4. The resulting systems were solvated in rectangular water boxes using the TIP3P water model [58] large enough to satisfy the minimum image convention in order to avoid artifacts by self-interaction, a standard and well-described technique [59]. Na+ and Cl- ions were added to neutralize protein charge and reach physiologic salt concentration ∼150 mM. All simulations were performed with Gromacs 2019.1 [57] using the CHARM22 plus CMAP force field for proteins [60] and orthorhombic periodic boundary conditions. All simulations were in full atomistic detail.

Energy minimization and equilibration

Generating equilibrated starting structures for the steered molecular dynamics simulations required four steps: (1) steepest descent energy minimization to ensure correct geometry and the absence of steric clashes; (2) 100 ps simulation in the constant volume (NVT) ensemble to bring atoms to correct kinetic energies. A temperature of 310 K was maintained by coupling all protein and nonprotein atoms in separate baths using the velocity rescaled thermostat with a 0.1 ps time constant [61]; (3) 100 ps simulation in the constant pressure (NPT) ensemble using Berendsen pressure coupling [62] and a 2.0 ps time constant to maintain isotropic pressure at 1.0 bar; and (4) production MD simulations were conducted for 100 ns with no restraints. Steps (2)-(3) used position restraints on all protein atoms. To ensure true NPT ensemble sampling during 100 ns production runs, the Nose-Hoover thermostat [63] and Parrinello-Rahman barostat [64] were used to maintain temperature and pressure with time constants 2.0 and 1.0 ps, respectively, utilizing the isothermal compressibility of water, 4.5-5 bar−1. Box size for equilibration was 10.627 × 7.973 x.10.685 nm3 with ∼48,000 water molecules, ∼300 ions, and ∼157,000 total atoms (Table 1). All simulations used the Particle Ewald Mesh algorithm [65], [66] for long-range electrostatic calculations with cubic interpolation and 0.12 nm maximum grid spacing. Short-range nonbonded interactions were cut off at 1.2 nm. All water bond lengths were constrained with SETTLE [67] and all other bond lengths were constrained using the LINCS algorithm [68]. The leap-frog algorithm was used for integrating equations of motion with a 2 fs time steps. After 100, 95, and 90 ns of production run, respectively, MD configurations for each peptide mutant were extracted and used as the three initial configurations for steered molecular dynamics simulations.

Table 1

Equilibration of mutant peptides.

TCR-pMHC	Peptide Amino Acid Sequence	Equilibration Water Molecules	Equilibration Ions	Equilibration Total Atoms
L1	LAGIGILTV	47998	164 Na, 144 Cl	157134
MART1	AAGIGILTV	48001	164 Na, 144 Cl	157134
GVA	GAGIGVLTA	47997	164 Na, 144 Cl	157110

Steered molecular dynamics (SMD)

The full protein structure was extracted for each peptide mutant (L1, MART1, GVA) at the above-described production MD configurations (100, 95, 90 ns). Following solvation and the addition of Na + and Cl- ions, as described earlier, all systems underwent (1) energy minimization (2) 100 ps NVT (3) and 100 ps NPT. Details of the SMD are provided in Table 2. During the separation (“pull”) of the TCR from the pMHC, the Nose-Hoover thermostat and Parrinello-Rahman barostat were again used to maintain temperature and pressure. In SMD, there are two ways to apply load: controlling pull force or pull rate. Since this level of detail is not known in vivo, neither approach is fundamentally superior. Both methods have been used in the literature [69] and the appropriate selection may depend on the experimental technique used to measure force-dependent kinetics (e.g., laminar flow chamber: constant force and biomembrane force probe: constant pull rate). Comparative effects of load on TCR- pMHC force-dependent dissociation kinetics cannot be evaluated when controlling the pull rate because the simulation timescale is fixed. Therefore, we selected to control pulling force because the focus of this study is, in part, comparing the force-dependent dissociation kinetics of TCR-pMHCs with known biological function. Previously, using a constant pulling rate of 2 nm/ns, a critical force of 400–500 kcal/mol/nm (2800–3500 pN) was attained to separate the TCR and pMHC [36]. Moreover, another previous study used a constant pull rate of 0.1 nm/ns with a spring constant of 70 pN/nm to separate the TCR and pMHC, however, the critical applied force was not measured [37]. Both these examples represent the current state-of-the-art in SMD, but both apply forces that exceed those present in vivo due to constraints on computational speed. Comparable to the magnitude of applied force in previous studies, we chose to apply a constant force of 500 pN to the center of mass (COM) of the TCR and pMHC in the x-direction, until the distance between the COMs reached 0.49 times the box size. This was more than adequate to fully separate the TCR and pMHC (i.e., interaction energy between the TCR and pMHC reached zero) and this force is similar in magnitude to previous studies. Nonetheless, while observing bond lifetimes at a force of 500 pN is computationally feasible, we understand that this is higher than the physiological force range. Thus, we demonstrated, for a series of configurations along the reaction coordinate, that the root mean square fluctuations of the TCR-pMHC interface are not different at a physiological separation force of 10 pN (Fig. S2). These force control simulations (i.e., 10 pN) were pulled from the COM and performed for 50 ns on multiple configurations along the reaction coordinate (i.e., 6.16 and 6.20 nm) of the 100 ns MART1 configuration 500 pN pull.

Table 2

Steered MD simulations of mutant peptides.

TCR-pMHC	Peptide Amino Acid Sequence	Equilibration Time (ns)	SMD Water Molecules	SMD Ions	SMD Total Atoms
L1	LAGIGILTV	100	118977	363 Na, 343 Cl	370469
		95	118993	363 Na, 343 Cl	370517
		90	118960	363 Na, 343 Cl	370418
MART1	AAGIGILTV	100	118968	363 Na, 343 Cl	370433
		95	118925	363 Na, 343 Cl	370304
		90	118934	363 Na, 343 Cl	370331
GVA	GAGIGVLTA	100	118977	363 Na, 343 Cl	370448
		95	119013	363 Na, 343 Cl	370556
		90	118973	363 Na, 343 Cl	370436

Although the TCR is embedded in the cell membrane from the termini, the TCR is also in complex with the CD3 coreceptor [70] and thus it remains a significant challenge to determine the precise effective location of applied load in vivo. The COM was chosen as the site of applied force because pulling from the TCR and MHC termini results in artificial unfolding (Fig. S1). The termini control simulation was performed on the 100 ns MART1 configuration, with equivalent force magnitude (i.e., 500 pN), and was pulled from the COM of the terminal chain residues (MHC⍺: 275, TCR⍺: 199 & TCRβ: 242) (Fig. S1). All simulation trajectories and selected frames visualized using the Pymol Molecular Graphics System (Schrodinger, LLC; New York, NY).

Endpoints and data analysis

Data analysis was performed by built in Gromacs functions, standard python packages for data handling and visualization (i.e., numpy, pandas, seaborn, matplotlib, statistics, and GromacsWrapper), and custom python scripts. The geometry of a LJ-Contact is defined as a distance<0.35 nm between atoms. The principal component analysis of the simulation trajectories was performed using the package MDAnalysis [71], [72]. Secondary structure was designated based on a hydrogen bond estimation algorithm [73], [74]. Custom scripts relevant to the production of figures have been made available on a Github repository: https://github.com/zrollins/TCR-pMHC.git.

Statistical analysis

All data was processed in python utilizing standard packages for data handling and visualization (i.e. numpy, pandas, seaborn, matplotlib, statistics, scipy, and pingouin). Results were presented as mean ± SEM. As indicated in figures, statistics were performed in python using scipy for Student’s t-tests, scipy for one-way analysis of variance (ANOVA), and pingouin for pairwise Tukey-HSD post-hoc tests. Detailed outputs of statistical analysis were written to excel and are provided as supporting information.

Data availability

Starting configurations and Supplementary Videos have been made available: https://doi.org/10.25338/B8FK8D. Custom scripts relevant to the production of figures have been made available on a Github repository: https://github.com/zrollins/TCR-pMHC.git.

Author contributions

ZAR performed the simulations, analyzed and interpreted the data, and wrote the manuscript. RF designed the experiments, analyzed and interpreted the data, wrote the manuscript, and secured computer time. SCG designed the experiments, analyzed and interpreted the data, wrote the manuscript, and secured the funding.