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

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 [12].) The central players in the immune system are the white blood cells (leukocytes) and it is important that these cells be 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 task of discrimination is complicated. An important element of the process of distinction between self and non-self is the activity of the class of leukocytes 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’ [4]. When it recognizes its antigen the T cell changes its behaviour and is said to be activated. In what follows we study a mathematical model for what happens in the first few minutes after a T cell recognizes its antigen.

In [1] 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 [5] the authors introduced a radically simplified version of the model of [1]. The new model includes the essential explanatory power of the old one while being much more transparent and tractable for analytical investigation. It also made some new predictions which were confirmed experimentally. In [5] 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 aim of the present paper is to obtain results about the qualitative behaviour of solutions of the model of [5] which are as general as possible. In Section 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 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\leq 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 Section 4. Let us consider the function $f(L_{1},\nu_{1})$, which expresses the degree of activation in terms of the parameters $L_{1}$ (concentration of agonist ligand) and $\nu_{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 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 $\nu_{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 [10].

The analysis of the response function is extended to cover the effects of the antagonist in Section 5. The last section is devoted to conclusions and an outlook.

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 an MHC molecule. Thus we talk about a pMHC complex as the object to be recognized. In the model of [5] 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 T cell receptor. When this happens certain proteins associated to the T cell receptor are phosphorylated, i.e. phosphate groups become attached to them. For simplicity we will describe this by saying that the receptor-pMHC complex is phosphorylated.

The reaction network for the model of [5] is shown in Figure 1.

The state variables will now be listed. The concentration of unphosphorylated complexes of the T cell receptor with the agonist will be denoted by $C_{0}$ and the concentration of unphosphorylated complexes of the T cell receptor with the antagonist will be 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 but it may be worth to note that a biologically reasonable value of $N$ could be as large as 20 while in the model of [1] we have $N=9$. $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 $\phi$ and dephosphorylation with rate constant $b$. They are also dephosphorylated by SHP-1 with rate constant $\gamma$ and dissociate with rate constants $\nu_{1}$ and $\nu_{2}$. Antigens bind to the receptor with rate constant $\kappa$. SHP-1 is activated by the singly phosphorylated complexes with rate constant $\alpha$ and deactivated with rate constant $\beta$. 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

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. Using these to eliminate the additional variables leads to the system (1)-(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 (1)-(7) this implies statements about the positivity of the quantities $S$, $C_{j}$ and $D_{j}$ and of the differences $S_{T}-S$, $R-\sum_{j=0}^{N}(C_{j}+D_{j})$, $L_{1}-\sum_{j=0}^{N}C_{j}$ and $L_{2}-\sum_{j=0}^{N}D_{j}$. Let us call the region where all these quantities are strictly positive the biologically feasible region. Note that due to the conservation laws this region is bounded. Let $\Sigma_{1}=\sum_{j=0}^{N}C_{j}$ and $\Sigma_{2}=\sum_{j=0}^{N}C_{j}$. Then it follows from (1)-(7) that

Note next that $C_{0}$ cannot be zero at an $\omega$-limit point. For if it were zero at such a point we could consider the solution passing through that point at a time $t_{0}$. The equation (2) would imply that $\dot{C}_{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 (3) with $j=1$ to show that $C_{1}$ cannot be zero at an $\omega$-limit point. This in turn allows us to prove using (1) that $S$ can never be zero at an $\omega$-limit point. In a similar way it can be concluded successively that $C_{2},\ldots,C_{N}$ and $D_{0},\ldots,D_{N}$ are positive at any $\omega$-limit point of a non-negative solution. This concludes the proof of the lemma.

The fact that all $\omega$-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\to\infty$ is strictly positive. When this last property holds the system is said to be persistent [2]. 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. [7], Chapter I, Theorem 8.2). Since 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 [5].

A question not addressed in [5] is whether there might exist more than one positive steady state for a fixed choice of parameters. In this section it will be 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 (8) and (9). Then $\frac{\partial f_{1}}{\partial\Sigma_{2}}$ and $\frac{\partial f_{2}}{\partial\Sigma_{1}}$ are negative and hence the system (8)-(9) is competitive. It follows that every solution of this system converges to a steady state as $t\to\infty$ [8].

A steady state $(\Sigma_{1}^{*},\Sigma_{2}^{*})$ of (8)-(9) satisfies the equations

We can solve for $\Sigma_{1}^{*}$ and $\Sigma_{2}^{*}$ as functions of $\Sigma_{1}^{*}+\Sigma_{2}^{*}$.

Hence

The function of $\Sigma_{1}^{*}+\Sigma_{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 $\Sigma_{1}^{*}+\Sigma_{2}^{*}=0$ and the opposite inequality holds for $\Sigma_{1}^{*}+\Sigma_{2}^{*}=\min\{L_{1}+L_{2},R\}$. Thus the equation has a unique solution for $\Sigma_{1}^{*}+\Sigma_{2}^{*}$ in the interval $[0,\min\{L_{1}+L_{2},R\}]$ From this it is possible to compute values of $\Sigma_{1}^{*}$ and $\Sigma_{2}^{*}$ which solve (10) and (11) and lie in the intervals $[0,\min\{L_{1},R\}]$ and $[0,\min\{L_{2},R\}]$, respectively. The quantities $\Sigma_{1}^{*}$ and $\Sigma_{2}^{*}$ are functions of the parameters $R$, $L_{1}$, $L_{2}$, $\kappa$, $\nu_{1}$ and $\nu_{2}$.

It can be concluded that the solution passing through an $\omega$-limit point of a solution of the original system satisfies a simplified system containing $\Sigma_{1}^{*}$ and $\Sigma_{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 is

This form of the equations is valid for $N\geq 3$. In the case $N=2$ it is still correct if it is taken into account that the condition $2\leq j\leq 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.

In order to get more information we will restrict in the remainder of this section to what we call the agonist-only case. This is obtained from the system (1)-(7) by setting $L_{2}$ and the $D_{i}$ to zero. There is a corresponding limiting system, which is obtained from (15)-(21) by setting $\Sigma_{2}^{*}$ and the $D_{i}$ to zero. In this case we write $\Sigma^{*}$ instead of $\Sigma_{1}^{*}$ for brevity. Consider the limiting system in the agonist-only case with $N=1$. This is

Solving the equation $\dot{S}=0$ for $C_{1}$ and substituting the result into the equation $\dot{C}_{1}=0$ gives the quadratic equation

Since 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 (24) that this root is less than $S_{T}$. Equation (23) implies that $C_{1}<\Sigma^{*}$ at a steady state and so these quantities can be completed to a steady state of the original system by defining $C_{0}=\Sigma^{*}-C_{1}$. The steady state is unique in this case.

In the case $N=2$ the equations are

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)=\sum_{k=0}^{N}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}<\Sigma^{*}$ and that the steady state is unique.

In the case $N=3$ the system is

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 $\nu_{1}$ in the system (28)-(31) and keep all others fixed. Note that these coefficients come from the parameters in the agonist-only case of (1)-(4). To obtain the desired variation of the coefficients we fix all parameters in (1)-(4) except $S_{T}$, $\nu_{1}$ and $\kappa$ and vary $\kappa$ in such a way that $\frac{\nu_{1}}{\kappa}$ does not change. This ensures that $\Sigma^{*}$ does not change. In fact we may simplify the calculations by setting $b=0$ since 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 $\nu_{1}$ depend on a parameter $\epsilon$ with both of them being positive for $\epsilon>0$. Suppose in addition that in the limit $\epsilon\to 0$ we have the asymptotic relations $S_{T}=\bar{S}_{T}\epsilon^{-1}+o(\epsilon^{-1})$ and $\nu_{1}=\bar{\nu}_{1}\epsilon^{4}+o(\epsilon^{4})$ for constants $\bar{S}_{T}$ and $\bar{\nu}_{1}$. Then we obtain asymptotic expansions $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}\epsilon^{3}+o(\epsilon^{3})$ for a constant $A_{0}<0$ and $a_{2}=A_{2}\epsilon^{-1}+o(\epsilon^{-1})$ for a constant $A_{2}<0$. Let $q(S)=\epsilon p(S)$. Then $q(1)$ converges to $A_{2}$ for $\epsilon\to 0$ and is thus negative for $\epsilon$ small enough. The same is true for $p(1)$. On the other hand

Hence for $\epsilon$ sufficiently small $p(\epsilon^{2})>0$. Putting these facts together shows that $p$ has three positive roots when $\epsilon$ 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}<\Sigma^{*}$ and defining $C_{0}=\Sigma^{*}-(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 since its derivative has no positive root. These results can be summed up as follows:

A concrete example of parameters for which there are three positive steady states is obtained by setting $\alpha$, $\beta$, $\gamma$, $\phi$, $L_{1}$ and $R$ equal to one and choosing $S_{T}=10$, $\kappa=2\times 10^{-4}$, $\nu_{1}=10^{-4}$. A computer calculation shows that the coordinates $(S^{*},C_{0}^{*},C_{1}^{*},C_{2}^{*},C_{3}^{*})$ of the steady states are approximately

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}$, see 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.

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\to\infty$. 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 (36). Thus if we make $\alpha$ and $\beta$ 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 $\beta$ large while fixing the other parameters we can arrange that the left hand side of (36) is greater than the right hand side. It follows that parameter values exist for which (36) holds. The reason why this is interesting is that the discriminant of the characteristic equation of the linearization is the sum of a term which vanishes when (36) holds and the expression $-4\alpha\gamma(S_{T}-S)C_{1}$. Thus when (36) 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 $\alpha=0$ and $S=0$. We refer to this as the kinetic proofreading system since it is closely related to McKeithan’s kinetic proofreading model [11]. 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 [13] 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 $\omega$-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 [3].

Consider now the full agonist-only system. Setting $\alpha=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 $\alpha=0$ which is on the boundary of the biologically feasible region and is a hyperbolic sink. Denote its coordinates by $(0,C_{j}^{*})$. For $\alpha$ small and positive there exists a hyperbolic sink which is a small perturbation of that for $\alpha=0$. It must be in the biologically feasible region since $C_{1}>0$ there and equation (1) would imply that $\dot{S}>0$ there if $S$ were negative. Thus for sufficiently small values of $\alpha$ there exists a positive steady state which is a hyperbolic sink $(S^{*}(\alpha),C_{j}^{*}(\alpha))$ close to $(0,C_{j}^{*})$. There exists a positive number $r$ such that for $\alpha$ sufficiently small, say $\alpha\leq\alpha_{0}$, $(S^{*}(\alpha),C_{j}^{*}(\alpha))$ is the only $\omega$-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 $\dot{h}\leq 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 $\dot{h}$ has a strictly negative maximum. We can consider the behaviour of the function $h$ for solutions of the system for positive $\alpha$. For small $\alpha$ it is still a Lyapunov function on the complement of a small ball about the steady state while there are no $\omega$-limit points except the steady state itself within that ball. Hence for $\alpha$ sufficiently small a solution can have no $\omega$-limit points other than the steady state. It follows that for $\alpha$ 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 $\omega$-limit point also has a unique steady state which is globally asymptotically stable for $\alpha$ sufficiently small. A similar argument applies in the case of the full system (1)-(7) since in that case the system obtained by setting $\alpha$ and $S$ to zero is just the product of two copies of the corresponding system in the agonist-only case.

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 $\nu_{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}$, $\nu_{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},\nu_{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 $\nu_{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 [5]. 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 [5].

The behaviour of the response function will be estimated in various parameter ranges. In order 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 (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 2 shows that for fixed parameters in (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 (8) in the form

The equation for $S$ can be solved to give the relation $S=S_{T}\frac{C_{1}}{C_{1}+C_{*}}$ with $C_{*}=\frac{\beta}{\alpha}$. Summing the expression for $C_{j}$ given in Lemma 2 over $j$ gives

The following equation relating $a_{-}$ and $a_{+}$ is equation (21) of [5].

Combining the last two equations gives

Having completed the necessary preliminaries we now proceed to study the qualitative behaviour of the reponse 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 (46) can be written in the form

Note that $\sum_{j=0}^{N}C_{j}\leq L_{1}$ so that this equation implies that

where $-\frac{1}{4}<q<0$. Using (48) it is possible to write down an explicit expression for $C_{N}$, namely

It follows from (49) that

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\to 0$. Now $C_{1}\leq\sum_{j=0}^{N}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\to 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}\to-\infty$.

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 2. To do this we compute the left hand side of the characteristic equation (39) for $r=\frac{\phi}{\phi+\nu_{1}}$. The result is $-\frac{\phi\nu_{1}(b+\gamma S)}{(\phi+\nu_{1})^{2}}<0$. It follows that $r_{-}<\frac{\phi}{\phi+\nu_{1}}$. Hence $1-r_{-}>\frac{\nu_{1}}{\phi+\nu_{1}}$. Substituting this into (49) gives $a_{-}>\frac{\nu_{1}}{\phi+\nu_{1}}\left(\sum_{j=0}^{N}C_{j}\right)$. Note that $\frac{S}{S_{T}}\geq\min\left\{\frac{C_{1}}{2C_{*}},\frac{1}{2}\right\}$. Hence a positive lower bound for $C_{1}$ implies a positive lower bound for $\frac{S}{S_{T}}$.

Next we will derive a lower bound for $\gamma S$ in the case that $S_{T}$ is large. This will be proved by contradiction. Suppose that $\gamma S\leq\rho$ for some $\rho>0$. Then it follows from the characteristic equation that $r_{-}\geq\frac{\phi}{\phi+\rho+\nu_{1}}$. Using this in the equation for $C_{1}$ gives $C_{1}\geq\frac{\phi\nu_{1}}{(\phi+\nu_{1})(\phi+\rho+\nu_{1})}\left(\sum_{j=0}^{N}C_{j}\right)$. It follows that

It is then clear that for a given value of $\rho$ 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 $\rho>0$ there is a lower bound for $S_{T}$ which implies that $\gamma S\geq\rho$. It is convenient to make the restrictions that $\kappa R\geq 1$ and $L_{1}/R\leq 1$ since then it is possible to replace $\sum_{j=0}^{N}C_{j}$ in (55) by $\frac{3L_{1}}{4(1+\nu_{1})}$ by using (51).