Transient Kinetic Proofreading

We propose a new stochastic model for understanding the transient kinetic proofreading mechanism in a T-cell. Our model indicates that a stochastic version of absolute ligand discrimination is a consequence of the finite number of receptors on the cell surface; thus, pointing to receptor number control as being critical to T-cell activation. We propose four different metrics to characterize the performance of kinetic proofreading mechanisms. We explore the numerical experiments that explore the trade-offs between speed, specificity, sensitivity, and robustness of T-cell activation as a function of the model parameters. We also consider the impact of receptor clustering on these trade-offs.

T-cells are a component of the adaptive arm of the immune system (Perelson and Weisbuch, 1997). 23 In contrast to the innate immune system, which attacks a broad range of pathogens based on com-24 mon cell-surface markers, each individual cell in the adaptive immune system is specially tailored 25 to sense and react to antigens associated with a specific type of pathogen. Among other things, 26 a T-cell is responsible for sensing and killing pathogens expressing the T-cell's target antigen on 27 their cell surface. T-cells are able to sense very small differences in the binding energy between 28 a specific "non-self" antigen and T-cell receptors (TCRs) and that of a "self" antigen and the TCR to 29 activate only in the presence of the "non-self" antigen. Kinetic proofreading has been proposed 30 as a mechanism that allows T-cells to discriminate with much better accuracy than that implied by 31 the difference in binding energy alone.

32
McKeithan (1995) proposed the first kinetic proofreading-based model for TCR activation. This model is described by the following set of ODEs François and Altan-Bonnet (2016). The existing models of absolute ligand distribution are both mean field and steady state. Therefore, 67 these models are not able to represent the trade-off between the delay and other performance 68 measures. Furthermore, these approaches are also not able to model the impact of stochastic Figure 1. The activation set  bounded by the vertical asymptote at min and the horizontal asymptote at min defining the ideal absolute ligand discrimination activation set  . The delay can be fixed and the probability of activation increased to move  up and to the right, and the probability of activation can be fixed and the delay can be increased to move  down and to the left. the cell activates within delay with probability at least for any value of the affinity .

81
(b) Specificity: The gap 0.95 − 0.05 between the asymptotes along the -direction for = 0.05 82 and = 0.95, i.e. the difference between the affinity values for which the cell almost never 83 activates (i.e. = 0.05) and almost always activates (i.e. = 0.95) in the limit of large density. 84 The affinity of the non-self antigen should be greater than 0.95 , and the affinity of the self-85 antigen should be less than 0.05 to guarantee that the mechanism makes very few errors. 86 Therefore, the smaller the gap, the more specific is the response of the mechanism. 87 (c) Fidelity: This metric quantifies the distance of the activation set  associated with a mech-88 anism, and that associated with ideal absolute ligand discrimination, i.e. the distance of the 89 curve from the two asymptotes. We provide one possible definition of this metric in Sec- ing the mechanism are perturbed. Robustness is an important metric since noise becomes 93 very important at low copy numbers. 94 Although versions of these metrics apply to any proofreading mechanism involved in switching, 95 we trade-off between these desirable properties in the context of a particular model for absolute 96 ligand discrimination. 97 We propose a new stochastic model where a T-cell activates when the total number of antigen-TCR-98 downstream complexes exceeds a certain threshold . We show how to efficiently compute the 99 probability that the T-cell activates with a certain delay , allowing us to investigate the impact of

139
The state of a T-cell in our model is described by two variables: the number of antigen-TCR complexes ≤ and the number of antigen-TCR-Lck complexes ≤ . We let ( , ) denote the state of the T-cell at any time ≥ 0. Since , ≥ 0, ≤ , and the cell activates when = , it follows that the state space is and the dynamics of {( , ) ∶ ≥ 0} is described by a two-dimensional continuous time Markov chain with an absorbing boundary at i.e.  refers to all cell states corresponding to an activated cell. See Figure 2 for a schematic of . 142 We assume that the T-cell must make a decision within a specified delay from the first presentation of the antigen, i.e. if the cell does not have Lck's recruited by time , the cell has missed its opportunity to sense the antigen in time. This is a departure from previous kinetic proofreading models, in which activation is a function of the steady-state concentrations and information theo- retic tools can be used to analyze the system (Cover, 1999;Ganti et al., 2020;Das et al., 2018). We want to characterize the set  of density-affinity pairs ( , ) for which the T-cell activates before with probability ∈ (0, 1). In order to characterize the operating characteristics of the T-cell, we define the first passage time ( , ) from the state ( , ) ∈ ∖ to a state ( , ) ∈  as follows: i.e. ( , ) is the first time at which a cell, initially with ≥ 0 antigen-bound TCRs and 0 ≤ ≤ antigen-TCR-Lck complexes, has antigen-TCR-Lck complexes. We emphasize that ( , ) is a function of the density-affinity pair ( , ) since the rate matrix (see Appendix for details) depends on ( , ). Since the random time to activation after the initial presentation of the antigen is given by In order to characterize  , we need to characterize the entire distribution of (0,0) ( , ) as opposed 143 to just the mean first passage time (MPFT), as was done in earlier works (Husain, 2018;Taylor et al., 144 2017). We are interested in characterizing sets  for a range of in order to understand the impact 145 of the probability threshold .

146
Since an analytical expression for the distribution of the first passage time (0,0) is unknown, we 147 compute and invert the Laplace transform  (0,0) ( ) = [ − (0,0) ] of (0,0) to compute the entire cumu- denote the lowest density for an antigen with affinity that ensures that the cell activates within 154 a delay . In Figure 3 we plot  ( ) vs for ∈ {0.05, 0.5, 0.95}. The set  is above and to the 155 right of the corresponding  -curve. We also plot the minimum density required for the mean 156 first passage time (MFPT) to be below the delay . The plots in this figure also clearly show the additional information in the collection of plots  for different values of . For example, the plots 158 show that the set  is much more sensitive to affinity as compared to the density . This 159 information is not available in the single plot corresponding to the MFPT.

160
The  curve is very closely approximated by two asymptotes: an asymptote in the -direction 161 at the affinity min , and an asymptote in the -direction at the density min . (See Figure 3 for the 162 asymptotes corresponding to the  0.05 ( ) curve.) The thresholds ( min , min ) are functions of the 163 parameters of the mechanism ( , ), the allowed delay , and the probability threshold . This 164 asymptotic behavior directly leads to a probabilistic version of the absolute ligand discrimination 165 described in Lalanne and François (2013). For antigens with sufficiently large affinity, the cell acti-

166
vates with very high probability by the cutoff time , and it does so for all antigen density above 167 some threshold min . For a fixed density above this lower threshold, one can decrease the affinity 168 until a critical affinity min is reached, at which point the probability of activation starts to dramati-169 cally decrease as the affinity is further decreased. Finally, this critical affinity is the same for almost 170 all densities above the density threshold.

171
We show that the asymptotes are an immediate consequence of the fact that is finite, and do 172 not exist if is allowed to be infinite. First consider the asymptote in the -direction. For large 173 enough densities , all of the TCRs are very quickly bound to an antigen, and therefore the Lck 174 dynamics, and hence the activation probability, are independent of density . We rigorously prove 175 this in section C in the appendix.

176
The dynamics that lead to the asymptote in the -direction are more subtle. Justifying the existence  for all points in time. Therefore, the dynamics of the proofreading network in this limit is effectively 198 independent of the affinity . The  curves will also have a horizontal asymptote i.e. a density 199 min below which no antigen will activate the cell in time even with extremely high affinity, and   (4)).
Plots of  ( ) for = 8 = 12 and delay = 10 3 . For each affinity , the minimum density was found using bisection over densities ranging from 10 −4 to 10 16 . The delay = 10 3 seconds was chosen because it represents a realistic amount of time for human immune cells to operate over (10 3 seconds ≈ 15 minutes). Horizontal and vertical asymptotes min and min displayed for = 0.05. Also included is the minimum density needed for a cell exposed to antigen with affinity to have a MFPT less than . Time to compute = 0.50 curve was 30 seconds, time to compute MFPT curve was 2 seconds. This factor of 15 difference comes from using 15 summands to approximate the Bromwich integral. This same ratio holds when testing larger and values (see section A in the appendix).

205
In the Introduction, we introduced four metrics that we argued effectively characterize the perfor-206 mance of a kinetic proofreading mechanism: specificity, sensitivity, fidelity and robustness. In this 207 section, we define these metrics for our proposed mechanism. In order to define these metrics, 208 we need to define the following function: i.e. ( ) is the minimum affinity at which the proofreading network activates with probability by 210 delay when the antigen density is . Note that min = (∞) is the -axis location of the asymptote 211 in the -direction. Next, we define the four metrics. hood of activation decreases quickly as the antigen affinity is decreased. Therefore, networks 217 with higher specificity are better at discriminating between two affinity values.

218
(c) Fidelity: The first two metrics refer to the asymptotes. For the ideal proofreading network, the 219  curve is the -shaped curve defined by the two asymptotes. Fidelity is a metric that mea-220 sures how far  is from the ideal -shape. We need to define some intermediate quantities 221 in order to define fidelity.

222
Fix > 0. We say that the curve  ( ) reaches within a factor (1 + ) of the asymptote in the -direction at = ((1 + ) min ), and reaches within a factor (1 + ) of the asymptote in the -direction at =  ((1 + ) min ). See Figure 4. Monotonicity properties of the model imply that the  ( ) curve lies in the following two -shaped regions: The set  would represent an ideal proof-reading curve provided the maximum of thewidth and -width ( − min ) is close to zero, i.e. = max min , − min ≈ 0. Similarly, define = max − min , min ), and = min{ , } to be the minimum of the dis- tance to ideal performance implied by either of the two covering sets. We define the fidelity ( , , , ) as the inverse of the minimum "distance" over all possible choices for , i.e.
Note that high fidelity implies that the transition is closer to ideal behavior.

223
(d) Robustness: Since the actual number of TCRs + on the cell surface may differ from because of transport noise, the corresponding  curve may no longer be able to discriminate between the self and non-self antigens. We define the robustness ( , , ) of the kinetic proofreading network with parameters ( , , ) to errors in as follows: i.e. ( , , , ) is the maximum allowed relative change in before the kinetic proofread-224 ing network is unable to differentiate between the self and non-self antigens.

225
The robustness in is defined as It is possible that the robustness ( , , , ) is very sensitive to . In order to study this effect, we define robustness ( , , , ) to simultaneous perturbations in and as follows: To summarize, a kinetic proofreading network that exhibits perfect absolute ligand discrimination 226 has zero specificity, zero sensitivity, and infinite fidelity. Next, we investigate the possible values of 227 sensitivity, specificity, and fidelity achievable by the proposed Markov chain model as a function of 228 the parmeters ( , , ).

229
We simulated the performance of the kinetic proofreading system for of points corresponds to values of ( , , ) that are Pareto optimal with respect to sensitivity and 235 delay, i.e. these points represent the optimal trade-off between these two metrics. Note that on 236 the Pareto front, the sensitivity min is a decreasing function of the delay , i.e. the proofreading 237 network can be designed to be more sensitive if a higher delay is tolerable.

238
Next, we want to visualize the ( , , ) Pareto front. There are often several close together ( , ) values that result in very different values for fidelity. In order to account for this, we take the maximum value of fidelity in a neighborhood of ( , ), and define The color of the point ( , ) in figure 5 denotes the value of ( , ). On the Pareto front, the fidelity 239 and sensitivity both increase as the allowed delay increases, and begins to level off between 240 = 100 and = 1000 seconds, a biologically realistic range of activation times. In the right-hand 241 plot in Figure 5 we see that, as one would expect, increasing the required specificity by increasing 242 self limits the delay, sensitivity, and fidelity with which the cell can operate.

243
In Figure 6 the points are colored by robustness. In the first row the points are colored by ro-244 bustness to perturbations in , i.e. defined in (7), in the second row the points are colored by 245 robustness to perturbations in , i.e. defined in (8)  To approximate the impact of clustering we modify the original model as follows. We assume that 260 the cell surface has two compartments, each initially having TCRs. The dynamics of the two com-261 partments are assumed to be independent. We let denote the antigen density on compartment 262 = 1, 2, and say that the cell activates once one of the compartments has recruited Lck molecules. 263 We introduce volatility in the antigen density by setting 1 = (1 − ) and 2 = (1 + ) for some 264 ∼ uniform (−1, 1), so that the average antigen density over the whole surface is still . We assume 265 that the cell is able to sense density and places (1 + ) receptors in the compartment where the 266 density is higher, and (1 − ) receptors in compartment with the lower density. We propose a kinetic proofreading model that exhibits a stochastic version of absolute ligand dis-285 crimination as a direct consequence of controlling the number of receptors. We propose that the 286 performance of a kinetic proofreading mechanism can be characterized by sensitivity, specificity, 287 fidelity, robustness to biochemical errors and delay. We explored the trade-off between these dif-288 ferent metrics within our proposed model. We also examined the impact of receptor clustering on 289 these trade-offs.

290
In Figure 9 we plot the predictions of our model as the number of receptors is changed. Our The full rate matrix for the model described in section 2.1 is defined as follows (Sigman, 2009).

380
The cumulative density function (CDF) is given by where is any real number to the right of all singularities of  (0,0) ( )∕ . Let = + , for , ∈ ℝ.
the Bromwich intergral in (13) with an infinite series Abate et al. (2000) lim Re( (0,0) ( + ℎ))cos( ℎ ), and appropriately truncate the series to approximate the integral. For ℎ and the number of terms 383 in the approximation, we use the specification in Abate et al. (2000). We found that the number of Throughout the text we implicitly assume that the delay is chosen to be large enough so that it's 397 possible for the cell to activate quickly enough with sufficiently high probability. In this section, we 398 compute a lower bound for the delay such that ℙ( (0,0) ( , ) ≤ ) ≥ for some ( , ).

399
Since the forward transitions, i.e. those that increase ( , ), are increasing in and the backward transitions, i.e. those that decrease ( , ), are decreasing in , it follows that the first passage time where Exp( ) denotes an exponential random variable with rate and Hypo( 1 , ..., ) is a hypoexponential distribution defined by the given rates. The same limit results from taking the limit → ∞ in the Markov chain ( ) that we discuss in Section D. Thus, in order to ensure that the network reaches  with probability at least , we need We ensure that ≥ * ( , , ) in all of our analyses so that min (∞) and min (∞) are well-defined. For brevity, let min ( ) denote min ( ; , , , ).

427
The intuition here is that in the limit of → ∞, all TCRs are instantaneously bound to anti-428 gens, and the activation is determined by the dynamics of the antigen bound TCRs acquiring Lck 429 molecules. In Lemma C.2 we rigorously establish the first result.  Proof. The rate out of each state ( , ) in the chain is upper bounded by + + + . We define the following four independent Poisson processes.
We think of the arrivals for the process ( ) as "eastward" movements from ( , ) to ( , +1), arrivals for ( ) as "westward" movements from ( , ) to ( − 1, ), arrivals for ( ) as "northward" movements from ( , ) to ( , + 1), and arrivals for ( ) as "southwestward" movements from ( , ) to ( − 1, − 1). We can see that Next, we have as → ∞, and  Next, we establish that dynamics of the 2D Markov chain in the limit of → ∞ is given by the 1D to an increase in .

448
By construction, defined above has the same distribution as the number of Lcks in the Markov 449 chain. Also by construction, ≤ on a sample-path basis, and therefore ℙ ( ,0) Additionally, − ≤ ( ( )) on a sample-path basis, with ( ) being the total amount of time the Markov chain spends in states with < for ∈ (0, ], and ( ( )) being the number of arrivals in during the periods of time in (0, ] when < . Thus, we have that where the second-to-last inequality uses the fact that ( ) is independent of . Taking limits, we see

467
To study the behavior of min ( ) as → ∞, we define where the last equality holds because the limit and probability can be exchanged for the same reason as in section C. In the limit of → ∞, the unbinding rate → 0, and therefore the rates confirms that min (∞) = ( ) min and that ( ) min can be computed much faster than min (∞) (∼10-100 times faster for moderate-sized and ). Computing the horizontal asymptote is faster for smaller 477 and because it only requires one pass through the problem data (of size ( )), whereas 478 computing the vertical asymptote requires two passes through the problem data (of size ( )).

479
The difference in complexity starts to show up once and/or become larger. In Sections C and D we established that the vertical and horizontal asymptotes of the ( ) curve 484 can be characterized in terms of simpler Markov chains whose first passage time distributions can 485 be computed much more quickly. In this section, we study the shape of the  curve for finite 486 values of and . We do this because being able to efficiently calculate ( ) and  ( ) is essential 487 for computing the fidelity metric discussed in section 2.3.

488
In order to efficiently simulate the entire  curve and investigate its "distance" from ideal absolute The Markov chain ( ) starts from the initial state ( ) (0) defined below. Before describing the 495 procedure for selecting ( ) (0), we describe the rationale for this decoupling of the ( , ) dynamics.

496
The T-cell activates when = . We would be able to efficiently estimate the quantities of interest if one could describe a Markov chain where the transition rates are only a function of the state . The rate from state to state − 1 in the 2-D Markov chain is given by , i.e. it only depends on , and so we keep this rate as it is. However, the rate from state to + 1 depends on the number of TCR-antigen complexes , and therefore cannot be written as a function of only the state . We know that ≥ , and so define a Markov chain ( ) with the state space consisting of { , … , } and an absorbing state  that represents the event → + 1. Note that in the 2-D Markov chain, it's possible that is much larger than when the chain makes a transition from → + 1, and so setting ( ) (0) = could significantly overestimate the first passage time. To correct this, we set the initial states ( ) (0) in the following recursive manner:  (17).

505
Inverting the Laplace transform for the birth-death process allows us to approximate and .

506
Since calculating the mean first passage time to  for each of the ( ) processes involves solving a 507 tri-diagonal system of equations, this step takes at most ( ) time, resulting in a total complexity 508 of ( ). The Laplace transform for the birth-death process can be computed in ( ) time. Thus, 509 the approximate model results in an ( ) speedup compared to using the model .