## Abstract

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 (TCR) 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 TCR 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 a maximal activation of the T cell. This suggests that mechanisms taking place in the first few minutes after activation and included in the model studied in this paper suffice to explain the optimal dissociation time seen in experiments. 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.

## 1. 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. (For background information on immunology, we refer to [1].) The central players in the immune system are the white blood cells (leucocytes) and it is important that these cells are able to distinguish between dangerous substances and host tissues. This is often referred to as the distinction between non-self and self. A failure to combat dangerous substances may lead to infectious diseases becoming life-threatening. On the other hand, if the immune system attacks host tissues, this can lead to autoimmune disease. The process of discrimination between self and non-self is complicated, involving numerous mechanisms. An important element of this process, which is investigated in the present paper, is the activity of the class of leucocytes called T cells. An individual T cell is supposed to recognize a particular substance (antigen) and take suitable action if that substance is dangerous. Recognition is based on the binding of the antigen to a molecule on the T-cell surface, the T-cell receptor (TCR). It is believed that the most important aspect of this process is the time the antigen remains bound before being released (the dissociation time), an idea which has been called the ‘lifetime dogma’ [2]. When it recognizes its antigen, the T cell changes its behaviour and is said to be activated. In this work, we study a mathematical model of the first few minutes of T-cell activation after the TCR binds to its antigen.

In [3], Altan-Bonnet and Germain introduced a model for the initial stage of T-cell activation. Simulations using this model gave results which fitted a number of experimental findings. On the other hand, it was too elaborate to be readily accessible to a mathematical analysis of its dynamics. In [4], the authors introduced a radically simplified version of the model of [3]. The new model includes the essential explanatory power of the old one while being much more transparent and tractable for analytical investigation. It also predicts features of experimental data which had not been explained previously, such as the fact that the response of a T cell can decrease as a function of the amount of antigen when the concentration of the phosphatase SHP-1 is sufficiently high. In [4], a number of interesting analytical calculations were performed, but the mathematical conclusions which can be drawn from these were not worked out in detail.

The relations between these two models will now be explained briefly. The TCR is associated with other proteins (CD3 and the *ζ*-chain), forming the TCR complex. These other proteins have cytoplasmic tails on which there are regions called immunoreceptor tyrosine-based activation motifs (ITAMs). Each ITAM contains two tyrosines on which it can be phosphorylated (i.e. phosphate groups can become bound to these tyrosines) separated by a few other amino acids. Phosphorylation of the ITAMs is a typical sign of T-cell activation. In the TCR complex, there are 10 ITAMs and thus a total of 20 phosphorylation sites of potential importance for the activation of the T cell. In a later step of the process, the protein ZAP-70 binds to the doubly phosphorylated ITAMs of the *ζ*-chains and itself becomes phosphorylated. There are two *ζ*-chains in the T-cell complex and each contains three ITAMs. Thus a total of six further phosphorylations are possible. The exact order in which all these sites are phosphorylated is not understood in detail and so this part of the system is treated in a rather schematic way in the models. In the model of [4], it is assumed that there are *N* sites which are phosphorylated sequentially, i.e. in a particular order. ZAP-70 is not included in the model. In the simulations, the choice *N*=5 is made. In the model of [3], the phosphorylation sites included are those of one *ζ*-chain and ZAP-70, leading to a total of nine. Both models include a negative feedback acting through the phosphatase SHP-1, which can dephosphorylate the sites just discussed. The importance of SHP-1 in controlling T-cell activation was pointed out in [5].

The other main difference between the models of [3,4] is the treatment of events downstream of the process of phosphorylation of the receptor complex. In [3], phosphorylation of ZAP-70 leads to a chain of events culminating in the activation by double phosphorylation of extracellular signal-regulated kinase (ERK). There is also a positive feedback loop from ERK through SHP-1 to the receptor complex. The positive feedback loop is absent from [4] and is thus seen to be unnecessary for explaining the main effects studied in [3]. In [3], it was found that the reactions linking phosphorylation of the receptor to the activation of ERK act as a switch: when the concentration of the phosphorylated receptor complex exceeds a certain threshold, ERK becomes activated. In the model of [4], this switch is incorporated in the form that when the concentration of the maximally phosphorylated form of the receptor complex exceeds a certain threshold, this is taken as the defining property of the T cell being activated.

The aim of the present paper is to obtain results about the qualitative behaviour of solutions of the model of [4] which are as general as possible. In §2, the model is defined and some of its basic properties are derived. The model describes a situation where both an agonist (the antigen which should be recognized) and an antagonist (a competing antigen) are present. Section 3 is concerned with the number of steady states and their stability. After some general results have been derived, the discussion turns to more detailed properties of the solutions in the case that the antagonist is absent and treats cases where the number *N* of phosphorylation sites included in the model is small. In particular, it is 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.

## 2. 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:
*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*,

### 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* *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*=*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 *L*_{1} and equation (2.9) implies that *L*_{2}. Summing (2.8) and (2.9) shows that *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 *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

## 3. 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

A steady state *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 *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 *C*_{0} and *D*_{0} can be eliminated from this system in favour of the other *C*_{j} and *D*_{j}. The result is
*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 *D*_{i} to zero. In this case, we write *Σ** instead of *N*=1. This is
*C*_{1} and substituting the result into the equation *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 are
*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 *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 is
*N*=3 analogous to those already done gives a quartic polynomial. Its coefficients are
*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 *a*_{4}=*A*_{4}, *a*_{3}=*A*_{3}+*o*(1), *a*_{1}=*A*_{1}+*o*(1) for positive constants *A*_{4}, *A*_{3} and *A*_{1}, *a*_{0}=*A*_{0}*ϵ*^{3}+*o*(*ϵ*^{3}) for a constant *A*_{0}<0 and *a*_{2}=*A*_{2}*ϵ*^{−1}+*o*(*ϵ*^{−1}) for a constant *A*_{2}<0. Let *q*(*S*)=*ϵp*(*S*). 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,
*ϵ* 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}, *C*_{1}+*C*_{2}+*C*_{3}<*Σ** and defining *C*_{0}=*Σ**−(*C*_{1}+*C*_{2}+*C*_{3}) 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 *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 *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 *α* 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 that *S* were negative. Thus for sufficiently small values of *α*, there exists a positive steady state which is a hyperbolic sink *r* such that for *α* sufficiently small, say *α*≤*α*_{0}, *ω*-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 *r* about the steady state the function *α*. 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.

## 4. 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 *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 *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
*C*_{j} defined by (4.1) satisfy the steady state form of equation (2.3). That *r*_{−}<1<*r*_{+} can be seen by noting that the function on the left-hand side of (4.3) is negative at *r*=1. The condition that the quantities *C*_{j} satisfy the equations (2.2)–(2.4) with *r*_{±} as in (4.2) and certain coefficients *a*_{−} and *a*_{+} together with the equations obtained by substituting (4.1) into the equations *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
*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
*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 left-hand 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
*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
*a*_{−} and *a*_{+} is equation (21) of [4]:

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 *C*_{N}, namely
*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 *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

Next we look at an intermediate regime where the amount of activated SHP-1 is well away from both zero and *S*_{T}. As a first step, we obtain an estimate for *r*_{−} which is sharper than that in lemma 4.1. To do this, we compute the left-hand side of the characteristic equation (4.3) for *r*=*ϕ*/(*ϕ*+*ν*_{1}). The result is −*ϕν*_{1}(*b*+*γS*)/(*ϕ*+*ν*_{1})^{2}<0. It follows that *r*_{−}<*ϕ*/(*ϕ*+*ν*_{1}). Hence 1−*r*_{−}>*ν*_{1}/(*ϕ*+*ν*_{1}). Substituting this into (4.13) gives *C*_{1} implies a positive lower bound for *S*/*S*_{T}.

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 *ρ* 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 *L*_{1}/4(1+*ν*_{1}) using (4.15).

From (4.2), it can be concluded that
*η*=(*ϕ*+*ν*_{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
*C*_{1}=*O*(*η*) and the expression relating *S* and *C*_{1} then shows that

These relations indicate that in leading order *r*_{−} is proportional to *S*. However, it is also the case that
*r*_{−} is proportional to *S*^{−1}. Hence
*r*_{−} gives
*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.

### Theorem 4.2

*Consider the response function* *for the steady states of the system* (*2.1*)–(*2.4*) *with L*_{2}=0 *and D*_{j}=0. *Choose fixed values for all parameters in the system except L*_{1} *and S*_{T}. *Suppose that κR*≥1. *Let ϵ*>0. *Then for any sufficiently small constant δ*>0, *the following holds. If* 0<*L*_{0}<*δ, there exists μ*>0 *such that if S*_{T}≥*μ the inequality*
*holds on the interval*

### 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 *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 *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].

## 5. 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 *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 *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:
*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:
*η*′ 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].

## 6. 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.

## Authors' contributions

A.D.R. conceived the study and contributed to the mathematical analysis. E.D.S. contributed to the mathematical analysis and carried out simulations. Both authors contributed to the writing of the text and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

A.D.R. work did not receive any specific funding. E.D.S. is supported, in part, by grants AFOSR FA9550-14-1-0060 and ONR N00014-13-1-0074.

- Received July 3, 2017.
- Accepted October 6, 2017.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.