The influence exerted by the TCR-pMHC affinity on the immune response was studied by performing computations from 12 phenotypic models with the specific aim of analyzing the relationship between affinity and response. To this end, the following phenotypic models were considered:

occupancy

Models (a)-(d) are reviewed in Lever et al. (4), model (e) is described in Lever et al. (4) and François et al. (38), model (f) in Dushek and van der Merwe (39), model (g) in Gálvez et al. (37), and model (l) in Lever et al. (23). In addition, and because the modular structure of phenotypic models allows to study the effect exerted on the response by adding assumptions and new parameters to a simpler model, we have formulated the new models (h)-(k) by combining hypothesis of models (c) and (d), (c) and (e), (g) and (c), and (g) and (d), respectively. Thus, for example, in developing model (h) (kpr with limited and sustained signaling), we have considered that assumptions for limited signaling [namely, that TCRs having reached the signaling competent state can only signal for a limited period because of movement into the immunological synapse or to become tagged for removal (4)] and those for sustained signaling [i.e., that signaling competent TCRs are able to maintain signaling for a prescribed period of time, even after pMHC unbinding (4)] could be involved together in determining the T-cell response. Schemes of these models and of the corresponding parameters are shown in Figures 1A, B.

All models have a common set of parameters. Also, there are parameters specific to each model. The common set of parameter values used for computation in this work are similar to those used in Lever et al. (4), Altan-Bonnet et al. (14), Gálvez et al. (37), and François et al. (38): number of TCRs, TT=2×104; kon=5×10-5s-1; kp=1s-1. There is no concentrations units: all concentrations in figures, tables and in rate constants are per cell. Thus, kon=5×10-5(molecule.s)-1 and k_off = 1/τ where τ is the dissociation time of the TCR-pMHC complex. Besides the common parameters, models (c)-(l) also include specific parameters that are defined in Figures 1A, B. Their numerical values are listed next and are similar to those used in the accompanying references:

model (c), kpr with limited signaling (4): ϕ = 0.09s⁻¹. This value of ϕ was used in all models including limited signaling [models (c), (h), (i) and (l)].

In turn, the number of steps leading to the productive signaling complex (this parameter does not apply for the occupancy model) was N = 10 for all models except in the kpr with induced rebinding model where N = 25 [due to the assumed values of the rebinding rate constants (39), induced rebinding has little effect on the response when N < 25 (37, 39) and thus, if plots in this work are computed with N = 10 the resulting R-values are very close to those obtained with the basic kpr model].

Finally, it is worth to note that recent works have placed an added emphasis on experimental measurements of 2D k_on and k_off reaction rate constants vs. the corresponding 3D values (3, 41–44). Thus, while 3D TCR-pMHC interactions have been widely studied by surface plasmon resonance (SPR) (3, 44), determination of the 2D binding parameters is much more challenging and more elaborate techniques, such as adhesion frequency and thermal fluctuation assays (41), single-molecule microscopy and fluorescence resonance energy transfer (FRET) between fluorescently tagged TCRs and their cognate pMHC ligands (42), and a laminar flow chamber to monitor at the single molecule level the 2D TCR-pMHC interactions (43), were used. For computation purposes we have to take into account that 2D on-rates are actually effective on-rates with different units to the corresponding 3D on-rates which precludes that a direct comparison between 2D and 3D on-rates can be performed. Conversely, 2D and 3D off-rates have the same units (s⁻¹), and a direct comparison is possible. Nevertheless, in Zarnitsyna and Zhu (3) we have that 2D on-rates span a broad 4-log range while 3D on-rates are compressed into a narrow range (see Figure 3B in that reference). The set of on-rates used for computation in this work spanned both narrow and broad log ranges (from 0 to 4) so that the full range of 2D and 3D on-rates values is included. As regards the off-rates, Figure 3C in Zarnitsyna and Zhu (3) shows that the 2D k_off-values span approximately between 2 and 10 s⁻¹ while this range for the corresponding 3D values is from 0.001 to 0.1s⁻¹. The values of k_off used for computation in this work are between 0.0001 and 10s⁻¹ so that the full range of 2D and 3D off-rates is also covered.

All models in Figures 1A, B were described by systems of ordinary differential equations (ODEs). Numerical solutions of the system of ODEs as a function of time, as well as all remaining calculations and plots were performed using Mathematica 11.2. The system of ODEs for the above models are given in Appendix and in Supplementary Material.

Responses were computed both in transient phase as a function of the activation time and in steady-state. For most models, analytical steady-state solutions can be derived, although some of them are quite complicated and must be solved numerically. For this reason, responses were always computed in transient phase from t = 0 until values of t sufficiently large. Under these conditions a steady-state for all models was reached, which was confirmed by testing that the responses obtained from transient phase solutions at t>>1 gave the same values than those computed from analytical (or numerical) steady-state solutions. This procedure has the advantage that, in addition to show the behavior of the response at longer times (steady-state solution), its characteristics at shorter times (transient phase solution) can be also studied.

Since a set of common parameters is used, not only qualitative but also quantitative responses from different models can be compared. This also allows us to reveal the influence exerted on the response by those parameters which are specific for each model. To this end, non-normalized responses are always provided instead of normalized values because for comparative purposes normalized responses could be misleading or even meaningless if the normalization factors are different. On the other hand, the modular structure of phenotypic models allows to show the effect exerted on the response by the new assumptions and parameters added to a simpler model. For example, we shall find that responses in models including sustained signaling are much larger, but the time required to reach the steady state is also much longer. For computational purposes it is interesting to note that solving the system of ODEs as a function of time always gives a single solution (which is the correct one) for all phenotypic models. However, for models including negative feedback [models (e) and (i)] the steady-state solutions are obtained by solving a polynomial equation with several solutions of which only one is a valid value. To determine the correct value the safest way is to compare the polynomial solutions with the transient phase values at t>>1. Finally, it is worthy of mention that in those models including limited signaling (see below), conclusions drawn from quantitative analysis of R-values at short times and under steady-state conditions (t>>1) can be different. In this regard, it has been suggested (45, 46) that time-scales in T-cell activation can be more relevant than responses obtained under steady-state where the R-values are time-independent (this should lead to questioning why then T-cell responses are normally measured and analyzed in steady-state although, more likely, is due to that experimental measurements and theoretical solutions in steady-state are much easier to deal with than those in transient phase).

This section provides a brief description of the phenotypic models (a)-(l) used in this work, their mathematical assumptions, and their levels of complexity, bearing in mind that models (h)-(k) are new models which have not been previously described in the literature. The reader interested in more details on the problems discussed in this section can consult the review articles (3, 4, 47, 48) and references (14, 49–52). Also, the following considerations should be taken into account: (a) in general, quantitative predictions from a given model rely on the premises on which the model is built. Thus, in the case of TCR-pMHC interaction if a model is aimed at describing the binding process, its quantitative predictions are necessarily limited to the formation of the TCR-pMHC complex. Hence, predictions from this type of models that go beyond the binding process, v.g. on the length and kinetics of the activation chain, negative feedback, rebinding, and limited and sustained signaling, among others, would be only assumptions since such predictions can not be quantified by computation; (b) any model aimed at describing quantitatively T-cell responses should be able to consider and quantify the following key features of this response: its impressive capacity of speed, sensitivity and discrimination that allows to detect foreign pMHCs at very low concentration among much more abundant self-pMHC ligands (37–39); (c) recent works have placed an added emphasis on experimental measurements of 2D k_on and k_off reaction rate constants, comparison with the corresponding 3D values, and suggestions of how these 2D rate constants influence the T-cell response (3, 41–44). Unfortunately, a complete and rigorous theoretical integration of 2D and 3D domains in cell signaling processes is at present an extremely problem to tackle because we have to develop spatio-temporal models involving partial differential equations, deal with the complex problem of “reduction of dimensionality” (3D → 2D → 3D), and consider mass transport processes that show trafficking dynamics of ligands, cell surface components, intracellular signaling molecules through the different domains and interfaces (2D and 3D) of the system, among others (52). In turn, this means that a much larger set of kinetic parameters (most of which are unknown) are involved. Because of these difficulties phenotypic models are developed that, although they do not consider in detail all signaling events and spatial domains, they can provide a reasonable overview of the overall process on the basis of a minimum set of assumptions and parameters. This has been the approach followed in the present manuscript.

Summarizing: the developing of these models show that to take into account a new effect we need to figure out how to incorporate this effect (and the kinetic steps in which it is involved) in the mathematical framework of the phenotypic model. Thus, the response, R, in each model refers to the assumptions and parameters on which the model is built. For example, the occupancy model only considers the first binding step, while in the basic kinetic proofreading model R takes into account both the first binding step and the activation proofreading chain.

TCR-pMHC affinity is given by A = k_on/k_off. Hence, the values of k_on and k_off should change proportionally so that affinity remains constant. In Figures 2A, B we have plotted the predicted T-cell responses given by R (see Appendix) computed from the 12 phenotypic models for three different systems with equal affinity and the following values of k_on and k_off:

where the value of k_on is given in subsection 2.2.

The responses computed from the occupancy model are shown in panel a of Figure 2A. As expected, the transient phase is shorter as the k_on-values increase, but at longer times, once the steady state is attained, the R-values become independent of the individual values of k_on and k_off. Also, and because the three systems have equal affinity the corresponding steady-state responses show the same value in agreement with Equation (A1) in Appendix. However, in the basic kpr model (Figure 2A, panel b) the responses are quite different despite the three systems have the same affinity. Thus, the largest response is obtained for the system with the lowest values of k_on and k_off (blue curve). On the other hand, the responses for the kpr with limited signaling model (Figure 2A, panel c) are more complicated. As a result, the largest responses appear as peak values in transient phase with system 1 showing the largest value of R. However, as signaling progresses, the R-values decrease so that in steady state the response for system 1 becomes the lowest of the three systems. In turn, panel d of Figure 2A displays the R-values obtained for the kpr with sustained signaling model. In this case the largest responses for the three systems are attained in steady state, with system 3 (the system with the largest values of k_on and k_off, black curve) producing the highest of the three R-values. Comparatively, we have that sustained signaling predicts much larger responses than those from others models¹, although the time required to reach a steady state is also much longer. This behavior also occurs in those models that include this effect (panel d of Figure 2A, and panels h and k of Figure 2B). The response for the kpr with limited signaling coupled to an incoherent feed-forward loop model shows similar characteristics to those of the simpler kpr with limited signaling model (compare panel l of Figure 2B, and panel c of Figure 2A), although the R-values computed from the kpr-iff model are smaller. As discussed in Appendix, the maximum steady state response for the kpl-iff model is not limited by T_T or P_T but by X_T. And although for comparative purposes we have considered P_T = X_T = 100, the lower responses of the kpl-iff model in comparison with those of the simpler kpl model are due to the presence of the iff loop and to the modulating effect of the activation chain C₀→C₁ → ⋯ → C_N on the signaling species X.

In any case, a clear pattern emerges from Figures 2A, B, namely: that, exception made of the occupancy model, systems with equal affinity show, however, quite different responses which would make affinity an unreliable estimate of the T-cell response. Further examples are given in Supplementary Material.

In Figures 3A, B we have illustrated the responses computed from the 12 phenotypic models in three systems that have different affinities and the following values of k_on and k_off:

where the value of k_on is given in subsection 2.2.

Panel a of Figure 3A shows that, as discussed in previous subsection and in Appendix, in the occupancy model responses increase always with affinity. Conversely, in the basic kpr model the system 3, with the smallest affinity, shows the largest response (panel b of Figure 3A, black curve), while in system 1, despite being the system with the largest affinity, practically no response is observed and its plot (blue curve) is almost coincident with the x-axis. Also, the following models: kpr with negative feedback (Figure 3A, panel e), kpr with induced rebinding (Figure 3A, panel f), kpr with stabilizing activation chain (Figure 3B, panel g), and kpr with negative feedback and limited signaling (Figure 3B, panel i), exhibit similar behavior to that observed for the basic kpr model, i.e., the largest response is obtained for the system with the lowest affinity (black curves) while null or very small values of R were found for the system with the largest affinity (blue curves). In turn, the kpr with limited signaling model displays a dual behavior (Figure 3A, Panel c): in transient phase the system with lowest affinity (black curve) shows the largest response, while in steady state the system with an intermediate value of A (red curve) gives the largest value of R. On the other hand, the system with the highest binding affinity practically shows no response at all times (blue curve). This behavior is also shown by the kpl-iff model (Figure 3B, panel l) although in this case the response for the system with the highest value of A coincides with the residual response in absence of kpr (see discussion in Appendix). In those models where sustained signaling is considered (panel d of Figure 3A and panels h and k of Figure 3B), we have that intermediate binding affinity gives the largest responses (red curves), the lowest affinity displays intermediate responses (black curves), while the system with the highest affinity (blue curves) shows the lowest R-values (although for the h and k models the R-values for the highest and lowest affinities are close, see Figure 3B).

Summarizing, we have that for the systems shown in this subsection only the occupancy model predicts a direct correlation between affinity and response independently of the values of k_on and k_off. In the remaining models that correlation does not exist and thus, there are systems with low or intermediate binding affinity exhibiting larger responses than those with a higher affinity, or systems with low and high affinities giving very close responses.

In Figures 4A, B we have displayed the responses computed from the 12 phenotypic models in three systems that have equal k_off and increasing values of k_on:

Because k_off is constant (τ = 10 s) the increase in affinity results only from the corresponding increase in the values of k_on. Hence and, unlike the two previous sections where the computed values of R as a function of affinity were due to variations in both k_on and k_off, the predicted responses in this section result only from the effect exerted by k_on (i.e., a single parameter), rather than by A. Therefore, since there are no trade-offs between k_on and k_off the data analysis of the corresponding responses is facilitated. Thus, in agreement with previous results, the occupancy model shows (Figure 4A, panel a) that under these conditions an increase in the values of k_on (and therefore, in affinity) produces an enhancement of the response. The same behavior is observed in the remaining panels of Figures 4A, B with the exception of the models in which negative feedback is involved (Figure 4A, panel e; Figure 4B, panel i), where systems with the largest values of k_on display the lowest responses (black curves). However, additional results for a larger value of τ (=100 s) given in Supplementary Material reveal that all models, including those with negative feedback, now show an increase of the response as k_on becomes larger. Further analysis of the results displayed in Figures 4A, B and in Supplementary Material also show that, whereas the remaining parameters are maintained constant, responses from all models become independent of the k_on-values, and therefore of affinity, for k_on>>1. In short, we have that all models (with the exception of the kpr with negative feedback model when τ = 10 s) predict that when k_off (or τ) is constant, the R-values increase with k_on and therefore, with affinity. But strictly speaking, what this data analysis really reveals is the effect exerted by variations in the values of the rate constant k_on on the predicted immune responses when the others parameters remain constant.

Immune responses exhibit a remarkable dependence on the values of k_off(= 1/τ) (4, 37, 38, 56) and, in fact, it has been suggested that the value of k_off is the best estimate of T-cell activation (19, 57–59). To explore in further detail this dependence we have displayed in Figures 5A, B the responses computed from the 12 phenotypic models for five systems having equal k_on and decreasing values of k_off (increasing values of τ) for a relatively narrow interval of values of τ and A:

In turn, Figures 6A, B display the plots obtained for systems with a much wider range of τ and A-values:

In Figures 5A–6B the value of k_on is constant so that increased affinity is only caused by the corresponding decreasing in the values of k_off. Thus, we should bear in mind that plots in these Figures display in fact the effect exerted by a single parameter, k_off (or τ), rather than by A, on the predicted responses (in examining Figures 5A–6B, note that for some models and values of τ the corresponding R-values are so small that their plots are not observed because they almost coincide with the x-axis).

Panels (a) in Figures 5A, 6A display the R-values computed according to the occupancy model for the narrow and wide ranges of values of τ. Plots in these panels illustrate that, as already discussed for this model, values of R always increase as affinity becomes larger since they do not depend on the individual values of k_on and k_off. The same behavior is observed for the basic kpr model [panels (b) in Figures 5A, 6A].

Note, however, that the behavior shown by the basic kpr model in Figures 5Ab, 6Ab, namely that an increased affinity as result of an increase in the τ-values when k_on is constant leads to larger responses, contrasts with the behavior described previously for the same model in subsection 3.2, where systems with larger affinities produced lower values of R (panel b in Figure 3A). This is due to the fact that affinity is not a parameter but a quotient of two parameters, so that variations in the values of A can occur in one of the following three ways: (a) changes in the values of k_on while k_off remains constant; (b) changes in the values of k_off while k_on remains constant; and (c) changes in the values of both k_on and k_off. However, these three ways of changing affinity are not equivalent and produce different responses. This is shown in the following simple example demonstrating that the R-values are quite different when a 5-fold increased affinity is caused by a 5-fold increase in the value of k_on while k_off remains constant, or when the same increased affinity is obtained by dividing by five the value of k_off while k_on is constant (other conditions as in Figure 2A):

Data analysis of R-values computed from other models as a function of τ is more complex. Thus, for the kpr with limited signaling model the R-values increase with affinity for the smaller τ-values both in transient phase and in steady state (panel c in Figure 5A). However, for the largest values of τ (panel c in Figure 6A) the R-values increase with affinity reaching peak values in the transient phase, but as the time activation progresses the curves cross over so that in steady state the system with the largest affinity exhibits the lowest response². The same behavior is also exhibited by the kpr with limited signaling coupled to an incoherent feed-forward loop model (panels l in Figures 5B, 6B). However, for the following three models, kpr with negative feedback, induced rebinding, and stabilizing activation chain (panels e, f in Figures 5A, 6A, and panel g in Figures 5B, 6B), the corresponding responses show a similar behavior to that observed for the simpler occupation and basic kpr models, i.e., that the R-values increase with affinity for all values of τ when k_on is constant. Regarding the kpr with sustained signaling model, the corresponding responses are more elaborate: the values of R increase with affinity for the smaller values of τ (panel d of Figure 5A), while conversely for the larger values of τ the opposite effect occurs (panel d of Figure 6A).

Finally, in the mixed models (h)-(k) the responses obtained are determined by the parameters and new assumptions added to the basic kpr model. For example, for the kpr with limited and sustained signaling model (where both effects modulate the response), the largest value of R in panel h of Figure 5B is not obtained for the system with the largest affinity (A = 20k_on, black curve) but for the system with A = 5k_on (magenta curve). This contrast with the behavior observed under the same conditions for the kpr with only limited signaling model (Figure 5A, panel c), and the kpr with only sustained signaling model (Figure 5A, panel d) where the largest response in both cases was obtained for the largest affinity (black curves). These mixed characteristics are also observed quantitatively: the values of R in the mixed model (panel h in Figure 5B) are much larger than in Figure 5Ac (only limited signaling), but lower than in Figure 5Ad (only sustained signaling). Responses exhibiting combined behaviors are also observed for the other mixed models (panels i, j, k) in Figures 5B, 6B.

It is also worthy of note that comparison of responses from this and previous section shows that the influence exerted by the rate constant k_off is more involved than that exerted by k_on. This is not surprising since Figures 1A, B show that the role of k_on in the reaction schemes of all models is limited to the first step of the activation chain, while conversely k_off is involved in most of the steps of the activation process. This is in agreement with experimental studies which had revealed the predominant influence exerted by τ on T-cell responses (19, 57–59).

Summarizing: results in this and previous sections allow to assess that the relationship between affinity and predicted T-cell responses is highly complex so that there is no a general and simple correlation between affinity (which is the quotient of two parameters) and response. Obviously, this does not exclude that for some systems with a particular set of parameters, or systems whose parameters are within a particular interval, a certain degree of correlation could be found.

TCR-pMHC affinity (= k_on/k_off) is independent of the ligand concentration. However, it is expected that quantitative values of T-cell responses depend on the values of P_T and T_T, the total amount of pMHC and TCR respectively. Therefore, the question arises as to whether data analysis from experimental studies aimed at establishing a correlation between affinity and T-cell response could be influenced by the values of P_T used in the assays. To this end, we have reproduced Figures from previous sections using a different P_T-value (PT=2×104 instead of P_T = 100), and the results are shown in Supplementary Material (SM). Thus, for systems with the same affinity but different values of k_on and k_off (Figures 2A, B in main text and the corresponding Figures 2AS, 2BS in Supplementary Material), in general we find that the behavior of the responses obtained from most models for both values of P_T were similar to those already discussed in subsection 3.1 of main text, although quantitative responses were much larger in the Figures displayed in SM. However, interesting exceptions were also found. For example, the kpr with sustained signaling model predicts that for the lower value of P_T responses for the three systems with the same affinity but different values of both rate constants are also different (panel d in Figure 2A in main text). However, the same model for PT=2×104 shows that for the three systems with equal affinity the values of R are now practically equal (panel d in Figure 2AS in Supplementary Material). A superficial observation could state that “systems with equal affinity give the same response.” This, however, is a false conclusion because what this model really predicts is that under these conditions (a high value of P_T) the three systems saturate and the maximum responses are attained independently of the values of affinity.

At the other end, we have that the kpr with negatived feedback model predicts extremely low or null responses for this high value of P_T (compare panels e of Figure 2A, Figure 2AS in main text and SM respectively). Again, this is not related to affinity but with the fact that for some values of the parameters (among them the P_T-values) this mechanism can act as a switch so that for parameter-values below this switch the signaling chain goes forward (larger signaling) while above goes backward (no response).

In turn, inspection of the responses for the kpr with limited signaling coupled to an incoherent feed-forward loop model shows (compare panels l in Figure 2B, Figure 2BS): (a) that for P_T = 100 the values of R in steady state for the three systems increase in this order blue → black → red, while for PT=2×104 we have red → black → blue, i.e., the R-values are reversed; (b) that unlike other models which predict much higher responses when PT=2×104, the R-values in panel l of Figure 2B (P_T = 100, main text) and in panel l of Figure 2BS (PT=2×104, SM) are always < 100. The reason, as shown in Appendix, is that responses for this mechanism are not bounded by P_T or T_T but by X_T, the total amount of species X. Hence, and because in our modeling we have assumed X_T = 100, responses for this model can't be larger than this value independently of the values of P_T, T_T, or affinity.

Regarding the models built by adding assumptions and parameters from others models (models h-k) we notice that, effectively, responses can be modulated by these assumptions and parameters. Thus, for example, for the kpr with negative feedback and limited signaling model we find that for PT=2×104 a null response is obtained (compare panels i of Figure 2B, Figure 2BS in main text and in SM) which reveals the strong effect exerted by the negative feedback process under these conditions (see above). Likewise, the kpr with stabilizing activation chain and sustained signaling model also shows the modulating effect exerted by the sustained signaling process (compare panels k of Figure 2B, Figure 2BS in main text and in SM).

Similar considerations apply also for the remaining Figures with PT=2×104 displayed in SM, and for this reason further discussion is not given here.

In previous sections we have highlighted that to study without ambiguity the influence of a given parameter on the T-cell response we should analyze that response by varying only this parameter while all other parameters of the system are kept constant. Thus, if T-cell responses as a function of affinity are measured under conditions in which other system parameters also vary there will be cross-effects of these parameters on the observed response that mask the influence exerted by affinity. This could lead to misinterpret the results obtained and lead to false conclusions. In this section we illustrate a representative example of this situation by considering three different systems that have equal values of k_on and k_off and, therefore, of affinity, but they have different values of the rate constant along the proofreading activation chain, k_p:

By assuming that we limit ourselves to consider T-cell responses and affinity, the responses from the three systems should be the same because they have equal affinity. For the occupancy model (Figure 7A, panel a) this is really the case because in this model only the binding process (which is independent of the rate constant k_p) is involved. Hence, and since the values of k_on, k_off and A are the same for the three systems, this model predicts that their responses are also equal and the three curves are overlaid. However, for the remaining models the predicted responses are quite different (Figures 7A, B, panels b-l) which, obviously, is due to the fact that computed responses were not obtained as previously in Figures 1–6, i.e., with all other system parameters remaining constant. This reveals that if T-cell measurements are analyzed without taken into account properly the parameter space it could lead to erroneous conclusions.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.