The mechanism that enables T cells to process signals so remarkably well has remained elusive-partially because equilibrium thermodynamic processes are insufficient to explain ligand discrimination [4,6,7]. To address the discrepancy, it was suggested that the TCR uses kinetic proofreading (KP) [8] to discriminate binding events of different duration [9]. Briefly, KP occurs when ligand binding to a receptor triggers a series of irreversible, energy-consuming biochemical modifications of the receptor until a final, signaling-competent state is reached, while intermediate states do not relay the signal. If ligand unbinding results in modifications being removed, this means that only long binding lifetimes result in receptor signaling. Applied to TCR activation through multisite phosphorylation, KP can account for TCR discrimination; KP models have been found to fit experimental TCR activation data better than models without proofreading steps [5,10,11]. For a discussion of the "lifetime dogma" of TCR activation, see Ref. [12].
Independently of the physical basis assumed for KP-aided ligand discrimination, all KP steps are downstream of initial ligand binding to the TCR. In other words, the TCR must first sense that ligand binding has occurred before KP-based ligand discrimination may occur, even if the modifications underlying KP can be made to unbound TCRs. Despite much research, the mechanism of TCR triggering is also still an unresolved question. Proposed models have either evoked conformational change or aggregation of TCRs to explain the ligand-induced shift in TCR phosphorylation or a mechanism referred to as kinetic segregation (KS) model of TCR triggering [4,7]. The KS model proposes that TCR phosphorylation and activation is determined by the TCR's residence time in so-called close contacts between T cells and APCs, when these cells bring their membranes into close proximity to allow for TCR-pMHC binding [13,14]. Due to the tight intermembrane spacing, these close contacts are depleted of deactivating phosphatases with large extracellular domains, such as CD45, but retain the TCR kinase Lck needed for TCR phosphorylation and activation [15][16][17]. TCRs will activate if they remain in a close contact for a sufficiently long time so that they are fully phosphorylated, with ligand binding only serving to enhance residence times by preventing TCRs from diffusing out of the close contact. Among other work supporting the KS model, TCR elongation was shown to prevent TCR signaling [18], while TCRs were observed to activate ligand independently within sufficiently large close contacts with a decreased CD45-to-Lck ratio [17].
Importantly, the two questions of how TCRs are triggered and how they discriminate between ligands are connected to each other: the assumptions made in any TCR triggering model will impact how KP could be physically implemented in TCR triggering. Yet, there is little literature over how a particular triggering model, e.g., TCR triggering modeled by KS, will affect KP-based ligand discrimination by the TCR. Classic formulations of KP in TCR signaling models usually assume that the TCR is "reset" to its initial, unphosphorylated state by ligand unbinding. When incorporating the KS model of TCR triggering, one needs to assume that phosphorylation can occur as long as the TCR is in the close contact, whether ligand bound or not, and that TCR dephosphorylation and thus resetting only occurs when a TCR leaves the close contact. The mathematical implementation of KP is thus directly linked to the triggering model assumed and may affect overall model performance. Moreover, modeling KS-driven TCR triggering also requires explicit modeling of kinase-enriched membrane areas, close contacts, and TCR diffusion (see also, e.g., Ref. [19]), which makes it hard to directly compare models of KS-based TCR triggering to other models that often do not model TCR diffusion and, if they do, not in the context of close contacts.
Here, we address this knowledge gap by developing a computational framework for describing a combined model for TCR triggering and ligand discrimination: a classic, multistep downstream signaling model implementing KP combined with an explicit model for TCR triggering inside the close contacts postulated by the KS triggering model. While this framework is tailored to the KS model of TCR triggering (our KS-KP model), we can easily modify the model to assume that instead of by KS, the TCR is strictly triggered by ligand binding only, e.g., based on conformational change of the TCR (our CC-KP model). Our new approach therefore allows us to directly compare CC-KP and KS-KP models of TCR activation and test how performance differs for those two models. We find that while both models predict similar trends of TCR activation for a wide range of parameters, the KS-KP model predicts a much higher "background" probability of activation than the CC-KP model in the absence of any activating ligands for realistic TCR phosphorylation rates. While both CC-KP and KS-KP predict similar patterns of TCR activation, we also find differences in several biologically important situations, such as when small numbers of activating ligands are present or for ligands with comparatively larger off rates. To reach similar discrimination than CC-KP, KS-KP needs a large number of KP steps; moreover, when looking at the trade-offs between discrimination and sensitivity, CC-KP outperforms KS-KP with the notable exception of slower phosphorylation rates. While we discuss some experimental evidence to support KS-KP over CC-KP, we also conclude that CC-KP and KS-KP models are difficult to distinguish through the mere comparison of downstream metrics of T cell activation in the absence of further biochemical data on signal transduction.
To allow a direct comparison between a conformational change-kinetic proofreading KP (CC-KP) and a kinetic segregation-kinetic proofreading (KS-KP) model (Fig. 1), we developed a mathematical framework and computational model that is described in detail in the Materials and Methods section. In this framework, we explicitly consider the close contacts formed between T cells and APCs. We allow individual TCRs to diffuse into and out of the close contact by Brownian motion, binding to and unbinding from both activating and nonactivating ligands while being in the close contact. (While this terminology is often used in the literature, it is possible that a nonactivating ligand can activate a TCR, although at a much lower rate than an activating ligand.) We assume that ligand bound TCRs are unable to leave the close contact, since we define the close contact as the region in which the cell membranes of the two contacting cells are in sufficiently close apposition to allow a pMHC-TCR complex to form across the two contacting cells.
In our model, we consider a T cell to be activated once all n proofreading steps are completed for at least one TCR. Consistent with recent work, we consider small values of n, which are known to be consistent with ligand discrimination [5,20]. Moreover, while we model these proofreading steps as phosphorylation reactions, we note that these steps are not intended to directly correspond to specific TCR phosphorylation sites but could also include phosphorylation sites on Zap70 or LAT [21], in line with recent literature suggesting that the number of proofreading step may not directly correspond to all steps in the biochemical network [5,22]. In the CC-KP model, ligand unbinding results in immediate reset of the proofreading cascade (full dephosphorylation of all sites) whereas in the KS-KP model, ligand unbinding resets proofreading if and only if the TCR has subsequently also diffused out of the close contact.
Unlike earlier related models [19,23], we modify the (de)phosphorylation rates depending on whether the ligand is bound or unbound and depending on which version of KP we use. In CC-KP, the phosphorylation rate of unbound TCR is zero, and the dephosphorylation of unbound TCR rate is effectively infinite, resulting in an immediate reset after ligand unbinding. Meanwhile, in KS-KP, the (de)phosphorylation FIG. 1. Overview of our computational models of T cell receptor (TCR) activation in close contacts. (top row) The conformational changekinetic proofreading (CC-KP) model is depicted with three phosphorylation sites (n = 3) in the example. A TCR can be phosphorylated when bound to nonactivating or activating ligands and is rapidly dephosphorylated upon ligand unbinding (scheme, left). The signaling active state (all three sites phosphorylated) is highlighted in blue. (middle and right columns) An unproductive trajectory (on which the TCR does not reach the active state before leaving the close contact) and productive trajectory, respectively. Note that in the unproductive trajectory, the TCR only binds to nonactivating ligands, with a rapid off rate preventing TCR phosphorylation. In the productive trajectory, a long period of time where the TCR is ligand bound results in all three sites being phosphorylated while the TCR is immediately fully dephosphorylated upon ligand unbinding. (bottom row) The kinetic segregation-kinetic proofreading (KS-KP) model is shown in the bottom half of the figure. The diagram on the left shows how, with the KS-KP model, phosphorylation can take place whether a ligand is bound or not and can only be reversed when the TCR leaves the close contact (depicted with dashed arrows). Accordingly, in both the unproductive and productive trajectories shown for KS-KP, phosphorylation events occur for both bound and unbound TCRs. KS-KP is still sensitive to ligand concentrations because ligand binding influences the time that a TCR spends in the close contact before leaving. In this figure , M representing the off-rates); and k p represents the phosphorylation rate. S1 for the parameters used in this figure . rate is independent of the binding state in the close contact, and the timescale for reset depends on binding only insofar as binding increases the TCR's residence time in the close contact, as effective dephosphorylation only occurs once a TCR has diffused out of the close contact. (As a check of the accuracy of our diffusive model, in the SI we show that the residence time in the close contact is dependent on the off rate of the activating ligand, indicating that the binding dynamics of the ligand modify its diffusion, as expected; see Fig. S1 of the Supplemental Material [24].
As we model only the very first steps of TCR activation, we assume that TCRs do not interact with each other. While formation of TCR clusters is known to play an important role in TCR signal integration [25,26], experimental data also suggest that triggering a single TCR is sufficient to activate a T cell [1,2], consistent with our assumption that the T cell is activated once TCR phosphorylation exceeds a certain threshold.
Finally, as a check for consistency of our framework with previous models, we calculated the probability of a single TCR remaining in the close contact as a function of time; we then compare it with the same quantity calculated from a partial differential equation (PDE) model adapted from Fernandes et al. [23]. The simulation results match reasonably well with those from the PDE model (Fig. S2 of the Supplemental Material [24]). We may therefore assume that the results we obtain are comparable with results obtained from earlier related models.
We first investigated how the number n of proofreading steps affects the activation probability for different activating or nonactivating (activating or self) ligand ratios in both CC-KP and KS-KP models at a fixed phosphorylation rate (k p = 1 s -1 [5]), using values for n in the range of 2-4 (Fig. 2).
We note that we use values of n in the range that was previously found to well explain T cell activation data [5] and that are optimal for ligand resolution in the presence of molecular noise [22].
Generally, the activation probability predicted by KS-KP was much more strongly dependent on n as that predicted by the CC-KP model. We can rationalize these findings in the light of the overall increased activation probability seen for the KS-KP model. This "background activation" stems from ligand-independent triggering of TCRs that, even in the absence of a ligand, show a close-contact residence time sufficient to complete all KP steps. Thus, it is logical that for KS-KP and a small number of KP steps, we observed higher levels of background activation than for the CC-KP model. Interestingly, however, very small changes to n (±1 step) and/or k p (factor of two or less) were sufficient for KS-KP and CC-KP to yield highly similar results (inset to the n = 4 panel of Fig. 2 and Fig. S3 of the Supplemental Material [24]). Taken together, these results suggest that CC-KP and KS-KP may be difficult to distinguish through downstream metrics of T cell activation, such as, e.g., in Ref. [5]. In contrast, quantification of the ligand-independent contribution to TCR triggering would allow us to further test our KS-KP model and to discriminate between CC-KP and KS-KP models.
When we varied the phosphorylation rate k p more broadly, we found that for a lower phosphorylation rate (k p = 0.5 s -1 ), KS-KP was able to distinguish low and high G values better than CC-KP for n = 3 and n = 4 (Fig. S4 of the Supplemental Material [24]). Meanwhile, at a higher phosphorylation rate (k p = 2 s -1 ), KS-KP was entirely unable to distinguish between different activating fractions at n = 2. This is again in line with the expectation that in KS-KP, ligand-independent phosphorylation of the TCR contributes to a higher sensitivities at slow rates for the KP steps, helping to "push" TCR across the signaling threshold. Conversely, for fast KP rates this ligand-independent phosphorylation results S1 for parameters used in this figure [24]. in all TCR being immediately phosphorylated and hence loss of discrimination.
Finally, we note that in our model, TCR phosphorylation times (for each step) are exponentially distributed (we assume first-order kinetics). As we assume that each of the n steps occurs at a similar rate, it is reasonable to approximate the total time (i.e., the time until full phosphorylation or activation time) as an Erlang-n distribution [see Materials and Methods; Eq. ( 25)]. We note that earlier work by Fernandes et al. [23] used a simplified model for TCR activation, in which any TCR that remained in the close contact for at least 2 s was immediately and fully activated, i.e., the effective probability density function for the activation time has a Dirac δ function at t = 2 s. In such a simplified model, activation probabilities for KS-KP are markedly underestimated at low fractions G of activating ligand, implying better discrimination by KS-KP than are likely realistic. Further information on this point is found in the Supplemental Material (SM) (Fig. S9 and SM text section S1 [24]).
Activating ligand off rates have been shown to correlate with TCR activation [7,27,28], a fact that has been used as evidence of KP. We therefore wanted to test how activating ligand off rates influence the performance of our models. We used activating ligand off rates of k - a ∈ {0.1 s -1 , 1 s -1 , 10 s -1 }, and we use a nonactivating (self) ligand off rate of k - n = 50 s -1 for all simulations; that is, the weakest binding activating ligand we examine (k - a = 10 s -1 ) has a fivefold smaller dissociation constant than the nonactivating ligand [23]. We first tested how activation probability changed as a function of activating fraction G for different close contact lifetime t c and different k - a . For a small t c , activation probability always varies with k - a , with smaller k - a giving rise to more activation as expected, but there is little difference between the activation predicted by either KS-KP or CC-KP models (Fig. 3 As before, these findings can be rationalized by considering the contribution from ligand-independent triggering to the activation probability, only present in the KS-KP model and significant for close contact lifetimes substantially longer (t c = 30, 120 s) than the time taken to activate (Erlang distribution around 2 s). For none of the close contact lifetimes t c tested, KS-KP or CC-KP can reliably discriminate ligands with the lowest k - a = 10 s -1 or a fivefold smaller dissociation constant than the nonactivating ligand for the ligand fractions tested (G 0.1). The ligand-independent contribution to activation is also readily visible when plotting activation probability as a function of close contact lifetime t c (Fig. 4): in contrast with CC-KP, KS-KP predicts nonzero TCR activation even in the absence of activating ligand (G = 0) that increases as t c increases, as expected independently of ligand off rate k - a [Fig. 4(a)]. Interestingly, the extent of this ligand-independent activation seen in the KS-KP model is similar to the level of activation predicted for a ligand of k - a = 1 s -1 at G = 0.01 in the CC-KP model [4(b)]. For larger G (e.g., G = 0.1) and slow k - a , predictions by CC-KP and KS-KP differ little but they remain marked for faster values of k - a [Fig. 4(c)], consistent with the previous results.
In summary, while we find similar overall behavior between the two models as we change the activating ligand off rate (k - a ), it is noteworthy that the KS-KP model predicts close contact lifetime-dependent ligand-independent activation probabilities (Fig. 4) in contrast with the CC-KP model.
While ligand-independent activation is qualitatively consistent with experimental observations of activation in the absence of activating ligands [17,18], we expect it to negatively affect ligand discrimination (by adding spurious activation in the absence of ligand) while sensitivity could be positively affected (increasing background phosphorylation of S1 for parameters used in this figure [24].
the TCR can lower the threshold for its activation). Thus, we tested the trade-offs made between discrimination and sensitivity by KS-KP and CC-KP models.
We define sensitivity as the derivative of P a with respect to
We make the dependence of P a on G explicit here, although in general we do not write out this dependence for the sake of simplifying notation. This quantity can be interpreted as the ability of the T cell to respond to introduced activating ligands. We define discrimination as
describing how well the T cell can discriminate between high and low activating ligand concentrations. If α is close to zero then the T cell will likely activate even at low activating ligand concentrations, while if α is close to one the T cell will not activate at low concentrations but will activate at high enough G values. See the Methods section for more detail. We investigated for both models how sensitivity and discrimination change as a function of KP steps (Fig. 5). As to be expected, increasing the number of KP steps n from one (i.e., no KP) to two always increases both sensitivity and discrimination. After this, discrimination keeps increasing further with every added KP step, but the trade-off between sensitivity and discrimination is different for KS-KP and CC-KP. At some point, for CC-KP, increasing n results in marginal improvements in already high discrimination for large losses in sensitivity. In contrast, the trade-offs for KS-KP continues to involve large gains in discrimination for small losses in sensitivity before nearing perfect discrimination and low sensitivity at larger n. KS-KP only outperforms CC-KP in sensitivity for slower k p rates: for k p = 0.5 and n = 3, KS-KP is more sensitive with moderate discrimination loss [Fig. 5(a), bottom row]. For k p 0.5, CC-KP is both more sensitive and discriminate than KS-KP at equal n and hence equal energetic cost (Fig. 5).
To compare our theoretical results for the CC-KP and KS-KP models to experimental data, we used published EC 40 and EC 50 values for several ligands with known on and off TCR binding rates [27,28] (see Table S2 [24]). In aggregate, these experiments provide information about the overall T cell activation probability for three nonactivating ligands and four activating ligands each in the absence of any other ligand. To compare the CC-KP and KS-KP models to these results, we ran simulations and calculated activation probabilities as before, using the corresponding experimental on and off binding rates for each ligand. For nonactivating ligands we used G = 0 and for activating ligands (weakly activating or otherwise) we used G = 1, while varying the overall ligand concentration, M. We estimated EC 40 and EC 50 values (according to the quantity available for each ligand, from here on referred to as EC 40/50 ) by finding the value of M for which the activation probability was 40% or 50% maximal, respectively. While we do not expect our proxy for T cell activation (P a ) to quantitatively predict experimental EC 40/50 values, it should be possible to qualitatively predict the correct order of the ligands. Both the CC-KP and KS-KP models were able to distinguish activating from nonactivating ligands, with lower theoretic EC 40/50 values predicted for all activating ligands (Fig. 6). However, remarkably, there is a difference between KS-KP and CC-KP in ordering of the ligands, and particularly, the KS-KP model differentiates between the activating and weakly activating ligands better than the CC-KP model. The ordering of, e.g., T72 and OVA is also more accurate by the KS-KP model. These differences can be rationalized by KS-KP results being dependent on the dissociation constant, FIG. 5. Trade-offs between discrimination and sensitivity in the number of proofreading steps are more pronounced for KS-KP than for CC-KP. Sensitivity and discrimination are defined in the Methods section, equations (37)-(39). Blue points are for CC-KP while orange points are for KS-KP. The leftmost points of the KS-KP and CC-KP points with the lowest discrimination are for n = 1 proofreading steps (i.e., no proofreading), while the rightmost points with the highest discrimination are for n = 8. For each column, with (a) k p = 0.5 s -1 , (b) k p = 1 s -1 , and (c) k p = 2 s -1 , the plot in the bottom row shows the trade-off region highlighted in the box within the top row. In the bottom row, the lowest and highest number of steps shown in the trade-off region are indicated. Below the smallest n in the trade-off region, increasing n results in both higher discrimination and sensitivity. For all three values of k p shown, the trade-offs for CC-KP in n involve marginal improvements in high discrimination for large losses in sensitivity. In contrast, trade-offs for KS-KP involve large gains in discrimination for small losses in sensitivity at lower n before nearing perfect discrimination and low sensitivity at large n. See Table S1 for the parameters used in this figure [24].
while CC-KP is more sensitive to the off rate. In other words: the duration of an active period for KS-KP only depends on FIG. 6. Comparison of experimental EC 40 and EC 50 values for ten ligands with theoretical predictions from the CC-KP and KS-KP models. The model EC 40 and EC 50 values are the ligand concentrations at which activation probability is equal to 0.4 and 0.5, respectively, with G = 0 for nonactivating ligands and G = 1 for (weakly) activating ligands to match experimental conditions where a single ligand type is present. Both the CC-KP and KS-KP models differentiate nonactivating from activating ligands, separating them into two clear clusters along the y axis. The KS-KP half-maximal activation probability ligand concentrations are lower than the CC-KP predictions, consistent with the results in previous sections. The ten ligands considered, along with the corresponding binding and unbinding rates, are listed in Table S2 [24]. See Table S1 for the parameters used in this figure [24]. how long a TCR stays in the close contact, which correlates to how often it is bound. CC-KP, however, solely depends on the duration of binding, which is controlled by the off rate.
While both CC-KP and KS-KP predict similar patterns of TCR activation, we also find differences in several biologically important situations, such as when small numbers of activating ligands are present or for ligands with comparatively larger off rates. To reach similar discrimination than CC-KP, KS-KP needs a large number of KP steps and when looking at the trade-offs between discrimination and sensitivity, CC-KP outperforms KS-KP with the notable exception of slower phosphorylation rates. We also conclude that CC-KP and KS-KP models are difficult to distinguish through the mere comparison of downstream metrics of T cell activation in the absence of further biochemical data on signal transduction. For example, recently published experimental work that has been interpreted as evidence for CC-KP [10,11] would also be consistent with KS-KP if, under the conditions used for these studies, phosphatases are excluded from regions of TCR-ligand interactions. Other experimental evidence indicates that TCRs may activate even with no activating ligands present [17,18], which is more likely under KS-KP than CC-KP, but the high specificity seen in Refs. [1,2] is more readily explained by CC-KP. An advantage of our computational approach is that both KS-KP and CC-KP might be viewed as extreme points on a spectrum of possible behaviors: for example if ligand binding changes k p but is not essential for it, such a model would exist between KS-KP and CC-KP, and could be evaluated in our framework in future studies.
In this work, we have developed a computational framework for describing a combined model for TCR triggering and ligand discrimination: a classic, multistep downstream signaling model implementing KP combined with explicitly modeling TCR triggering inside the close contacts postulated by the KS triggering model. We have used this to compare the KS model of TCR triggering (our KS-KP model) to one where the TCR is strictly triggered by ligand binding only (our CC-KP model). Considering existing models in the literature, our KS-KP model is generalized so that phosphorylation rates while bound can be different from those while unbound (compared with Ref. [19]), and has a more accurate model of activation time (compared with Ref. [23]), while applying a diffusion framework to allow direct comparison between CC-KP and KS-KP. Integrating diffusion with conformational change-based KP in the CC-KP model captures aspects of T cell activation not found in KP models ignoring diffusion. If TCR complexes cannot leave the close contact when bound to a ligand, trajectories that happen to spend more time in a bound state will be more likely to remain in the close contact longer than trajectories that are in the bound state less often. This has a selective effect where trajectories with higher individual activation probabilities, P s (t ), are more likely to remain in the close contact longer (Fig. S10 [24]). However, because most TCR complexes rapidly exit the close contact, the overall probability of a single TCR complex activating is greatly reduced in the CC-KP model with diffusion, as compared with a model ignoring diffusion. This complex coupling between diffusion and activation necessitates a model explicitly including both aspects.
Over a wide parameter range, CC-KP and KS-KP models make similar predictions. This is due to the fact that the kinetic segregation mechanism in KS-KP also implements a form of kinetic proofreading, where a reset to the fully dephosphorylated state occurs upon the TCR complex diffusing out of the close contact (instead of upon ligand dissociation, as in CC-KP). Our modeling of CC-KP and KS-KP in a close-contact zone context suggest that all else being equal, KS-KP leads to higher activation probabilities than CC-KP along with sizable, nonzero activation probabilities without activating ligands present. This positive minimal activation probability results from any TCR increasing the activation probability as long as k p is sufficiently large to drive TCR phosphorylation while in a close contact, and it results in worse differentiation between activating ligands with varied off rates in the KS-KP than CC-KP model. However, we note that for lower phosphorylation rates (k p = 0.5 s -1 ) KS-KP performed better than CC-KP in terms of sensitivity. For k p = 0.5 s -1 , both sensitivity and discrimination were better for KS-KP, as for this rate, full phosphorylation of the TCR (and hence activation) cannot be reached before ligand dissociation in CC-KP. Conversely, at fast phosphorylation rates (k p = 10 s -1 ), KS-KP was entirely unable to distinguish between different activating fractions for all n tested as at this rate all TCRs are fully phosphorylated before they leave close contacts, even for residence times expected for free TCRs. For other parameters, KS-KP and CC-KP are generally similar at low t c , low-to-moderate k - a , or moderate-to-high G, while significant differences are seen at high k - a and low G, with these differences being largely due to ligand-independent activation present in the KS-KP model.
In this light, contemporary experimental evidence observed for different TCR signaling models is hard to interpret unambiguously but points in favor of KS-KP. For example, recently published experimental work that has been interpreted as evidence for CC-KP [10,11] would also be consistent with KS-KP if, under the conditions used for these studies, phosphatases are excluded from regions of TCR-ligand interactions. Other experimental evidence indicates that TCRs may activate even with no activating ligands present [17,18], which is more likely under KS-KP than CC-KP, but the high specificity seen in Refs. [1,2] is more readily explained by CC-KP. An advantage of our computational approach is that both KS-KP and CC-KP might be viewed as extreme points on a spectrum of possible behaviors: for example if ligand binding changes k p but is not essential for it, such a model would exist between KS-KP and CC-KP, and could be evaluated in our framework in future studies.
Throughout our modeling, we assume that the timescale of immune synapse rearrangement (minutes) is different to the dynamics in the KS model that result in initial TCR triggering (milliseconds to seconds). Close contacts, i.e., the region of phosphatase exclusion, are expected to grow on the timescale of TCR triggering and diffusion, and the effects of modeling this growth similar to Ref. [23] on CC-KP and KS-KP predictions could thus be explored in future work.
In conclusion, our results may partly explain the ambiguity surrounding the models for TCR activation, as without direct evidence for the number of phosphorylation steps and the rate of phosphorylation and without experiments that recreate the close contact context in which TCRs are activated, our results suggest that distinguishing between the KS-KP and CC-KP models through downstream measures of TCR activation may be difficult. On the other hand, our results point to clear differences that carefully designed experiments could measure. For example, if TCR activation and discrimination was measured for different phosphorylation rates, e.g., making use of kinetic mutants of Lck [29,30], ligand discrimination should be much more affected for CC-KP than KS-KP. These promising potential directions of future experimentation, combined with consistency between our model and experimental findings on TCR activation, is encouraging to future experimental work. We hope that our results will facilitate further investigation into TCR signaling mechanisms, perhaps leading to a synthesis of two views of T cell activation (KS-KP and CC-KP) that are often viewed as being mutually exclusive.
Our goal is compare the consequences of conformational change and kinetic segregation TCR triggering mechanisms. We are able to directly compare the two mechanisms by simultaneously modeling TCR diffusion in a disk, representing a close contact between a T cell and an APC, and TCR activation. For both triggering mechanisms, we assume that TCRs can only trigger when inside the close contact and that they cannot leave the close contact while bound to a ligand. Ultimately, we wish to calculate the probability that at least one TCR will trigger, given the lifetime of the close contact, the TCR concentration, the ligand concentration, and the phosphorylation and binding rates.
We accomplish this task by first simulating TCR diffusion and ligand binding-unbinding. The diffusion-bindingunbinding dynamics are identical for both triggering mechanisms, so we can use the same simulations for both. Although many TCRs can enter and leave a close contact through diffusion, we assume that TCRs diffuse independently before activation, allowing us to combine the results of single TCR simulations. We initialize a diffusion-binding-unbinding simulation with an unbound TCR close to the close contact boundary. Using a fixed time step, we choose the next location of the TCR according to a normal distribution and the binding state of the TCR according to a Markov chain determined by nonactivating and activating binding and unbinding rates. If the TCR is in a bound state at a time step, we disallow it from leaving the close contact by diffusion in the next time step. This coupling between diffusion and binding complicates modeling the conformational change mechanism using partial differential equations. For KS-KP, the total time a TCR spends in the close contact determines the probability of triggering, while for CC-KP the history of ligand binding times is needed to calculate triggering probability. The diffusion-binding coupling prevents simple calculation of the binding time distribution from PDE calculations, which we instead use to check that the stochastic simulations are correct.
We can numerically calculate the KS-KP and CC-KP triggering probabilities for each simulation replicate and average over all replicates to find the probability that a single TCR will trigger. We can then calculate the probability that at least one TCR will activate in the close contact lifetime, given TCR concentration and the size of the close contact.
We wish to directly compare the behavior of a typical KS-KP kinetic segregation model with that of a CC-KP model in the context of a T cell-APC close contact. In KS-KP models, a single TCR can be phosphorylated as long as it is in the close contact, regardless of whether it is bound to a ligand or not, and dephosphorylation occurs only upon leaving the close contact. In CC-KP models, TCRs can be phosphorylated only when bound to a ligand and are rapidly dephosphorylated upon ligand unbinding. To capture the behavior of both KS-KP and CC-KP models, we used stochastic simulations of a TCR in a close contact. We simultaneously modeled Brownian motion of the TCR and its ligand binding dynamics, where interaction between the two processes has significant effects on activation probability for CC-KP and no effect on KS-KP activation beyond extending close contact lifetimes. We explain this point further below. The output of each simulation described in this section is a trajectory of the two TCR spatial coordinates and the binding state, stored in a 3 × N T matrix, where N T is the number of sampled times.
Within the close contact, we model TCR motion according to
where x(t ) is the two-dimensional Cartesian coordinates of the TCR at time t, t is the time step used, and η( t ) is total random motion undergone by the TCR over t. The random motion η( t ) is independent between successive time steps and between dimensions and is normally distributed with
and
for dimensions i, j ∈ {1, 2}, where D is the diffusion coefficient and δ i j is the Kronecker δ function.
We write the biochemical state of the TCR at time t as
where b u is the unbound state, b n is the nonactivating bound state, and b a is the activating bound state. Using the same time step t, we concurrently simulated binding dynamics according to the transition-probability matrix among b u , b n , and b a , given by
using a Monte Carlo simulation. We obtained the transitionprobability matrix from the rate matrix
where the entries in the matrix are in the order given by (6). The concentrations of nonactivating and activating ligand are written as M n and M a , where the overall ligand concentration is set as M and M n and M a are determined by the activating fraction G
M a = GM. (10) We modeled the close contact as a two-dimensional disk with radius R = 220 nm [23,31] centered at the origin,
In keeping with previous KS models of TCR dynamics, we assume that the TCR can only leave the close contact when in the unbound state, so that when the TCR was in the b n or b a biochemical states, the Brownian motion step was repeated until the TCR remained within at the next time. When in the b u state, a Brownian motion step that led to the TCR leaving ended the simulation run. We initialized all simulations with x(0) = [0 nm, 210 nm] and b(0) = b u , in the unbound state. We chose the initial TCR position to be close to the edge of the close contact to reflect recent entry into the close contact, as starting simulations in the unbound state on the close contact boundary would result in the TCR immediately leaving the close contact. We ran each simulation until either the TCR left the close contact or until 10 s of model time had elapsed. The 10 s maximum simulation time is justified by the convergence of P s (t ) well before 10 s for even the slowest activating parameter set (Fig. S11 [24]). In calculating P a (t c ), any value of P s (t ) for t > 10 s was set to P s (10). We repeated each simulation 10 000 times. We performed these simulations using random sampling with the Numpy package in Python [32].
We independently calculated the probability of a TCR remaining in a close contact after time t by adapting a previously published model of TCR diffusion, employing a coupled system of partial differential equations (PDEs) [23].
As with the stochastic simulations, we model the close contact as a disk of fixed radius R = 220 nm [23,31] [as in Eq. ( 11)], which individual TCRs can diffuse in and out of. First, we consider the diffusion of an individual TCR from the instant it has entered the close contact. Starting at this instant, we are interested in calculating the probability that the TCR then exits the close contact after a given amount of time has passed. We write locations within the close contact as x = [x 1 , x 2 ] under the constraint (x 2 1 + x 2 2 ) 1/2 R. A TCR in the close contact can exist in unbound, nonactivating bound, and activating bound form. The twodimensional diffusion coefficients can in principle differ between unbound (D T ) and bound (D C ) TCRs but here are identical. Once a TCR reaches the boundary of the close contact, it is absorbed (leaves the close contact) if unbound to a ligand or is reflected if ligand bound. We write the boundary as
Our goal is to find the probability distributions of unbound, nonactivating bound, and activating bound TCRs over as a function of time since TCR entry to the close contact. We can find these distributions by solving a coupled system of PDEs for the probability density of unbound, T (x, t ), nonactivating bound, C n (x, t ), and activating bound, C a (x, t ), TCR at x at time t
Here, t is the time elapsed from the entry of the TCR to the close contact and k + s (k - ns ) are the on/off rates for nonactivating (activating) ligands binding to TCR.
The initial conditions for this system of PDEs for are
C n (x, 0) = 0 for all x ∈ , (
C a (x, 0) = 0 for all x ∈ , (
where 1 S (x) is the indicator function of the specified set S.
The set A is defined as
with inner radius r i and outer radius r o such that r i < r o < R.
The sets ∂A i and ∂A o are the inner and outer boundaries of A defined as
The boundary conditions are
where n(x) is the outward-facing unit normal vector at x ∈ ∂ . These boundary conditions represent the ability of unbound TCRs to leave the close contact and the inability of bound TCRs to leave the close contact. We used the py-pde Python package [33] to numerically solve the system of PDEs in Eqs. ( 13)-( 15) using finite difference methods.
From a set of replicate trajectories generated from stochastic simulations, we can calculate the probability that over the lifetime of the close contact at least one TCR will reach the fully phosphorylated state. We perform this calculation after the simulations described in Sec. IV B are completed, using the trajectories produced from those simulations. While a TCR is in the close contact, the functional TCR complex can assemble and be phosphorylated. In our model, we ignore complex assembly and only focus on phosphorylation. If the phosphatase concentration is low enough in the close contact, we can assume that phosphorylation occurs through n irreversible steps. Then, the probability that a single TCR complex is fully phosphorylated by time t in n irreversible steps, each with the same rate k p , follows the Erlang distribution
If we were to allow for dephosphorylation, we could write the model as a linear system of ordinary differential equations where the fully phosphorylated state is absorbing. However, in this work we assume that dephosphorylation in the close contact occurs at a slow enough rate to be ignored. In Fig. S9 [24], we also calculate activation probabilities for
as was used in Fernandes et al. [23]. This approach approximates the case of large k p and n in Eq. ( 25), chosen so that the mean of the distribution is 2 s.
To calculate the probability of a single TCR complex reaching the fully phosphorylated state before leaving the close contact for KS-KP, given the fraction of ligands that are activating, we calculate
with
where τ i is the close contact residence time of the TCR in trajectory i. Equation ( 27) is the probability of a single TCR complex reaching the fully phosphorylated state before leaving the close contact, averaged over all trajectories. We calculated the standard deviation of P s (t ) for KS-KP according to
given the independence of τ i between different trajectories.
To calculate the probability of a TCR complex reaching the fully phosphorylated state with CC-KP, we use the same trajectories but instead calculate
where the product is over m j bound periods, indexed by i, for a single trajectory indexed by j. Each term of the product is the probability that in trajectory j during the ith bound period the TCR did not reach the fully phosphorylated state. The length of the ith bound period of trajectory j is
which allows us to consider arbitrary times, including those which fall in the middle of bound period for trajectory. The product is the probability that none of the bound periods in trajectory j resulted in the TCR reaching the fully phosphorylated state, so that one minus this quantity is the probability that at least one of the bound periods resulted in full phosphorylation. This quantity is then averaged over all N trajectories. As with P s (t ) for KS-KP, we calculated the standard deviation of P s (t ) for CC-KP according to σ (P s (t )) = σ 1 -
Equations ( 27) and (30) provide the probability of a single TCR complex reaching the fully phosphorylated state by time t. However, in an actual T cell multiple TCRs can enter and leave a close contact, and we are interested in the probability that any one of those TCRs become fully phosphorylated in the lifetime of the close contact. We can write the probability that no TCR will enter the close contact and become fully phosphorylated over the close contact lifetime by assuming a constant rate of TCR entry into the close contact and by considering m -1 time intervals of length t c /m, for positive integer m. The entry rate of TCR entry into the close contact is [23] κ = 4π DT m ln(A/(π R 2 ) -1)
where T m is the bulk membrane TCR density far from the close contact, A is the cell surface area, and R is the radius of the close contact. Using Eq. ( 33) and discretizing time into intervals of length t c /m, the probability that a TCR will not activate is
1,...,m-1} [1 -P s (t clt c /m)] κt c /m . (34) Taking the limit of P na (t c ; m) as m → ∞, corresponding to continuous time, we have the geometric integral [34] lim m→∞ P na (t c ; m) = exp κ t c 0 ln [1 -P s (t ct )]dt . (35)
Because Eq. ( 35) provides the probability that no TCR will activate in t c , we can then easily find the probability that at least one TCR will activate in t c as
We again note that for t > 10 s in calculation of P a (t c ), we set P s (t ) to P s (10) as justified by convergence of P s (t ) with time (Fig. S11 [24]). Equation (36) is the activation probability shown in Figs. 234. While the notation of our model obscures the dependence of P a (t c ) on G, we emphasize that with a fixed value of M, G fully determines the nonactivating and activating concentrations M n and M a and through the stochastic simulations influences P a (t c ). To find the standard deviation of P a (t c ), we used linear error propagation theory implemented in the Uncertainties Python package [35] using Eqs. (32) and (29) for the standard deviation of P s (t ) for CC-KP and KS-KP, respectively.
Parameters and parameter ranges for our models were taken from the literature [5,23,27,28,31,36]. See Tables S1-S3 for parameter values and sources for specific values [24].
For activation probability P a as a function of G, we define sensitivity as
In practice, we calculate sensitivity as S ≈ P a (0.005) -P a (0) 0.005 -0 .
(38)
We define discrimination as α = 1 -P a (G = 0, t c ) (39) so that if α = 0 the TCR cannot discriminate nonactivating ligands from activating ligands, while if α = 1 the TCR does not activate without activating ligands. Our definitions of sensitivity and discrimination are related to previous definitions but differ in several ways. To make these differences explicit, we contrast our definitions to those in Pettmann et al. [5], who use P 15 (K D ), the concentration of ligand leading to activation of 15% of T cells as a function of ligand dissociation constant in their definitions of discrimination and sensitivity. By fitting a power law to P 15 ,
sensitivity is defined as
log 10 P 15 (1), (41) while discrimination is α = d log 10 P 15 d log 10 K D . (42)
A smaller value of C means higher sensitivity, as 15% activation at lower concentrations means that the TCR is more sensitive to the introduction of small amounts of ligand.
Larger α means more discrimination, since α measures how much an increase in log 10 K D will change log 10 P 15 . Note that if P 15 (K D ) is smaller, sensitivity as we have defined it (S) will be larger. Our definition of sensitivity is similar to that of Pettmann et al., although it is not defined with a specific K D . Rather than measuring how well the TCR can discriminate between ligands with arbitrary K D values as in Eq. (42), our definition measures the extent to which a TCR can mistake a specific nonactivating ligand for a specific activating one.