Asymmetric framework motion of TCRαβ controls load-dependent peptide discrimination

Mechanical force is critical for the interaction between an αβT cell receptor (TCR) and a peptide-bound major histocompatibility complex (pMHC) molecule to initiate productive T-cell activation. However, the underlying mechanism remains unclear. We use all-atom molecular dynamics simulations to examine the A6 TCR bound to HLA-A*02:01 presenting agonist or antagonist peptides under different extensions to simulate the effects of applied load on the complex, elucidating their divergent biological responses. We found that TCR α and β chains move asymmetrically, which impacts the interface with pMHC, in particular the peptide-sensing CDR3 loops. For the wild-type agonist, the complex stabilizes in a load-dependent manner while antagonists destabilize it. Simulations of the Cβ FG-loop deletion, which reduces the catch bond response, and simulations with in silico mutant peptides further support the observed behaviors. The present results highlight the combined role of interdomain motion, fluctuating forces, and interfacial contacts in determining the mechanical response and fine peptide discrimination by a TCR, thereby resolving the conundrum of nearly identical crystal structures of TCRαβ-pMHC agonist and antagonist complexes.


30
The TCR ( TCR) consists of the heterodimeric receptor TCR formed by and chains each 31 containing the pMHC-binding variable (V) and constant (C) domains (Figure 1A), and the noncova- iments. We refer to a simulation as either low or high load based on the average load, which was 83 around the physiological 10-20 pN range (Table 1). To our knowledge, the present study is the first 84 to elucidate the dynamic mechanism of the A6 complex harboring different peptides under load. 85 We found that differences between the WT and peptide mutants lie in dynamic responses to 86 applied load. In the WT, physiological level load stabilized the TCR -pMHC interface as well as present results suggest that the conserved TCR framework motion is leveraged when determin- 94 ing the mechanically matched pMHC, a mechanism that is broadly applicable to different TCR 95 systems.

97
We first study the WT-based systems to gain insight into the functional implications of TCR -pMHC . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint Load stabilizes WT TCR -pMHC interfacial contacts 102 We assessed the effect of load on WT first by counting high-occupancy contacts with pMHC ( Fig-103 ure 2A; Contact analysis). WT low had the least number of contacts, followed by WT 0 and WT high , indi-104 cating low and high loads may have opposite effect on the interfacial stability. V -pMHC without 105 the C-module formed the most contacts. This indicates that without a proper load, the C-module 106 is detrimental to the stability of the interface with pMHC, as we found previously for the JM22 TCR 107 (Hwang et al., 2020). 108 Time-dependent changes in the interfacial contacts were monitored using the Hamming dis-109 tance  (Hamming, 1950; Hwang et al., 2020).  is the number of the initial high-occupancy con-110 tacts (those with greater than 80% average occupancy during the first 50 ns) that are subsequently 111 lost during the simulation. A low  means that such contacts persist while a high  means the 112 corresponding number of initial high-occupancy contacts are missing. Consistent with the contact 113 count,  remained low for WT high and V -pMHC ( Figure 2B). In WT 0 and WT low ,  increased af-114 ter about 500 ns and 900 ns, respectively. Thus, the relatively high number of interfacial contacts 115 for WT 0 (Figure 2A) is due to the formation of new contacts rather than by maintaining the initial  Figure 2D,E). Differences in the interfacial contacts also manifest 121 in their location. We displayed the backbone C atoms of V and V residues that form contacts 122 with pMHC with greater than 80% average occupancy (averaging was done after the initial 500 ns; 123 Figure 2F). In WT 0 , contacts are spread apart, and in WT low they lie mostly along the length of the 124 peptide. These layouts potentially make interfacial contacts more prone to break via easier access 125 by water molecules. In WT high and V -pMHC, high-occupancy contacts form more compact clus-126 ters. Exposure to water of the TCR residues involved in those contacts was measured by their 127 buried surface area (BSA), which follows the same trend as the number of high occupancy contacts 128 (Figure 2A vs. G). Furthermore, this trend also applied to the root-mean square fluctuation (RMSF) 129 of Tax peptide backbone C atoms, even though RMSF values were small ( Figure 2H). 130 Experimentally, the WT complex has a relatively strong affinity as a TCR (about 1 M) (Ding 131 et al., 1999), which may be why the interface with pMHC involved more contacts in WT 0 compared 132 to WT low . At low load, the short distance between restraints on the ends of the complex ( Figure 1A) 133 allows wider transverse motion that in turn generates a shear stress or a bending moment at the 134 interface. Such a transverse stress will be less for WT 0 where the end moves freely, and for WT high 135 where lateral motion is suppressed. The high stability of V -pMHC and WT high agree well with the 136 results for JM22 . 137 138 We analyzed the motion between V and V (V -V motion) to find its effect on the TCR -pMHC in-139 terface. Compared to the unliganded systems (V and T ; Table 2), the number of high-occupancy 140 V -V contacts increased slightly in V -pMHC ('V' in Figure 3A), while it decreased in full TCR -141 pMHC complexes ('0', 'Low', and 'High' in Figure 3A). This shows that the V -V interface is difficult 142 to organize with the restrictions imposed by the bound pMHC, except in the absence of the con-143 stant domains. The number of V -V contacts in the liganded systems is the smallest for WT low , 144 similar to the case for the number of contacts with pMHC (Figure 2A), again reflecting a destabiliz-145 ing effect with low load. . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Influence of pMHC and load on V -V motion
The copyright holder for this preprint this version posted September 13, 2023.   directions was measured by the absolute value of the dot product between PCs as 18-dimensional 152 unit vectors (for the 6 arms from two triads) in different systems. A value of 1 corresponds to the 153 same direction, and 0 means an orthogonal direction (Figure 3-figure Supplement 1C). For PC1, a 154 high degree of similarity was observed between T and V , which is consistent with their similar-155 ity in the number of V -V contacts ( Figure 3A) and PC amplitudes (Figure 3-figure Supplement 1A). 156 Among triad systems with bound pMHC, WT low differed significantly in the PC1 direction compared 157 to others (Figure 3-figure Supplement 1C, darker colors). The dot products varied more for PC2 158 and PC3, which capture finer motions with smaller amplitudes (Figure 3-figure Supplement 1C). 159 To determine how the V -V motion influences the interface with pMHC, we measured the dis-

169
PCA decomposes the V -V motion into mutually orthogonal directions. We made 2-dimensional 170 histograms of each of these projections versus the corresponding CDR3 distance (Figure 3-figure   171 Supplement 1D). If any of the PC modes is strongly correlated with the CDR3 distance, the corre-172 sponding histogram would exhibit a slanted profile. However, no clear correlation could be seen 173 (Figure 3-figure Supplement 1D), suggesting that the changes in CDR3 distance may depend on 174 combinations of PCs. We addressed this possibility by considering angles between matching arms 175 of the two triads ( Figure 3B). The 1 -1 angle, herein called ∠ 1 (and similarly define ∠ 2 and ∠ 3 ; 176 Figure 3B), can change either by the 1 arms turning up and down ("flap") or in and out of the page 177 ("twist") in Figure 3B. Angles ∠ 2 and ∠ 3 depend primarily on rotation indicated by dashed arrows 178 in Figure 3B ("scissor") (Hwang et al., 2020). 179 Histograms of the three angles ( Figure 3C) show a clearer difference than individual PCs among 180 the systems tested, and the CDR3 distance varies with the angles (Figure 3D). The wider distribu-181 tions for angles of WT low and WT high reflect their higher PC amplitudes (Figure 3-figure Supple-182 ment 1A). The symmetric distributions of ∠ 2 and ∠ 3 can be seen from the two peaks for WT low , 183 which is due to the reciprocal behavior of the scissoring motion involving the two angles ( Figure 3B, 184 side view, and Figure 3C, open diamonds). The two peaks are also related to the changes in the 185 CDR3 distance (Figure 3D), which reflects an agitating effect of the mild load on the scissor motion.

186
Given the definitions of the angles, the CDR3 distance will increase (dashed arrows in Figure 3B) 187 with larger ∠ 1 or ∠ 3 , or with smaller ∠ 2 ( Figure 3D). WT 0 , despite the apparent stability of the 188 interface, had a larger CDR3 distance than WT high and V -pMHC, again indicating a disrupted in-189 terface. These results show that the CDR3 motion is coupled to the V -V motion, especially the 190 scissoring motion.

191
Asymmetric V-C motion influences the load response of the complex 192 We next analyzed the motion of the V-module relative to the C-module (V-C motion). The number 193 of high-occupancy C -C contacts did not vary significantly (in the 31-34 range) and they were 194 more than the number of V -V contacts, similar to the case for the JM22 TCR (Hwang et al., 2020). 195 The C-module thereby influences the V-module as a single unit. The V-C motion was analyzed by

202
We noticed that V bends more compared to V , as can be gleaned from the longer PC ar-  Figure 4C). Compared to T , binding of pMHC increases the -chain motion, which is the greatest 207 in WT low (Figure 4C, PC1 in top row). The greater degree of V -C motion is consistent with the 208 smaller number of V -C contacts compared to V -C (Figure 4-figure Supplement 1B,C). 209 The asymmetry was further analyzed by measuring hinge angles ∠TCR and ∠TCR ( Figure 4A). 210 Distributions of ∠TCR varied more compared to ∠TCR (Figure 4D). A wide distribution of ∠TCR 211 for WT low is related to the increase in the CDR3 distance and concomitant changes in the triad arm 212 angles later during the simulation (Figure 3-figure Supplement 2C). For WT low and WT high , CDR3 213 distance decreases with increasing hinge angles, especially with ∠TCR ( Figure 4E), which suggests 214 that unbending of the V-module under load helps with bringing the CDR3 loops closer together.

215
In WT low , this state is not maintained and ∠TCR decreases (more bending) as the CDR3 distance 216 increases ( Figure 4E) which happens after the increase in  (Figure 2B). These results suggest a 217 7 of 27 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint mechanism by which the asymmetric response of the whole TCR to load affects the binding with

Effects of point mutations on the peptide 222
In the WT crystal structure, the side chain of Y5 in the Tax peptide is located between the CDR3 loops 223 of V and V while V7 mainly contacts CDR3 ( Figure 1C). P6 makes one contact with CDR3 . The . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint    ber of contacts with pMHC for modified agonists ( Figure 5A) despite an increase in  suggests an 238 altered binding rather than maintaining the initial contacts.

239
In contact heat maps, the Y5 residue of the WT peptide forms a hydrogen bond with S31 and 240 nonpolar contacts with a few residues in both V and V ( Figure 2C-E, Figure 5C). In Y5F, the hy- . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint     positions; its absence would make the nonpolar contacts more prone to break under load. The 247 relative stability of Y5F 0 can also be seen by the similarity in the locations of high-occupancy con-248 tact residues between WT and Y5F 0 ( Figure 2F vs.  The modified agonists had more V -V contacts than WT high while the antagonists had fewer, 263 except for Y8A high (Figure 6-figure Supplement 1A). While the amplitude of V -V motion was  Supplement 1B), the CDR3 distance was larger for all mutant systems except for P6A, which had 266 a weak dependence on triad angles (Figure 6A-B). The angles in turn varied among systems and 267 loading conditions (Figure 6-figure Supplement 2). These results suggest that point mutations to 268 the peptide cause alterations in the load-dependence of the interface and the V -V motion.

269
The mutants affected the average V-C BOC similarly as that for dFG high (Figure 6C,D vs. Appendix 270 1- Figure 1C). Among them, Y8A high had an average BOC approaching that of WT high , which aligns 271 10 of 27 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint with the comparable number of contacts with pMHC ( Figure 5A). However, the location of its H 272 differed (Figure 6D), and the CDR3 distance was larger ( Figure 6B). To quantify deformation of the 273 average BOC, we measured displacements of centroids from the corresponding ones in WT high .

274
They were overall greater for the chain than the chain (Figure 6-figure Supplement 3A,B). 275 Consistent with this, the mutants had fewer V -C contacts than WT high and a similar number of 276 V -C contacts (Figure 6-figure Supplement 3C,D). 277 Similar to the WT systems, the greater motion of the chain than the chain was observed in 278 the mutant systems, as seen from the differences in V-C PC1 amplitudes (Figure 6- Figure 4B). Thus, point mutations on the WT peptide can affect the con-  To probe the dynamic relation between the TCR -pMHC (intermolecular) interface and intra-TCR 291 (intramolecular) interfaces formed between subdomains of the complex, we calculated the total 292 occupancy of the high-occupancy contacts in respective cases (Figure 7A,B). For the intramolecular 293 contacts, we excluded the C -C interface contacts since they are larger in number compared to 294 other interfaces and did not differ significantly across different systems, i.e., the C-module moves 295 mostly as a single unit (Hwang et al., 2020). 296 For WT 0 , the intermolecular contact occupancy stayed at around 20 (WT in Figure 7A, horizontal 297 bar on the bottom) and for WT low , it decreased later in simulation (WT in Figure 7A, darkening of 298 circles without outline). In comparison, the intramolecular contact occupancy remained relatively 299 constant for both WT 0 and WT low (WT in Figure 7B, horizontal bar on the bottom and circles without 300 outline). For WT high , the intermolecular contact occupancy was steady even with wider fluctuation 301 in force (WT in Figure 7A, outlined circles), and the intramolecular occupancy also remained high, 302 indicating the subdomains are held together tightly (WT in Figure 7B, outlined circles). For dFG low , 303 the intermolecular contact occupancy stayed low and intramolecular occupancy was relatively con-304 stant (dFG in Figure 7A,B, circles without outline). In dFG high , the contact occupancy with pMHC 305 increased (dFG in Figure 7A, outlined circles), but the intramolecular contact occupancy became 306 low (dFG in Figure 7B, outlined circles), which suggests that the complex is not as tightly coupled 307 compared to WT.

308
For modified agonists, the no load and low load cases had overall higher occupancy, both with 309 pMHC and within TCR , but occupancy fluctuated more as can be seen by the changes in colors in 310 the occupancy trajectories (Y5F and V7R in Figure 7A,B, horizontal bars on the bottom and circles 311 without outlines). Under high load, intermolecular contact occupancy decreased over time (Y5F 312 and V7R in Figure 7A, darkening of outlined circles) while intramolecular contact occupancy either 313 increased (Y5F high ) or decreased (V7R high ) relative to the respective low load cases. For antagonists, 314 both occupancy measures were lower than the WT, and further reduction could be seen over time 315 in some cases (P6A and Y8A in Figure 7A,B, darkening of colors in outlined circles).

316
The stability of the TCR -pMHC interface also manifested into their relative motion, which   dFG low , the peptide changed orientation by more than 20 • , and for dFG high , it stabilized, but at 322 a higher value than WT high , which also was reached in dFG low later during simulation, suggesting 323 a more orthogonal binding (dFG in Figure 7C). For modified agonists, similar to the behaviors of 324 the total intra-and intermolecular contact occupancy, the peptide angle was affected more under 325 high loads, again becoming more orthogonal compared to WT high (Y5F and V7R in Figure 7C). For

332
The present study elucidates how the load-dependent TCR framework motion influences the 333 dynamics of the TCR -pMHC interface (Figure 8). A main feature of the framework is the smaller 334 number of contacts for the V -C compared to the V -C interface. This causes an asymmetric V-C 335 motion, primarily bending, where V moves more compared to V relative to the C-module, which 336 serves as a base. This in turn generates relative motion between V and V , which can destabilize 337 the contacts with pMHC, especially by affecting the distance between CDR3 loops that play the 338 most direct role for sensing the bound peptide (Figure 8A,B). Applying a physiological level force 339 stabilizes the interface by straining the whole complex into a more tightly coupled state, as can be 340 seen by the increase of both inter-and intramolecular contacts in WT high (Figure 7 and Figure 8C). 341 The CDR3 distance of WT high (10.3±0.3 Å; Figure 3-figure Supplement 2C) was shorter than 342 that of WT 0 or WT low (Figure 3D), and it is also shorter than the 10.9-Å CDR3 distance in the crystal 343 structure (PDB 1AO7). The applied load slightly increases the spacing between pMHC and TCR , 344 which provides room for the CDR3 loops to adjust as well as allow other contacts to 'lock' into 345 more stable states with higher and more persistent occupancy. Absence of load or low load do not 346 12 of 27 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint properly channel the framework motion and thereby increase exposure to water (Figure 2F,G) and 347 destabilize the interface.

348
The C FG-loop stabilizes the V -C interface, thereby contributing to the asymmetric V-C mo-349 tion. It also controls the relative orientation between V and C , hence it affects the orientation 350 of the CDR loops of the V-module with respect to the loading direction (Appendix 1- Figure 1C). 351 Consistency in these findings between the present study and our previous simulations using JM22 352 TCR (Hwang et al., 2020) underscores that the proposed mechanism based on the asymmetric 353 framework motion is applicable to other TCR systems.

354
After engagement with a cognate pMHC under load, reversible transition to an extended state 355 is possible which has been observed both in in vitro single-molecule experiment using TCR and 356 on cell displaying the full TCR holoreceptor (Das et al., 2015; Banik et al., 2021). Since V -pMHC 357 lacking the C-module forms a more stable binding (Figure 2) that was also observed in our previous 358 simulations of the JM22 TCR (Hwang et al., 2020), the C-module likely undergoes partial unfolding 359 in the extended state. Thus, while the folded C-module serves as the base for the asymmetric  plex. Yet, for T-cell based cancer immunotherapy, mechanistic knowledge of the mechanosensing 372 through a TCR has a greater practical significance (Reinherz et al., 2023). 373 A recent study using a laminar flow chamber assay fit the measured bead survival distribution 374 using Bell's equation to estimate the zero-force off rate off and the force sensitivity distance 375 (Pettmann et al., 2022). They found a negative correlation between off and , to conclude that me-376 chanical forces impair antigen discrimination. However, the force range tested was up to 100 pN, 377 where even systems exhibiting catch bond in the 10-20-pN range will switch to a slip bond behavior. . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint than the length of a single covalent bond. They also performed steered MD simulation that applies 382 hundreds of pN forces, which is inadequate for studying behaviors of the system under loads in 383 the 10-20-pN range (Hwang et al., 2020). Use of a coarse grained model without appropriately 384 incorporating atomistic properties of the TCR further makes it difficult to compare their simulation 385 with experiment. 386 We earlier proposed that the residues of the antigenic peptide play a role more as "teeth of a 387 key" for screening the TCR -pMHC interaction fitness rather than bearing applied loads (Hwang   388   et al., 2020; Reinherz et al., 2023). The present study confirms this through simulations of mutant 389 systems, where several contacts across the interface with pMHC were impaired due to a single-390 residue mutation on the peptide in ways that reflect the functional outcome of the mutation. In 391 considering how a T-cell may respond to an unknown peptide, the pMHC motion and the asymmet-392 ric V-C motion are two points of guidance (Figure 8). Stabilization of the inter-and intramolecular 393 interfaces throughout the whole complex under 10-20-pN load would indicate a cognate TCR -394 pMHC interaction (Figure 7). Since these features are based on overall TCR -pMHC complex dy-395 namics, rather than changes to specific contacts or a particular conformational change, they can 396 be used to predict fitness of other TCR -pMHC combinations. Since such tests involve performing 397 many all-atom MD simulations and trajectory analyses, an in silico method would be needed that 398 efficiently predicts dynamic properties of the complex based on sequence and structural data only.

399
Atomistic insights gained from the present study will be helpful for developing such a method in 400 future studies. were introduced as noted in the PDB file. Histidine protonation sites were determined based on 411 the 1QSE crystal structure to promote hydrogen bond formation with neighboring residues. Where 412 neighboring residues were unlikely to hydrogen bond, we assigned the water-facing nitrogen of his-413 tidine as charged. This led to protonation of the N atom for all histidine residues except for MHC 414 H263 and 2m H84, where the N atom was protonated. For truncated structures, crystal waters 415 within 2.8 Å from the protein atoms were kept in the initially built system. For full structures, all 416 crystal waters were kept. 417 We extended the termini of the TCR -pMHC complex as handles for applying positional re-418 straints (Figure 1, "added strands") (Hwang et al., 2020). For MHC, we used the sequence from (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint C domain residue coordinates to 1AO7. After this, we performed a brief energy minimization on 431 the added domain while fixing positions of all other atoms of 1AO7. For missing residues in the C 432 domain, we used backbone C atoms of two residues each before and after the missing part to align 433 1QSE to 1AO7 and filled in coordinates, followed by a brief energy minimization of the added part 434 in 1AO7. In this way, the TCR -pMHC interface of the original 1AO7 is preserved. By comparison, 435 previous simulations mutated PDB 1QRN back to WT (Ayres et al., 2016), which corresponds to 436 to the A6P in silico WT system (Appendix 2), or converted a high-affinity variant of A6 (PDB 4FTV) 437 by mutating -chain residues, in particular nearly the entire CDR3 loop (Rangarajan et al., 2018). 438 Compared to our approach, those preparation methods thereby introduce more perturbation to 439 the interface with pMHC.

440
The 2m residues C67 and C91 were reverted (C67Y, C91K) based on UniProt P61769 referenced 441 in PDB 1AO7. These agree with the 2m sequence in other structures.

3QFJ (Y5F):
There were no missing residues. We made the D204N conversion in TCR .

448
WT truncated complexes: For truncation, we used the constructed 1AO7 complex.

449
• V : the last residues were D111 and E116.

452
• WT 0 : WT complex without the added C-terminal strands, as for T .

453
• dFG: residues L218-P231 removed from the corresponding WT complex. G217 and V232 454 were covalently joined.  Table 1. For solvation, we used the TIP3P water. Water molecules with 463 their oxygen atoms less than 2.8 Å from protein heavy atoms were removed. Neutralization of the 464 system was done using Na + and Cl − ions at about 50 mM concentration. Crystal water molecules 465 were kept in this procedure.

466
After neutralization, a 5-stage energy minimization was applied where protein backbone and 467 side chain heavy atoms were progressively relaxed (Hwang et al., 2020). This was followed by heat- Production runs were performed using OpenMM (Eastman et al., 2017). We used the CHARMM 478 param36 all-atom force field (Huang and MacKerell Jr, 2013) and the particle-mesh Ewald method 479 to calculate electrostatic interactions. We used an Ewald error tolerance of 10 −5 which is 1/50 of 480 the default value in OpenMM, for accuracy. The cutoff distance for nonbonded interactions was 481 12 Å, and the Nose-Hoover integrator of OpenMM at 300 K was used, with a 2-fs integration time 482 step. We ran OpenMM on GPUs with mixed floating point precision. Below are specific steps of 483 the MD protocol relevant to individual systems in Table 1 and Table 2. for the next 4-ns run. The process continued to yield 4-6 extensions.

495
After each extension run, we truncated the water box such that the length of the box was 496 12 Å larger than the maximum span of the complex on each side, and re-neutralized the system.

497
A representative water box size is 218×97×90 Å 3 for WT high , containing 187,250 atoms. Since the 498 system was already equilibrated from the previous run, we used a simpler energy minimization 499 scheme where backbone and side chain heavy atoms were restrained by 10-kcal/[mol⋅Å 2 ] and 5-500 kcal/[mol⋅Å 2 ] harmonic potentials, respectively, and 200 steps of steepest descent followed by 200 501 steps of adopted-basis Newton-Raphson energy minimization was performed. Heating, equilibra-502 tion, and the initial 2-ns dynamic runs with positional restraints were carried out as explained 503 above. We then carried out 60-100 ns production runs for each extension and selected two or 504 three extensions to continue for longer than 1000 ns.

505
Selecting extensions 506 We measured the average force on the complex during each 60-100-ns simulation, then selected 507 two extensions where the average force generated was representative of a "low" (around 10 pN) 508 and "high" (over 15 pN) load on the TCR. These values were based on the experimental 10-20-pN 509 catch bond activation force range (Das et al., 2015; Liu et al., 2016). 510 In some cases, in particular at low extensions, the flexible added strand either folded onto itself 511 or made contacts with the C-module of TCR , effectively shortening the span of the complex. Fac-512 tors such as this, together with differences in conformational behaviors of the complex, affected 513 the average force for a given extension. Thus we had to test and choose among different exten-514 sions for each system. We also ran 1-2 replicate simulations of comparable length (∼1 s) at given 515 extensions except for systems involving dFG and in silico mutants. However, even with nearly the 516 same extensions used, measured forces in replicate simulations varied. The final selection and 517 average forces are in Table 1. 518 Other systems . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint below about 0.5 Å in both WT low and WT high , which were P185-T187, L201-Y209, F241-V247, and 524 T258-H263.

525
V -pMHC 526 We applied a 0.01-kcal/[mol⋅Å 2 ] harmonic restraint to the backbone C atoms of the MHC 3 do-527 main (residues P185-L276) to prevent the complex from turning transversely in the orthorhombic 528 box. The restraints are 20 times weaker than those used for TCR -pMHC complexes mentioned 529 above. This was because V -pMHC is smaller in both size and aspect ratio. The FG-loop deletion was done after initially preparing (solvation and neutralization) the WT com-535 plex in the extended water box. After deletion, the system was re-neutralized. Subsequently, lad-536 dered extension, selecting extensions for high and low load cases, and longer production runs 537 were performed as explained above.

538
In silico mutants 539 Each in silico mutation (Appendix 2- Table 1) was performed for low and high load extensions of 540 the complex. To use similar extensions as in the original complexes, we used the last frame of 541 the 4-ns laddered extension simulation. After introducing the in silico mutation, we inspected the 542 structure to ensure there was no steric clash with neighboring residues or water molecules. We 543 performed a short energy minimization to relax the modified residue while keeping coordinates 544 of all other residues except for residues immediately before and after the mutated one on the 545 peptide. We then truncated the water box and re-neutralized the system, after which steps from 546 the initial energy minimization up to the final production run followed the same procedure as 547 explained above.

549
Coordinates were saved every 20 ps (0.02 ns) during production runs, resulting in 50,000 coordi-550 nate frames for 1000 ns. We excluded the initial 500 ns when calculating averages and standard 551 deviations in the number of contacts, CDR3 distance, BSA, PCA values, and angle data. Since all 552 systems were simulated for a minimum of 1 s, this leaves at least 25,000 frames. We report data 553 prior to 500 ns in trajectory plots and contact occupancy heat maps (e.g., Figure 2B-E). 554 Calculating force 555 Force on a restrained atom or the center of mass of the C-terminal atoms of the added strands 556 in TCR was calculated based on the deviation of its average position from the center of the har-557 monic potential, multiplied by the spring constant used (Hwang et al., 2020). Average force in Ta-558 ble 1 was computed from 500 ns to the end of the simulation. Instantaneous force in Figure 7A RMSF for backbone C atoms of a domain was calculated by aligning the C atoms to the structure 563 at the beginning of the production run. Coordinate frames after the initial 500 ns were used. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint Contact analysis 569 We used our previously developed method (Hwang et al., 2020). Briefly, H-bonds (including salt 570 bridges) were identified with the 2.4-Å donor-acceptor distance cutoff. Nonpolar contacts were 571 identified for atom pairs that are within 3.0 Å and both have partial charges less than 0.3 ( = 1.6 × 572 10 −19 C) in magnitude. The average occupancy was measured as the fraction of frames over which 573 a bond is present during the measurement period. Instantaneous occupancy was measured as a 574 40-frame (0.8-ns) rolling average. The average occupancy of a contact represents its abundance 575 during the simulation period while the instantaneous occupancy represents its temporal intensity.

576
For counting the number of contacts (e.g., Figure 2A, Figure 3A, Figure 4-figure Supplement 1B,C), 577 we used contacts with the average occupancy greater than 50% and at least an 80% maximum in-578 stantaneous occupancy after the initial 500 ns. Contact occupancy heat maps (e.g., Figure 2C-E) 579 report those with the overall average occupancy greater than 30%, and the maximum instanta-580 neous occupancy during the simulation greater than 80%.

581
The Hamming distance  (e.g. Figure 2B) was measured using contacts with greater than 80% 582 average occupancy during the first 50 ns.

584
For the BSA calculation (e.g. Figure 2G), we used residues in the V-module with the maximum 585 instantaneous contact occupancy with pMHC greater than 80%. We calculated the surface area for 586 the selected residue contacts and added them to get the total BSA. Per-residue BSA is the total BSA 587 divided by the number of residues forming the contacts in the given time interval. The reported 588 values (e.g. Figure 2G) are respective averages after 500 ns. domain, which is parallel to the -strands and points to the CDR3 loop ( Figure 3B). The 1 arm was 599 assigned by taking the direction from the center of masses of the selected atoms from the inner 600 to the outer -sheets of each variable domain and making it perpendicular to 3 . The 2 arm was 601 then determined as 2 = 3 × 1 .

602
PCA was performed on the trajectory of the two triads using a custom FORTRAN95 program 603 (Hwang et al., 2020). The PC amplitude (e.g. Figure 3- Supplement 1C), the absolute value of the dot product between them was calculated, which ranges 607 between 0 and 1. To project the V -V triad for a given frame to a PC direction (Figure 3-figure   608 Supplement 1D), the average triad calculated after the initial 500-ns was subtracted from the triad, 609 then a dot product was formed with the PC vector. The V-C BOC (Figure 4A) was assigned based on the method we developed previously (Hwang et al.,   612   2020). For beads representing the V-module, centroids of the two triads were used. For the C-613 module, the center of mass of backbone C atoms of the following residues in each domain were 614 used: for C , A118-R123, V132-D137, Y153-T158, S171-S176, and for C , T143-A148, L158-N163, 615 S192-V197, F209-Q214. We used N114 for H , and for H , the center of mass between D117 616 and L118 was used, which had large RMSF in WT high . (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint We aligned coordinate frames for all simulations to the first frame of WT high based on atoms 618 used to assign beads for the C-module. In this way, motion of the V-module relative to the C-module 619 can be analyzed. Also, by using a common reference structure (first frame of WT high ), average 620 BOCs can be compared, as in Figure 6C,D. PCA of the V-C BOC was performed using the 6 beads 621 representing the centroids and hinges. PCA for the V-module triads was done separately. Since 622 the reference of motion is the C-module, directions of PCs for the V-module triads indicate motion 623 of the V-module relative to the C-module (arrows on triad arms in Figure 4A), which complements 624 the direction of the V-module centroids obtained from PCA of the V-C BOC (arrows on centroids in 625 Figure 4A). 626 Time-dependent behavior 627 For the total occupancy in Figure 7A,B, we only considered contacts with greater than 50% overall 628 occupancy and over 80% maximum instantaneous occupancy during the entire simulation period.

629
In this way, changes in high-quality contacts under fluctuating force for a given trajectory can be  The altered conformation of dFG causes the interface with pMHC to tilt as observed in our previous study of JM22, which is detrimental to the stability of the complex (Hwang  et al., 2020). The increased motion of H may also deliver more agitation to the interface with pMHC. Thus, even though full dissociation with pMHC was not observed within the simulation time, the dFG-pMHC complex is likely to be less stable compared to the WT complex. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint 2- Figure 1B, top), while for in silico WT their range became narrower, to 11-12 Å which is similar to that for the modified agonist Y5F (Figure 6A vs. Appendix 2- Figure 1B, bottom row). The CDR3 distance also stabilized over time in all in silico WT, which was less so for in silico antagonists (Appendix 2- Figure 2D). However, similarly as the number of V -V contacts, there was no consistent load-dependence between the CDR3 distance and triad arm angles. On the other hand, there was a stronger load dependence in the average V-C BOC. The in silico antagonists that were built based on WT bent towards those of the corresponding antagonists, though the extent was not large (Figure 6D vs. Appendix 2- Figure 1C). The in silico WT in low load had average BOCs similar to those of the original antagonists whereas average BOCs of high-load in silico WT approached those of the actual WT (Appendix 2- Figure 1D). For Y8A WT, this happened even though the forces experienced at the two extensions were only marginally different (10.0 vs. 11.5 pN; Appendix 2- Table 1). The V -C and V -C contacts were respectively lower for in silico WT than in silico antagonists, suggesting effects of the in silico mutations of the peptide did not propagate sufficiently to the whole TCR during the simulation time (Appendix 2- Figure 2E). The lower number of V-C contacts in in silico WT would have made it easier to unbend under higher load or extension. The above results suggest that the in silico mutants behave like the target system to varying extents. This is likely because the rearranged interface between the V-module and pMHC of the base system cannot immediately be adjusted upon in silico mutation in loaded states.  (Figure 6C,D). . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint  . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint

877
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint (C,D) Antagonists. The same cutoff criteria were used to calculate initial contacts as in Figure 2B. Data after 500 ns were used for histograms on the right of each panel.

879
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint  (Figure 2F), contacts are overall unevenly distributed or dispersed. The same occupancy cutoffs as in Figure 2F were used for selecting residues. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint  Figure 3C is included in all panels (without markers) for comparison. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint

884
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint and peptide (blue) from the crystal structure overlaid with MHC 2 helix (green) and peptide (magenta) of the last rendered frame from panels A-D.

885
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted September 13, 2023. ; https://doi.org/10.1101/2023.09.10.557064 doi: bioRxiv preprint