Multiple steady states and the form of response functions to antigen in a model for the initiation of T cell activation

The aim of this paper is to study the qualitative behaviour predicted by a mathematical model for the initial stage of T cell activation. The state variables in the model are the concentrations of phosphorylation states of the T cell receptor complex and the phosphatase SHP-1 in the cell. It is shown that these quantities cannot approach zero and that the model possesses more than one positive steady state for certain values of the parameters. It can also exhibit damped oscillations. It is proved that the chemical concentration which represents the degree of activation of the cell, that of the maximally phosphorylated form of the T cell receptor complex, is in general a non-monotone function of the activating signal. In particular there are cases where there is a value of the dissociation constant of the ligand from the receptor which produces an optimal activation of the T cell. In this way the results of certain simulations in the literature have been confirmed rigorously and some important features which had not previously been seen have been discovered.


Introduction
In humans and other vertebrates, the immune system is of crucial importance for protecting an individual from dangers such as pathogens, toxins and cancer.  shown that when N = 3, there are parameters for which three positive steady states exist (theorem 3.1). A numerical calculation reveals that for a specific choice of these parameters two of the steady states are stable, while the third is a saddle. For N ≤ 2, there is a unique steady state and in the case N = 1 it is proved to be globally asymptotically stable. There are parameter values for which the approach to this steady state is oscillatory.
The qualitative behaviour of the steady-state concentration of the maximally phosphorylated state, which expresses the degree of activation of the T cell, as a function of the antigen concentration and the dissociation time, is investigated in the case where only the agonist is present in §4. Let us consider the function f (L 1 , ν 1 ), which expresses the degree of activation in terms of the parameters L 1 (concentration of agonist ligand) and ν 1 (reaction rate for the dissociation of the ligand from the receptor, i.e. the reciprocal of the dissociation time). It is shown that the dependence exhibits certain types of non-monotone behaviour in some cases. The results obtained include both rigorous results on general features of the function f (theorem 4.2) and simulations which reveal more detailed features. In particular, it is found that are values of the parameters in the model for which the function f has a maximum as a function of ν 1 for fixed L 1 . In other words, there is a value of the dissociation time which is optimal for T-cell activation. Thus, the model studied here is able to reproduce this fact which has been experimentally observed [6].
The analysis of the response function is extended to cover the effects of the antagonist in §5. The last section is devoted to conclusions and an outlook.

Definition of the model
In the introduction, it was stated that a T cell recognizes an antigen. In more detail, the molecule concerned is a peptide (a small protein) which is bound to a host molecule called a major histocompatibility complex (MHC) molecule. Thus, we talk about a pMHC complex as the object to be recognized. In the model of [4], two types of pMHC complexes are considered. The first, called an agonist, represents the case where the antigen comes from a pathogen and should activate the T cell. The second, called an antagonist, represents the case of a self-antigen, which should not activate the T cell. Detection takes place through the binding of a pMHC complex to the TCR. As explained in the introduction, when this happens certain proteins associated with the TCR are phosphorylated, i.e. phosphate groups become attached to them. For simplicity, we describe this by saying that the receptor-pMHC complex is phosphorylated.
The reaction network for the model of [4] is shown in figure 1. The state variables will now be listed. The concentration of unphosphorylated complexes of the TCR with the agonist is denoted by C 0 and the concentration of unphosphorylated complexes of the TCR with the antagonist is denoted by D 0 . C j and D j are the corresponding quantities for the case of j phosphorylations, up to a maximum value N. The specific value of N has little influence in what follows. In some of our results, we choose N small so as to obtain the simplest possible mathematical setting. The number of phosphorylation sites relevant to the models of [3,4] have been discussed in the introduction. R, L 1 and L 2 are the total concentrations of receptors and the two ligands, i.e. the agonist and antagonist. Another important element of the system is SHP-1. This substance is a phosphatase which means that when active it can remove phosphate groups from the receptor-pMHC complex. It contributes a negative feedback loop to the system. S is the concentration of active SHP-1. The receptor complexes are subject to phosphorylation with rate constant φ and dephosphorylation with rate constant b. They are also dephosphorylated by SHP-1 with rate constant γ and dissociate with rate constants ν 1 and ν 2 . Antigens bind to the receptor with rate constant κ. SHP-1 is activated by the singly phosphorylated complexes with rate constant α and deactivated with rate constant β. All the rate constants are assumed positive. S T is the total concentration of SHP-1. It is assumed that all reactions exhibit mass action kinetics and this leads to the following system of equations: 3)  Figure 1. The model considered in this paper. The species R is the T-cell receptor (TCR), and L 1 and L 2 are the two ligands, i.e. the agonist and antagonist. The species C 0 are unphosphorylated complexes of the TCR with the agonist, and the C j 's are the j-phosphorylated complexes. The D j 's are the analogous complexes for the antagonist. The phosphatase SHP-1 provides a negative feedback, and is represented by S. The different reactions represent receptor complex phosphorylation with rate constant φ and dephosphorylation with rate constant b, as well as receptor complex dephosphorylation by S with rate constant γ and dissociation rate constants ν 1 and ν 2 . Antigens bind to R with rate constant κ, and S is activated by the singly phosphorylated complexes with rate constant α and deactivated with rate constant β.Ḋ In a direct formulation of the system as arising from the reaction network, it is necessary to include the concentrations of free ligands, free receptors and inactive phosphatase. This extended system has four conservation laws corresponding to the total amounts of ligands, receptors and phosphatase. The explicit form of the conserved quantities is where L 1,U , L 2,U and R U are the concentrations of unbound ligands and receptors and S I is the concentration of the inactive form of SHP-1. Using these conservation laws to eliminate the additional variables leads to the system (2.1)-(2.7). The right-hand sides of the equations are Lipschitz and so there is a unique solution corresponding to each choice of initial data. To have a biologically relevant solution, the quantities in the extended system should be non-negative. It is a well-known fact for reaction networks of this type that data for which all concentrations are positive give rise to solutions with the same property and that data for which all concentrations are non-negative give rise to non-negative solutions. In terms of (2.1)-(2.7), this implies statements about the positivity of the quantities S, C j and D j and of the differences S T − S, Let us call the region where all these quantities are strictly positive the biologically feasible region. Note that owing to the conservation laws, this region is bounded. Let Σ 1 = N j=0 C j and Σ 2 = N j=0 D j . Then it follows from (2.1) to (2.7) thaṫ (2.9) Lemma 2.1. Consider a solution (S(t), C 0 (t), . . . , C N (t), D 0 (t), . . . , D N (t)) in the closure of the biologically feasible region. Then if (S * , C * 0 , . . . , C * N , D * 0 , . . . , D * N ) is an ω-limit point of this solution it is also in the biologically feasible region. In particular, any steady state is in the biologically feasible region.
Proof. If S * = S T we can consider the solution starting at that point at some time t 0 . Since the ω-limit set of a given solution is invariant, the solution under consideration lies entirely in the ω-limit set of the original solution. In particular, it is contained in the closure of the biologically feasible region. The solution starting at the point with S * = S T satisfiesṠ(t 0 ) < 0 because the first term on the right-hand side of (2.1) is zero for t = t 0 and the second term negative. Hence, the solution starting at the ω-limit point satisfies the inequality S(t) > S T for t slightly less than t 0 , a contradiction to the fact that the original solution was in the biologically feasible region. In a similar way, equation (2.8) implies that N j=0 C * j cannot attain the value L 1 and equation (2.9) implies that N j=0 D * j cannot attain the value L 2 . Summing (2.8) and (2.9) shows that N j=0 C * j + N j=0 D * j cannot attain the value R. Note next that C 0 cannot be zero at an ω-limit point. For if it is zero at such a point, we can consider the solution passing through that point at a time t 0 . As the inequalities already proved imply that the first term in equation (2.2) is positive for t = t 0 that equation implies thatĊ 0 (t 0 ) > 0 and that C 0 (t) < 0 for t slightly less than t 0 , a contradiction. Once the positivity of C 0 has been proved we can use equation (2.3) with j = 1 to show that C 1 cannot be zero at an ω-limit point. This, in turn, allows us to prove using (2.1) that S can never be zero at an ω-limit point. In a similar way, it can be concluded successively that C 2 , . . . , C N and D 0 , . . . , D N are positive at any ω-limit point of a non-negative solution. This concludes the proof of the lemma.
The fact that all ω-limit points of solutions in the closure of the biologically feasible region are in the biologically feasible region, together with the fact that the closure of that region is compact, implies that the infimum of the distance of a given solution to the boundary in the limit t → ∞ is strictly positive. When this last property holds, the system is said to be persistent [7]. Note in addition that the closure of the biologically feasible region is convex and hence homeomorphic to a closed ball in a Euclidean space. It follows from the Brouwer fixed point theorem that a steady state exists (cf. [8], ch. I, theorem 8.2). As steady states on the boundary have already been excluded, we can conclude that there is at least one steady state in the biologically feasible region for any fixed choice of parameters. That this is the case was assumed implicitly in [4].

Multiplicity of steady states
A question not addressed in [4] is whether there might exist more than one positive steady state for a fixed choice of parameters. In this section, it is shown that for some values of N and the reaction constants this is the case. The aim is to find any parameter values with this property while not worrying for the moment how biologically relevant this choice of parameters is. Let f 1 and f 2 denote the right-hand sides of equations (2.8) and (2.9). Then ∂f 1 /∂Σ 2 and ∂f 2 /∂Σ 1 are negative and hence the system (2.8)-(2.9) is competitive. It follows that every solution of this system converges to a steady state as t → ∞ [9].
A steady state (Σ * 1 , Σ * 2 ) of (2.8)-(2.9) satisfies the equations and We can solve for Σ * 1 and Σ * 2 as functions of Σ * 1 + Σ * 2 : and The function of Σ * 1 + Σ * 2 on the left-hand side of this equation is decreasing on the interval [0, L 1 + L 2 ]. The function on the right-hand side is increasing on the interval [0, R]. Their graphs can intersect in at most one point. We already know that they must intersect since a positive steady state of the full system exists. That they intersect can also be seen directly. For in all cases, the left-hand side is greater than the right-hand side for Σ * 1 + Σ * 2 = 0 and the opposite inequality holds for and [0, min{L 2 , R}], respectively. The quantities Σ * 1 and Σ * 2 are functions of the parameters R, L 1 , L 2 , κ, ν 1 and ν 2 .
It can be concluded that the solution passing through an ω-limit point of a solution of the original system satisfies a simplified system containing Σ * 1 and Σ * 2 as parameters. C 0 and D 0 can be eliminated from this system in favour of the other C j and D j . The result iṡ This form of the equations is valid for N ≥ 3. In the case N = 2, it is still correct if it is taken into account that the condition 2 ≤ j ≤ N − 1 is never satisfied so that the equations containing that condition are absent. The sum from j = 3 to N is zero in that case. The case N = 1 is exceptional from the point of the notation.
To get more information, we restrict in the remainder of this section to, what we call, the agonist-only case. This is obtained from the system (2.1)-(2.7) by setting L 2 and the D i to zero. There is a corresponding limiting system, which is obtained from (3.6) to (3.12) by setting Σ * 2 and the D i to zero. In this case, we write Σ * instead of Σ * 1 for brevity. Consider the limiting system in the agonist-only case with N = 1. (3.14) Solving the equationṠ = 0 for C 1 and substituting the result into the equationĊ 1 = 0 gives the quadratic equation As the quadratic polynomial has positive leading term and is negative for S = 0, it is clear that it has a unique positive root. It follows from (3.15) that this root is less than S T . Equation (3.14) implies that C 1 < Σ * at a steady state and so these quantities can be completed to a steady state of the original system by defining C 0 = Σ * − C 1 . The steady state is unique in this case. In the case N = 2, the equations arė Proceeding in a manner analogous to what we did in the case N = 1 it is possible to get a cubic equation for S in the case N = 2, which we can write schematically in the form p(S) = N k=0 a k S k . We have The sequence of signs of the coefficients a i is either (−, −, +, +) or (−, +, +, +). There is precisely one change of sign and thus by Descartes' rule of signs the polynomial has precisely one positive root. Once a value of S is given, the values of C 1 and C 2 at the steady state can be determined successively. Following the arguments in the case N = 1, we see that S < S T , C 1 + C 2 < Σ * and that the steady state is unique. In the case N = 3, the system iṡ A calculation for N = 3 analogous to those already done gives a quartic polynomial. Its coefficients are The coefficient a 0 is negative, while a 3 and a 4 are positive. Unless a 1 > 0 and a 2 < 0 Descartes' rule of signs implies that the polynomial only has one positive root. Otherwise, the rule implies that it has one or three positive roots (counting multiplicity), but does not decide between these two cases.
It will now be shown that in the case N = 3, there are values of the coefficients for which the polynomial p(S) has three positive roots. To do this, we vary the coefficients S T and ν 1 in the system (3.19)-(3.22) and keep all others fixed. Note that these coefficients come from the parameters in the agonist-only case of (2.1)-(2.4). To obtain the desired variation of the coefficients, we fix all parameters in (2.1)-(2.4) except S T , ν 1 and κ and vary κ in such a way that ν 1 /κ does not change. This ensures that Σ * does not change. In fact, we may simplify the calculations by setting b = 0 because if three positive roots can be obtained in that case the same thing can be obtained for b small and positive by continuity. Suppose that S T and ν 1 depend on a parameter with both of them being positive for > 0. Suppose in addition that in the limit → 0, we have the asymptotic relations S T =S T −1 + o( −1 ) and for constantsS T andν 1 . Then we obtain asymptotic expansions a 4 = A 4 , a 3 = A 3 + o(1), Then q(1) converges to A 2 for → 0 and is thus negative for small enough. The same is true for p (1). On the other hand, Hence for sufficiently small p( 2 ) > 0. Putting these facts together shows that p has three positive roots when is small. For each of these roots, the values of C 1 , C 2 and C 3 at the steady state can be determined successively. S < S T , gives a steady state of the original system. It has already been noted that p cannot have more than three positive roots. There are parameter values for which the positive steady state is unique. To see this, it is enough to assume that S T is small while keeping the other parameters fixed. Then a i > 0 for all i > 0 and the polynomial can have no more that one positive root because its derivative has no positive root. These results can be summed up as follows: Theorem 3.1. The agonist-only case of the system (2.1)-(2.7) has exactly one positive steady state for N = 1 and N = 2. In the case N = 3, there are parameters for which it has three positive steady states and it can never have more than three.
A concrete example of parameters for which there are three positive steady states is obtained by setting α, β, γ , φ, L 1 and R equal to one and choosing S T = 10, κ = 2 × 10 −4 , ν 1 = 10 −4 . A computer calculation shows that the coordinates (S * , It shows in addition that while the first and second of these steady states are asymptotically stable the third is a saddle with a one-dimensional unstable manifold. A plot of the steady states as a function of the parameter L 1 (figure 2), suggests that there is a fold bifurcation.
For higher values of N it is possible to derive a polynomial equation of degree N + 1 for S. There is no obvious reason why this polynomial should not have an arbitrarily large number of positive roots for N arbitrarily large. A simple upper bound is that the polynomial can have no more than N positive roots for N odd and no more than N + 1 for N even. Simulations indicate that in the case N = 5 there are parameters for which three steady states exist but no parameters were found for which there are more than three for any value of N.
In general, it is difficult to obtain information about the stability of the steady states by analytic methods. In the case N = 1, the vector field defining the dynamical system has negative divergence and so by Dulac's criterion und Poincaré-Bendixson theory, all solutions converge to the steady state as t → ∞. The system can exhibit damped oscillations as will now be shown. To do this, we choose parameters so that For fixed values of the quantities R and S T , the quantities C 1 and S are bounded uniformly in the quantities appearing in (3.27). Thus, if we make α and β small while fixing the other parameters, we can arrange that the left-hand side is smaller than the right-hand side. If starting from there, we make β large while fixing the other parameters we can arrange that the left-hand side of (3.27) is greater than that of the right-hand side. It follows that parameter values exist for which (3.27) holds. Why this is interesting is that the discriminant of the characteristic equation of the linearization is the sum of a term which vanishes when (3.27) holds and the expression −4αγ (S T − S)C 1 . Thus when (3.27) holds, the linearization has eigenvalues with negative real part and non-zero imaginary part and there are damped oscillations. An interesting limiting case of the agonist-only system is obtained by assuming that α = 0 and S = 0. We refer to this as the kinetic proofreading system because it is closely related to McKeithan's kinetic proofreading model [10]. In fact, McKeithan only considered the case b = 0, but this makes no essential difference for the analysis which follows. It was observed by Sontag [11] that the deficiency zero theorem of chemical reaction network theory can be applied to McKeithan's system to conclude that there is a unique steady state in each stoichiometric compatibility class and that this solution is asymptotically stable in its class. Strictly speaking, chemical reaction network theory is applied to the extended system which includes free receptors and free ligand as variables. To show that the steady state is globally asymptotically stable, it suffices to show that there are no ω-limit points on the boundary. That this is the case can be proved just as we did for the full system above. The steady state is hyperbolic as follows from appendix C of [12].
Consider now the full agonist-only system. Setting α = 0 gives a system where the kinetic proofreading system is coupled to a system describing the decay of S. The steady state of the kinetic proofreading system gives rise to a steady state of the agonist-only system with α = 0 which is on the boundary of the biologically feasible region and is a hyperbolic sink. Denote its coordinates by (0, C * j ). For α small and positive, there exists a hyperbolic sink which is a small perturbation of that for α = 0. It must be in the biologically feasible region because C 1 > 0 there and equation (2.1) would imply thaṫ S > 0 there if S were negative. Thus for sufficiently small values of α, there exists a positive steady state which is a hyperbolic sink (S * (α), C * j (α)) close to (0, C * j ). There exists a positive number r such that for α sufficiently small, say α ≤ α 0 , (S * (α), C * j (α)) is the only ω-limit point of any solution in the open ball of radius r about that steady state.
Let h(C j ) be the Lyapunov function in the proof of the deficiency zero theorem. It is known from general arguments thatḣ ≤ 0 with equality only for C j = C * j . It follows that on the complement of the ball of radius r about the steady state the functionḣ has a strictly negative maximum. We can consider the behaviour of the functionḣ for solutions of the system for positive α. For small α, it is still a Lyapunov function on the complement of a small ball about the steady state, while there are no ω-limit points except the steady state itself within that ball. Hence for α, sufficiently small a solution can have no ω-limit points other than the steady state. It follows that for α small the steady state is globally asymptotically stable. Of course, this means that the limiting system obtained from the agonist-only system by passing to a solution through an ω-limit point also has a unique steady state which is globally asymptotically stable for α sufficiently small. A similar argument applies in the case of the full system (2.1)-(2.7) because in that case the system obtained by setting α and S to zero is just the product of two copies of the corresponding system in the agonist-only case.

The response function
This section is concerned with the agonist-only system. From a biological point of view, the essential input parameters to the system are the ligand concentration L 1 and the binding time of the ligand to the receptor, which in the model corresponds to ν −1 1 . The latter is a measure of the signal strength. The essential output is the value of C N which is a measure of the activation of the T cell. Given values of L 1 , ν 1 and the other parameters, we can consider the value of C N in a steady state. In fact, it is more convenient to use the quantities log C N and log L 1 . This leads to a response function log C * N = F(log L 1 , ν 1 ). If there is more than one steady state for a given choice of the parameters, this has to be thought of as a multi-valued function. It might naively be thought that F should be an increasing function of L 1 and a decreasing function of ν 1 : more antigen leads to more activation of the T cell and a longer binding time leads to more activation. This turns out not to be the case and the function F is not a monotone function of its arguments. This was observed in the case of the dependence on L 1 in the simulations of [4]. It is possible to understand intuitively how this situation can arise. An increase in the stimulation of the T cell leads to activation of SHP-1 and that in turn has a negative effect on the activation of the T cell. Many of the calculations in this section are guided by those in [4].
The behaviour of the response function will be estimated in various parameter ranges. To do this, it is useful to parametrize the solutions in a certain manner which will now be described. In the case of a steady state, the equation (2.3) is a linear difference equation for the C j with constant coefficients. This suggests looking for power-law solutions, an idea which motivates the following result.

Lemma 4.1. Steady-state solutions of equations (2.2)-(2.4) in the agonist-only case can be parametrized in the form
where the coefficients r ± and a ± are positive and depend on S. The quantities r + and r − are given by and satisfy r − < 1 < r + .
Proof. Note first that the quantities r ± in (4.2) are the roots of the characteristic equation coefficients a − and a + together with the equations obtained by substituting (4.1) into the equationsĊ 0 = 0 andĊ N = 0. The explicit form of these last equations is (4.4) and It follows from the discussion in §3 that N j=0 C j , which was denoted there by Σ * 1 , is uniquely determined for fixed values of the parameters in (2.2)-(2.4) and fixed S. Thus for fixed values of these parameters and S, all quantities in (4.4) and (4.5) except a − and a + are known. It will now be shown that these equations have a unique solution for a − and a + . Note that as can most easily be seen by multiplying out the left-hand side of this equation and substituting for r + r − and r + + r − , which are the sum and product of the roots of the characteristic equation (4.3). Thus equation (4.5) gives a positive expression for a + /a − . Note also that (4.6) implies that the factors in the product on the left-hand side of that equation have opposite signs. As r − < r + , the first factor is positive and the second negative. Substituting the expression for a + /a − into (4.4) gives an equation of the form whose right-hand side is positive. Here and It follows from the fact that the first factor on the left-hand side of (4.6) is positive that the first factor in the expression for A is negative and hence that A itself is positive. In addition, a straightforward computation shows that A > B. If B were not positive, then the quantity in square brackets on the lefthand side of (4.7) would be positive. If B is positive, then the fact that r − < r + implies that the quantity in square brackets is again positive. Hence in any case, (4.7) can be solved to give a unique positive value of a − . Then a + can be determined in such a way that (4.4) and (4.5) hold. This completes the proof of lemma 4.1.
Lemma 4.1 shows that for fixed parameters in (2.2)-(2.4) and a fixed value of S the steady-state values of all the C j are determined, but this does not yet give expressions for the C j which can be directly applied to study the properties of the response function. For the purposes of what follows, it is convenient to rewrite (2.8) in the form (4.10) The equation for S can be solved to give the relation S = S T (C 1 /(C 1 + C * ) with C * = β/α. Summing the expression for C j given in Lemma 2 over j gives N j=0 C j = a + r N+1 The following equation relating a − and a + is equation (21) of [4]: Combining the last two equations gives Having completed the necessary preliminaries, we now proceed to study the qualitative behaviour of the response function in different regimes. When L 1 is small, it is to be expected that the concentration of the phosphatase is small and that the response function resembles that of the kinetic proofreading model. It will now be shown that when L 1 is small, the leading term in the function F depends linearly on log L 1 with slope one and the additive constant in this linear function will be determined. The equation (4.10) can be written in the form (4.14) Note that N j=0 C j ≤ L 1 so that this equation implies that where − 1 4 < q < 0. Using (4.12), it is possible to write down an explicit expression for C N , namely It follows from (4.13) that (4.17) Combining these equations gives The function of r − and r + in this equation defines a function of S. This function of S tends to a positive limiting value as S → 0. Now C 1 ≤ N j=0 C i = O(L 1 ) and S = O(C 1 ). Hence for R fixed, we can replace the function of r + and r − in the above expression by its limiting value for S → 0. If the resulting relation is plotted logarithmically, it gives a straight line of slope one as the leading order approximation in the limit log L 1 → −∞.
Next, we will derive a lower bound for γ S in the case that S T is large. This will be proved by contradiction. Suppose that γ S ≤ ρ for some ρ > 0. Then it follows from the characteristic equation that r − ≥ φ/(φ + ρ + ν 1 ). Using this in the equation for C 1 gives C 1 ≥ (φν 1 /(φ + ν 1 )(φ + ρ + ν 1 ))( N j=0 C j ). It follows that It is then clear that for a given value of ρ and fixed values of the parameters other than S T , this leads to a contradiction if S T is sufficiently large. In other words, given any ρ > 0 there is a lower bound for S T which implies that γ S ≥ ρ. It is convenient to make the restrictions that κR ≥ 1 and L 1 /R ≤ 1 since then it is possible to replace N j=0 C j in (4.19) by 3L 1 /4(1 + ν 1 ) using (4.15). From (4.2), it can be concluded that and where η = (φ + ν 1 )/(b + γ S). This gives approximate expressions for the roots of the characteristic equation if (φ + ν 1 )/(b + γ S) is small. As a consequence of these equations Taking the expression for C 1 supplied by Lemma 2 and using (4.12), (4.13) and (4.15) gives This implies that C 1 = O(η) and the expression relating S and C 1 then shows that S S T = O(η). In fact, These relations indicate that in leading order r − is proportional to S. However, it is also the case that which indicates that in leading order r − is proportional to S −1 . Hence and where η = max{η, b/(γ S)}. Combining these two relations gives Substituting this back into the equation for r − gives This means that where η = max{η , L 1 /R}. Choosing L 1 small enough makes L 1 /R small. With L 1 fixed, making S T large enough makes η small. Thus, η can be made as small as desired by choosing L 1 sufficiently small and S T sufficiently large.    Proof. To obtain the conclusion of the theorem, it suffices to show that under the given assumptions η can be made as small as desired. That this is possible follows from the above discussion.
Note that this theorem implies, in particular, that for N > 2 and suitable values of L 1 and S T there exists a range of L 1 in which the response function is decreasing. The theorem also implies that in this regime, the response function can be an increasing function of ν 1 . This effect was not captured by the calculations of [4] because there ν 1 /κR was assumed to be so small as to be negligible.
Finally, we examine the regime where L 1 /R is small, but the phosphatase is close to being completely activated. This means that S/S T is close to one. This holds provided C 1 is sufficiently large compared to C * . It remains to check that such a regime actually occurs for some values of the parameters. It is possible to make N j=0 C j large while keeping L 1 /R constant. This can be done by making R large. This makes a − large without making r − small. Hence it makes C 1 large and hence S close to S T . In this regime, the function of r + and r − occurring in the expression for C N can be replaced by its limit for S → S T and we again get a region where the slope of the graph of log C N as a function log L 1 is one but the line has been shifted compared to that obtained for L 1 /R small.
In [4], these types of behaviour were exhibited numerically in the case N = 5 with biologically reasonable choices of the parameters. We found that changing these parameters a little allows similar observations to be made in the case N = 3. In the plot shown in figure 3, the three regimes can be seen together with a fourth regime where L 1 /R is no longer small. It is clear that a regime of this type must exist because the response function is globally bounded.
We now turn to the dependence of the response function on ν 1 . It has been suggested in [13] that the kinetic proofreading model with negative feedback as studied here is not able to explain the presence of an optimal dissociation time, a biological effect confirmed by the experimental work of [6]. The plots of the response as a function of the dissociation time in that type of model in [13] show that it is increasing. Having an optimal dissociation time would require that there be a region where this function is decreasing. The response function being increasing as a function of the dissociation time corresponds to its being decreasing a function of ν 1 . Here, we have given an analytical proof in theorem 4.2 that there exist parameters for which the response is an increasing function of ν 1 , in contrast with the plots in [13]. As the theorem is of limited help in finding explicit parameters for which this happens, we also did a numerical search and identified parameters of this type. The results are displayed in figure 4, where it is seen that F has a maximum as a function of ν 1 for fixed L 1 , which corresponds to an optimal dissociation time. The conclusion of both the analytical and the numerical work is as follows. The claim that the kinetic proofreading model with feedback can only produce a response which is a decreasing function of the parameter ν 1 is dependent on the parameter values chosen to do the simulations and not a general property of the model. This means that the model of [4] can reproduce the observation of an optimal dissociation time and that as a consequence that phenomenon could arise by the mechanisms taking place in the first few minutes of activation which are included in the model of [4].

Including the antagonist
When the antagonist is included, the output variable expressing the degree of activation of the T cell is C N + D N . Now asymptotic expressions for this quantity will be derived. It has already been shown that for a steady state of the system (2.1)-(2.7), the quantities N j=0 C j and N j=0 D j can be expressed in terms of the parameters. The equation for S can be solved to give the relation S = S T ((C 1 + D 1 )/(C 1 + D 1 + C * )). C j solves the same difference equation as in the agonist-only case and D j solves the difference equation obtained from that one by replacing ν 1 by ν 2 . The quantities r − , r + , a − and a + differ in the two cases. We can, nevertheless, proceed as in the former case to see that the solutions for C j and D j allow parametrizations in terms of these quantities as before. Note that using the equations (2.8) and (2.9), it is possible to eliminate the D j from the equation for C 0 and the C j from the equation for D 0 . Thus, we have coupled equations for the C j and D j which can be analysed just as in the agonist-only case to express C 1 and D 1 as functions of S and the parameters. We can also write C N and D N as functions of Σ 1 and Σ 2 , respectively. Proceeding as in the agonist-only case, we get an expression for C N + D N in the kinetic proofreading regime. The multiple of L 1 obtained there as leading term is replaced by a linear combination of L 1 and L 2 .
Next the intermediate regime will be considered. For this, it is necessary to define a new parameter η = max{(φ + ν 1 )/(b + γ S), (φ + ν 1 )/(b + γ S)}. There are asymptotic expressions for r − and r + where the leading terms are just as in the agonist-only case. In particular, they are the same for C j and D j . Two asymptotic expressions for the quantity C 1 + D 1 can be obtained: This gives an expression for r − in terms of S. As in the agonist-only case, this gives an expression for r − where the dependence on S has been eliminated in leading order: where η is defined in terms of η as in the agonist-only case. Following the steps used in the agonist-only case leads to an expression for C N + D N which is the same as that previously obtained for C N except that the expression κRL 1 /(κR + ν 1 ) is replaced by κRL 1 /(κR + ν 1 ) + κRL 2 /(κR + ν 2 ). This leads in the end to an asymptotic expression for C N + D N under a suitable assumption on L 1 and L 2 . The assumption made in the agonist-only case can naturally be written as an assumption on κRL 1 /(κR + ν 1 ) and in the present case it is replaced by an assumption on κRL 1 /(κR + ν 1 ) + κRL 2 /(κR + ν 2 ). This implies that under certain circumstances, C N + D N increases when L 2 increases and L 1 is held fixed. An increase in the amount of self-antigen can lead to a decrease in the response to a foreign antigen. Note that the restriction needed to make this result hold is that first L 1 and L 2 are sufficiently small and then, with upper limits for these quantities having been fixed, that second S T is sufficiently large. It follows that these conditions can be achieved in situations where L 1 /R and L 2 /R are as small as desired and hence the competition of the antagonist with the agonist for occupancy of the receptor is negligible. Hence the effect by which more antagonist leads to a decrease in the response to an agonist is, in general, owing to the influence of SHP-1. This gives a rigorous confirmation of a fact already observed in [4].

Conclusion and outlook
In this paper, some properties of the solutions of the model of [4] for T-cell activation were proved.
A new discovery was that already in the case of three phosphorylation sites (N = 3), there can exist more than one positive steady state for given values of the parameters. Another new observation is that damped oscillations can occur. It was also proved that, as suggested by the calculations in [4], the output variable C N (concentration of the maximally phosphorylated receptor) is a decreasing function of the concentration L 1 of antigen in some parts of parameter space. In an analogous way, it was proved that under some circumstances the activation in response to an agonist can be decreased by increasing the concentration of the antagonist. It was proved that it can also happen that C N is an increasing function of the dissociation constant ν 1 . This abstract result was given a concrete illustration by a plot showing that C N can have a local maximum as a function of ν 1 . The stability of the steady states was only determined analytically in the very special cases N = 1 and α close to zero. For N = 3, numerical calculations showed the occurrence of two stable steady states for certain values of the parameters. It was proved that damped oscillations occur, but can there also be sustained oscillations (periodic solutions)? It is, thus, clear that there remain several aspects of the dynamics of this system which would profit from further investigations, analytical and numerical.
In immunology, it is important to describe diverse situations including the course of different types of infectious disease, the development of autoimmune diseases and the destruction of tumour cells by the immune system. It would be unreasonable to expect that a simple mechanism could be the key to describing all these situations. One strategy to try to obtain more understanding is to choose one mechanism and to investigate which types of situations it suffices to describe. This may be done by combining mathematical models with experimental data. What are the restrictions under which the type of model studied in this paper might be appropriate? The first assumption is that in the situation to be explained the distinction between self and non-self takes place within an individual T cell. In other words, it is assumed that it is not necessary to consider the population dynamics of the T cells involved or even the interaction of their population with that of other types of immune cells such as regulatory T cells or dendritic cells. A quite different type of mathematical model, where population effects are considered, can be found in [14]. In that case, in contrast with the lifetime dogma, the response depends on the rate of change of the antigen concentration. The second assumption which is important for the models studied here is that the distinction between self and non-self takes place on a sufficiently short time scale, say three minutes. On longer time scales, there may be essential effects related to the spatial distribution of molecules on the T cell surface (formation of the immunological synapse) so that a description by means of ordinary differential equations may be insufficient. It may also happen that some TCRs become inactive on a longer time scale (limiting signalling model, cf. [6]).
In this paper, we have concentrated on studying the mathematical properties of a particular model for the biological phenomenon of T-cell activation with arbitrary values of the parameters. A complementary question is to what extent known experimental data on the parameters may further constrain the dynamics in this model. In addition, it is important to know whether this model is consistent with all biological data and how it compares to other possible models for the same biological system. For a discussion of this, we refer to [6,13,15]. It was indicated in [6] that the situation where C N is a decreasing function of ν 1 cannot be reproduced using the model of [4]. Our results indicate that a failure of the model to reproduce this effect must depend not only on the model itself but on the choice of parameters used for simulations. At the same time, it may be that this effect only occurs in experiments where the measurements are done on long time scales (many hours) and not on the time scale of the initial activation (a few minutes) for which the models of [3,4] were primarily intended. We plan to investigate these questions further elsewhere.
Data accessibility. This article has no additional data.