\documentclass[reqno]{amsart}
\usepackage{hyperref}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2010(2010), No. 115, pp. 1--10.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2010/115\hfil Model for T cells and macrophages]
{Analysis of a mathematical model for interactions between T cells and
macrophages}
\author[A. D. Rendall \hfil EJDE-2010/115\hfilneg]
{Alan D. Rendall}
\address{Alan D. Rendall \newline
Max-Planck-Institut f\"ur Gravitationsphysik\\
Albert-Einstein-Institut\\
Am M\"uhlenberg 1\\
14476 Potsdam, Germany}
\email{rendall@aei.mpg.de}
\thanks{Submitted August 5, 2010. Published August 19, 2010.}
\subjclass[2000]{37N25, 92C37}
\keywords{Dynamical systems; qualitative behaviour; immunology}
\begin{abstract}
The aim of this article is to carry out a mathematical analysis
of a system of ordinary differential equations introduced by
Lev Bar-Or to model the interactions between T cells and
macrophages. Under certain restrictions on the parameters of
the model, theorems are proved about the number of stationary
solutions and their stability. In some cases the existence of
periodic solutions or heteroclinic cycles is ruled out.
Evidence is presented that the same biological phenomena could
be equally well described by a simpler model.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\section{Introduction}
Autoimmune diseases result in a great deal of suffering for
affected individuals and huge costs for society. Notable examples
are multiple sclerosis, rheumatoid arthritis and type I diabetes.
In these diseases the ability of the immune system to distinguish
between self and non-self is compromised, with the result that
host tissues are attacked and damaged. It is important to get a
better understanding of the processes involved and one way to do
so is to introduce theoretical models of the immune system, in
particular mathematical models.
There has been a lot of work on mathematical modelling of
interactions of the immune system with pathogens. See for instance
\cite{perelson99} and the book of Nowak and May \cite{nowak} which
concentrate on the case of HIV and present various models. Some
deeper mathematical analysis of these models can be found in
\cite{deleenheer03}, \cite{korobeinikov04} and \cite{egami09}.
Autoimmune diseases do not need to involve any pathogens, although
pathogens might contribute to them indirectly, for instance by
molecular mimicry \cite{oldstone}. For this reason it would be
interesting to have models for the intrinsic workings of the
immune system where non-self antigens play no direct role, in the
sense that they are not included among the dynamical variables.
Apparently few models of this type exist in the literature. One
example is introduced in a paper of Lev Bar-Or \cite{levbaror}. It
is a system of four ordinary differential equations which
describes the interactions of T cells and macrophages by means of
the cytokines they produce. The aim of this paper is to
investigate what can be said about the properties of the solutions
of this system on the level of mathematical proofs.
Some relevant concepts from immunology will now be introduced. For
a comprehensive introduction to the subject see \cite{murphy}. T
cells are white blood cells which mature in the thymus. One type
of T cell, the T helper cell, helps to direct the activity of
other immune cells. These cells are also known as CD4${}^+$ since
they carry the surface molecule CD4. T cells communicate with
other cells by secreting soluble substances known as cytokines.
Another important type of blood cell is the macrophage which
ingests pathogens and cell debris through phagocytosis.
Macrophages also secrete cytokines. Both T cells and macrophages
react in various ways to the cytokines which are present in their
surroundings. The list of known cytokines is long and each of them
has its own characteristics in terms of which types of cells
secrete it and what effects it has on cells which detect its
presence. An idea of the complexity of this signalling system can
be obtained from \cite{levbaror}.
It is common to distinguish between two types of T helper cells,
known as Th1 and Th2, according to the cytokines they produce.
There may be overlaps but roughly speaking it may be supposed that
there is one set of cytokines (e.g. interferon $\gamma$) which are
called type 1 and are produced by Th1 cells and another (e.g.
interleukin 4) called type 2 which are produced by Th2 cells.
Macrophages produce cytokines of both types. The basic quantities
in the equations of \cite{levbaror} are average concentrations
corresponding to type 1 and type 2 cytokines produced by T cells
and macrophages. These four concentrations are denoted by $C^T_1$,
$C^T_2$, $C^M_1$ and $C^M_2$. The populations of different cell
types do not occur directly in the system. The quantities which
have a chance of being measured directly are
$C_1=\frac12(C^T_1+C_1^M)$ and $C_2=\frac12(C^T_2+C_2^M)$.
There are cases where the immune response is dominated by either
Th1 or Th2 cells. This may be important in order to effectively
combat a particular pathogen. For instance a sufficiently strong
Th1 response is necessary for containing or eliminating a
tuberculosis infection \cite{wigginton}, \cite{marino}. An
inappropriate balance between these two states can also contribute
to autoimmune disorders. It has, for instance, been suggested that
multiple sclerosis is associated with an immune response which is
biased towards Th1. See \cite{lassmann} for a critical review of
this idea. In the context of the model it is said that there is
Th1 or Th2 dominance if $C_1>C_2$ or $C_2>C_1$, respectively.
Another important function of macrophages which plays a role in
the model of \cite{levbaror} is the presentation of antigens.
Small peptides which result from the digestion of material taken
up by a macrophage are presented on its surface in conjunction
with MHC II molecules. (Major histocompatibility complex of class
II.) This can stimulate T cells which come into contact with the
macrophage.
It has not been possible to give a complete analysis of the
dynamics of the system of \cite{levbaror}. A certain inequality on
the parameters of the system gives rise to a regime where there is
a unique stationary solution which acts as an attractor for all
solutions as $t\to\infty$. This is proved in Theorem \ref{thm1}. For a more
restricted set of parameters it is possible to show (Theorem \ref{thm2})
that each solution converges to a stationary solution which in
general depends on the solution considered. With further
restrictions on the parameters it is shown that there are between
one and three stationary solutions and the subsets of parameters
for which different numbers of stationary solutions occur are
described. This is the content of Theorem \ref{thm3}.
In Theorems \ref{thm2} and \ref{thm3}
the coefficients describing antigen presentation are set to zero.
One situation where information can be obtained on the dynamics
including antigen presentation is analysed in Theorem \ref{thm4}.
The analysis which has been done has uncovered no evidence that
the inclusion of the effect of macrophages makes an essential
difference to the behaviour of solutions. The types of qualitative
behaviour which have been proved to occur in this paper and those
which are shown in the figures in \cite{levbaror} can be found in
a truncated system which only includes the effect of T cells. This
is discussed in section \ref{Tcellonly}.
\section{Analysis of the dynamical system}
In what follows it will be convenient to use a notation which is
more concise than that of \cite{levbaror}. Let $x_1$, $x_2$,
$x_3$, $x_4$, $z_1$ and $z_2$ denote $C^T_1$, $C^T_2$, $C^M_1$,
$C^M_2$, $x_1+x_3$ and $x_2+x_4$ respectively. The basic dynamical
system is:
\begin{equation}\label{basic}
\frac{dx_i}{dt}=-d_ix_i+g(h_i);\quad i=1,2,3,4.
\end{equation}
The $d_i$ are positive constants. The function $g$ is given
by
\begin{equation}
g(x)=\frac12 (1+\tanh (x-\theta))
\end{equation}
where $\theta$ is a constant. It satisfies the relations
$g(x+\theta)+g(-x+\theta)=1$ and $g'(x)=2g(x)(1-g(x))$. Hence
$g(\theta)=1/2$ and $g'(\theta)=1/2$. The functions $h_i$
are defined by $h_i=\sum_{j}a_{ij}x_j$ for some constants
$a_{ij}$. The coefficients in \eqref{basic} satisfy the following
conditions
\begin{enumerate}
\item $a_{ij}=b_{ij}+c_{ij}$, $i=1,2$, for some coefficients
$b_{ij}$ and $c_{ij}$
\item $b_{1j}=-b_{2j}$ for all $j$.
\item $b_{11}>0$, $b_{13}>0$, $b_{12}<0$ and $b_{14}<0$
\item The ratio $c_{2j}/c_{1j}$ is independent of $j$ and positive.
\item $c_{11}>0$ and $c_{13}>0$
\item $a_{3j}=-a_{4j}$ for all $j$
\item $a_{31}>0$, $a_{33}>0$, $a_{32}<0$, $a_{34}<0$
\end{enumerate}
The sign conditions encode the fact that the effect of type 1
cytokines on cells is to increase their production of type 1
cytokines and to decrease their production of type 2 cytokines,
while the effect of type 2 cytokines is exactly the opposite. The
coefficients $c_{ij}$ encode the effects of antigen presentation.
This is the only role of antigens in the model and the intensity
of stimulus due to antigen is not a dynamical degree of freedom in
the model. No assumption is made on the signs of $c_{12}$ and
$c_{14}$. There are eighteen parameters in the model which are
only constrained by some positivity conditions. The quantities
$z_1$ and $z_2$ represent total concentrations of type 1 and type
2 cytokines. The biologically relevant region $\mathcal B$ is that
where all the $x_i$ are non-negative. The function $g$ is strictly
positive. Hence if one of the variables $x_i$ vanishes at some
time its derivative at that time is strictly positive. It follows
that $\mathcal B$ is positively invariant under the evolution
defined by the system. If some $x_i$ is greater than or equal to
$d^{-1}_i$ on some time interval then $x_i$ is decreasing at a
uniform rate during that time. It follows that all solutions exist
globally to the future and enter the region ${\mathcal B}_1$
defined by the inequalities $x_i\le d^{-1}_i$ after finite time.
Thus in order to study the late-time behaviour of any solution
starting in $\mathcal B$ it is enough to consider solutions
starting in ${\mathcal B}_1$. In fact any solution enters the
interior of ${\mathcal B}_1$ after finite time.
\begin{theorem} \label{thm1}
For any value of the parameters the system \eqref{basic} has at
least one stationary solution. If
\begin{equation}\label{contractive}
\sup_i\sum_j |a_{ij}|< 2\inf_i d_i.
\end{equation}
then there is only one stationary solution and all solutions
converge to it as $t\to\infty$.
\end{theorem}
\begin{proof}
That the system always has at least one stationary solution
follows from the Brouwer fixed point theorem, cf.
\cite[Theorem I.8.2]{hale}. A stationary solution of the system
satisfies
\begin{equation}
x_i=d_i^{-1}g(h_i).
\end{equation}
To prove the second part of the theorem consider the estimate
\begin{equation}
\begin{aligned}
|d_i^{-1}g(h_i(y))-d_i^{-1}g(h_i(x))|
&\le (2d_i)^{-1}|h_i(y)-h_i(x)|\\
&\le (2d_i)^{-1}\sum_j |a_{ij}||y_j-x_j|.
\end{aligned}
\end{equation}
The first of these inequalities uses the mean value theorem and
the fact the derivative of $g$ is nowhere greater than one half.
If \eqref{contractive} holds then the mapping $\{x_i\}\mapsto
\{d_i^{-1}g(h_i)\}$ maps ${\mathcal B}_1$ to itself and is a
contraction in the maximum norm. Hence it has a unique fixed
point. It follows that when the coefficients of the system satisfy
the restriction \eqref{contractive} the system has exactly one
stationary solution. If $x(t)$ and $y(t)$ are two solutions then
under the assumption \eqref{contractive} it can be shown that
$|x-y|$ decays exponentially as $t\to\infty$. Thus all solutions
converge to the unique stationary solution as $t\to\infty$. A
related statement is that if $x_*$ is the unique stationary
solution then $|x-x_*|^2$ is a Lyapunov function.
\end{proof}
A limiting case of \eqref{basic} is that where all
$c_{ij}$ vanish. This will be called the zero MHC system. The
pattern of signs in the coefficients in the zero MHC system is
such that the change of variables $\tilde x_i=(-1)^{i+1}x_i$ leads
to a system $\frac{d\tilde x_i}{dt}=f(\tilde x_j)$ satisfying
$\frac{\partial f_i}{\partial \tilde x_j}> 0$ for all $i\ne j$.
This means that the dynamical system for the $\tilde x_i$ is
cooperative \cite{hirsch}. This pattern of signs corresponds to
what is referred to as a \lq community with limited
competition\rq\ or \lq competing subcommunities of mutualists\rq\
in \cite{smith86}. It implies, using the Perron-Frobenius theorem,
that the linearization at any point of the vector field defining
\eqref{basic} has a real eigenvalue of multiplicity one which is
greater than the modulus of any other eigenvalue. Notice that in
the special case of the zero MHC system where all $d_i$ are equal
to $d$ the linearization has two eigenvalues equal to $-d$. As a
consequence all the eigenvalues of the linearization must be real
in this case. It is not clear how these facts about the eigenvalue
structure can be used to help to understand the dynamics. The
condition on the coefficients $d_i$, which says that the rate of
degradation of different cytokines is exactly equal, would not be
true with biologically motivated parameters. Nevertheless it is
not unreasonable to assume that these coefficients are
approximately equal and that the model with equal coefficients is
not a bad approximation to the situation to be modelled. An
assumption of this type is made in the model of \cite{marino}.
Call the system obtained by assuming $c_{ij}=0$, $d_1=d_2$,
$d_3=d_4$ and $\theta=0$ in \eqref{basic} System 2 while
\eqref{basic} itself is System 1. The condition $\theta=0$ implies
that $g(x)+g(-x)=1$. Adding the equations for $i=1$ and $i=2$ then
gives
\begin{equation}
\frac{d}{dt}(x_1+x_2)=-d_1(x_1+x_2)+1.
\end{equation}
Similarly
\begin{equation}
\frac{d}{dt}(x_3+x_4)=-d_3(x_3+x_4)+1.
\end{equation}
These equations can easily be analysed. An invariant manifold
$S_1$ is defined by $x_1+x_2=d_1^{-1}$ and $x_3+x_4=d_3^{-1}$.
Substituting this back into the full system gives a
two-dimensional system written out below which will be called
System 3. The $\omega$-limit set of any solution of System 2 is
contained in the invariant manifold $S_1$. Passing from System 1
to System 3 means setting some parameters to zero and then
restricting to a two-dimensional invariant manifold of the
resulting system. System 3 may not include all the most
interesting dynamics exhibited by solutions of the system of
\cite{levbaror} but can be perturbed to get information about
cases where the effect of the MHC is small but non-zero. The phase
portrait of System 2 follows immediately from that of System 3.
The explicit form of System 3 is
\begin{gather}\label{system3}
\frac{dx_1}{dt}=-d_1 x_1+g((a_{11}-a_{12})x_1
+(a_{13}-a_{14})x_3+a_{12}d_1^{-1}+a_{14}d_3^{-1})\\
\frac{dx_3}{dt}=-d_3 x_3+g((a_{31}-a_{32})x_1
+(a_{33}-a_{34})x_3+a_{32}d_1^{-1}+a_{34}d_3^{-1})
\end{gather}
Note that the coefficients of $x_1$ and $x_3$ in these equations
which are linear combinations of the $a_{ij}$ are all positive.
The constant terms are negative. This is a cooperative system.
This fact together with the fact that the dimension of the system
is two implies that each solution converges to a limit as
$t\to\infty$ \cite{hirsch}. In other words the $\omega$-limit set
of each solution is a single point. Furthermore, there are no
homoclinic orbits or heteroclinic cycles. What is not clear in
general is how many stationary solutions there are. Some of the
conclusions of the above discussion can be summarized as follows.
\begin{theorem} \label{thm2}
If $c_{ij}=0$ for all $i$, $j$, $d_1=d_2$, $d_3=d_4$ and
$\theta=0$ then every solution of \eqref{basic} converges to a
stationary solution satisfying $x_1+x_2=d_1^{-1}$ and
$x_3+x_4=d_3^{-1}$ as $t\to\infty$. There are no homoclinic orbits
or heteroclinic cycles.
\end{theorem}
Consider now the special case of System 3 obtained by setting
$a_{1j}=a_{3j}$ for all $j$ and $d_1=d_3$. Call it System 4. This
corresponds to the assumptions that the T cells and macrophages
have identical properties with respect to their death rate and
their interactions with cytokines. This system consists of two
copies of a single equation. The quantity $x_1-x_3$ decays
exponentially. The further assumption $x_1=x_3$ on the initial
data, which means that there are equal numbers of T cells and
macrophages, gives rise to a single ODE, call it the toy model. It
is the same equation which occurs twice in System 4. It is of the
form
\begin{equation}\label{pretoy}
\frac{dx}{dt}=-d_1x+g(ax-b)
\end{equation}
where $a=a_{11}-a_{12}+a_{13}-a_{14}$ and
$b=-(a_{12}+a_{14})d_1^{-1}$. Note that $a$ and $b$ are positive.
Here the zero MHC condition has been used. This can be simplified
by defining $x'=d_1 x$, $t'=d_1 t$ and $a'=a/d_1$.
Suppressing the primes leads to
\begin{equation}\label{toy}
\frac{dx}{dt}=-x+g(ax-b)
\end{equation}
Denote the right hand side of \eqref{toy} by $h(x)$, suppressing
the parameter dependence. Stationary solutions of the toy model
are given by solutions of the equation
\begin{equation}\label{toystat}
g(ax-b)=x.
\end{equation}
Since $02$. For $a>2$ call the zeroes of the derivative $x_1$ and $x_2$
with $x_12$ there is precisely one value of $b$ for which the
equation $-x_1+g(ax_1-b)=0$ holds. Call it $b_1(a)$. Define
$b_2(a)$ similarly, replacing $x_1$ by $x_2$. Then $b_1(a)$ and
$b_2(a)$ are the only values of $b$ for which there is a
simultaneous solution of $g(ax-b)=x$ and $ag'(ax-b)=1$. For all
$a>2$ the inequality $b_1(a)a/2$ it satisfies $x>1/2$. If $a>2$ and there are three
stationary points then the central (unstable) one satisfies
$x<1/2$ for $b1/2$ for $b>a/2$.
When there is only one solution it is a hyperbolic sink. When
there are three solutions the two outermost are hyperbolic sinks
while the intermediate one is a hyperbolic source. If the
parameter $b$ is held fixed at some value less than one and $a$ is
varied then at the value of $a$ where $b_1(a)=b$ there is a fold
bifurcation. To see this note that
\begin{gather}
\frac{\partial h}{\partial x}=-1+ag'(ax-b) \\
\frac{\partial h}{\partial a}= xg'(ax-b)=2xg(ax-b)(1-g(ax-b)) \\
\frac{\partial^2 h}{\partial x^2}=a^2g''(ax-b)
=4a^2g'(ax-b)(1-2g(ax-b))
\end{gather}
From the first equation it follows that at a critical point $x$
the quanitity $ag'(ax-b)$ must be equal to one. Hence the
derivative with respect to $a$ does not vanish there. The second
derivative with respect to $x$ can only vanish there if
$g(ax-b)=1/2$, which implies that $x=\frac{b}{a}$, $x=1/2$ and
$b=1$. The situation for $b>1$ is similar to that for $b<1$. To
understand the case $b=1$ note that there is a cusp bifurcation in
$(a,b)$ at the point $(2,1)$. This can be shown using the facts
that
\begin{gather}
\frac{\partial h}{\partial b}=-g'(ax-b), \\
\frac{\partial^2 h}{\partial x\partial b}=-ag''(ax-b), \\
\frac{\partial^2 h}{\partial x\partial a}=axg''(ax-b)
+g'(ax-b), \\
\begin{aligned}
\frac{\partial^3 h}{\partial x^3}
&=g'''(ax-b) \\
&=4a^3g''(ax-b)(1-2g(ax-b))-8a^3(g'(ax-b))^2.
\end{aligned}
\end{gather}
Hence $\frac{\partial^3 h}{\partial x^3}$ and $\frac{\partial
h}{\partial a}\frac{\partial^2 h}{\partial x\partial b}-
\frac{\partial h}{\partial b}\frac{\partial^2 h}{\partial
x\partial a}$ are both non-vanishing and
\cite[Theorem 8.1]{kuznetsov} applies.
If we remove the condition $x_1=x_3$ then $x_1-x_3$ decays
exponentially. The qualitative nature of the dynamics of System 4
with the given restrictions on the parameters is then clear. This
in turn gives information about the phase portrait of System 2
with these values of the parameters. Converting back to the
original variables has the effect of replacing $a$ by $a/d_i$ in
these criteria. Some of these results are summarized in the
following theorem.
\begin{theorem} \label{thm3}
If in Theorem \ref{thm2} it is additionally assumed that
$a_{1j}=a_{3j}$ for all $j$ and $d_1=d_3$ then the number of
stationary solutions of \eqref{basic} is between one and three for
any values of the parameters. It is one whenever
$a_{11}-a_{12}+a_{13}-a_{14}\le 2$. There is a non-empty open set
where it is three. The set where it is two is a union of two
smoothly embedded curves which is the boundary between the sets
where it is one and three.
\end{theorem}
Consider a choice of parameter set for \eqref{basic} for which
there are three hyperbolic stationary points. It follows by a
simple stability argument that there is an open neighbourhood of
this parameter set in parameter space such that for all parameters
in this neighbourhood there are precisely three hyperbolic
stationary points. If the original set of stationary points
consists of two sinks and one saddle then there is an open
neighbourhood where this property persists.
In the context of this theorem the question of Th1 or Th2
dominance can be examined for stationary solutions. It is of
interest to know whether changing the parameters in the system can
cause a switch from one type of dominance to the other. If this
happens there must be some values of the parameters for which
there is a stationary solution with $z_1=z_2$. Under the
hypotheses of the Theorem \ref{thm2}, $x_2=d_1^{-1}-x_1$ and
$x_4=d_1^{-1}-x_3$ for any stationary solution. Hence
$z_2=2d_1^{-1}-z_1$. If $z_1=z_2$ then this implies that
$z_1=d_1^{-1}$ and $x_1=\frac12 d_1^{-1}$. In terms of the new
variable introduced in \eqref{toy} this means that $x=1/2$.
The above discussion gives information about Th1 or Th2 dominance
for various parameter values satisfying the restrictions in the
statement of Theorem \ref{thm3}. What is of most interest for the
applications is the position of the stable stationary points. In
particular it is clear that there exist parameter values where
there are stationary solutions with both types of dominance.
Now another set of simplified versions of the system will be
studied. The phase portraits given in Fig. 2 in \cite{levbaror}
relate to one special case of this kind. In the notation used here
it is defined by the following restrictions on the parameters:
$d_i=1$ for all $i$, $\theta=0$, all coefficients $b_{ij}$ with
$j=1$ or $j=3$ are equal, all coefficients $a_{ij}$ with $j=2$ or
$j=4$ are equal, the coefficients $c_{ij}$ are equal for all $i$
and $j$. This reduces the number of free parameters to three. To
simplify the notation let $A=b_{11}$, $B=-b_{12}$, $C=c_{11}$.
These are all positive constants. For Fig. 2 in \cite{levbaror}
two sets of values for the coefficients are considered. In both
cases $C=0.5$. In Fig. 2a $A=0.4$ and $B=0.5$ while in Fig. 2b
$A=0.6$ and $B=0.65$. With these assumptions System 1 becomes:
\begin{gather}%\label{ABCsystem}
\frac{dx_1}{dt}=-x_1+g((A+C)(x_1+x_3)+(-B+C)(x_2+x_4))\label{abc1}\\
\frac{dx_2}{dt}=-x_2+g((-A+C)(x_1+x_3)+(B+C)(x_2+x_4))\label{abc2}\\
\frac{dx_3}{dt}=-x_3+g(A(x_1+x_3)-B(x_2+x_4))\label{abc3}\\
\frac{dx_4}{dt}=-x_4+g(-A(x_1+x_3)+B(x_2+x_4))\label{abc4}
\end{gather}
Call this System 5. If $A+B+2C<1$ then there is at most one
stationary point of this system. This can be proved in the same
way as Theorem \ref{thm1}. Adding the equations in pairs leads to a closed
system for the two experimentally accessible quantities $z_1$ and
$z_2$. These are the variables which are plotted in Fig. 2 of
\cite{levbaror}.
\begin{gather}
\frac{dz_1}{dt}=-z_1+g((A+C)z_1+(-B+C)z_2)+g(Az_1-Bz_2)\label{z1evol} \\
\frac{dz_2}{dt}=-z_2+g((-A+C)z_1+(B+C)z_2)+g(-Az_1+Bz_2)\label{z2evol}
\end{gather}
Call this System 6. Denote the right hand sides of \eqref{z1evol}
and \eqref{z2evol} by $f_1(z_1,z_2)$ and $f_2(z_1,z_2)$
respectively. For this system uniqueness of the stationary
solution follows from the assumption that the coefficients satisfy
$A+B+C<1$. If in addition the coefficient $C$ is assumed to vanish
then this reduces to
\begin{gather}
\frac{dz_1}{dt}=-z_1+2g(Az_1-Bz_2)\label{z1evol3}\\
\frac{dz_2}{dt}=-z_2+2g(-Az_1+Bz_2)\label{z1evol4}
\end{gather}
Call this System 7. Note that the system obtained by setting $C=0$
in System 5 is equivalent to a special case of System 2 and that
Theorems \ref{thm2} and \ref{thm3} apply to it.
The parameters are related by
$a=2(A+B)$, $b=2B$. The special cases in \cite[Fig. 2]{levbaror}
correspond to $a=1.8$, $b=1$ and $a=2.5$, $b=1.3$. It follows from
\eqref{z1evol3} and \eqref{z1evol4} that
$\frac{d}{dt}(z_1+z_2)=-(z_1+z_2)+2$ and that $z_1+z_2$ tends to
two as $t\to\infty$. Thus the dynamics is controlled by that on
the invariant manifold $z_1=2-z_2$. With the choices which have
been made it coincides with the invariant manifold $S_1$
introduced earlier.
Now consider what happens when $C\ne 0$. Provided
$C\le\min\{A,B\}$ the system consisting of \eqref{z1evol} and
\eqref{z2evol} is competitive. Thus each solution $(z_1,z_2)$ must
converge to a stationary solution as $t\to\infty$ \cite{hirsch}.
It can be concluded that the corresponding solution
$(x_1,x_2,x_3,x_4)$ tends to a stationary solution. Information
about stationary points of the system can be obtained by stability
considerations, as was mentioned following Theorem \ref{thm3}. This
discussion is summed up in the following theorem.
\begin{theorem} \label{thm4}
Any solution $x_i$ $(i=1,2,3,4)$ of \eqref{abc1}-\eqref{abc4}
with parameter values satisfying $C\le\min\{A,B\}$ converges to a
stationary solution as $t\to\infty$. For fixed values of $A$
and $B$ and $C$ sufficiently small the system has at most three
stationary solutions and at most two stable stationary solutions.
\end{theorem}
The results which have been obtained are unfortunately not
sufficient to give a rigorous confirmation of the qualitative
behaviour shown in Fig. 2a and Fig. 2b of \cite{levbaror}. The
condition $C<\min\{A,B\}$ is not satisfied by the parameter values
in Fig. 2a. It is satisfied in the case of Fig. 2b of
\cite{levbaror} but no information is obtained about the number of
solutions.
\section{Dynamics in the absence of macrophages}\label{Tcellonly}
In this section it will be shown that all the dynamical behaviour
which has been shown to occur in system \eqref{basic} also occurs
in a truncated model where the influence of macrophages is
ignored. The system has two unknowns $x_1$ and $x_2$ with the same
interpretation as before. The effect of the macrophages is turned
off by setting $a_{ij}=0$ when $i$ is one or two and $j$ is three
or four. This leads to the following closed system for $x_1$ and
$x_2$:
\begin{gather}
\frac{dx_1}{dt}=-d_1x_1+g(a_{11}x_1-a_{22}x_2) \\
\frac{dx_2}{dt}=-d_2x_2+g(-a_{11}x_1+a_{22}x_2)
\end{gather}
All the coefficients $d_1$, $d_2$, $a_{11}$ and $a_{22}$ are
assumed positive, as before. Restricting further to the case that
$d_1=d_2=d$ and $\theta=0$ allows the analysis of the dynamics to
be reduced to that of a scalar equation as in the previous
section. The scalar equation is
\begin{equation}
\frac{dx_1}{dt}=-dx_1+g((a_{11}+a_{22})x_1-d^{-1}a_{22}).
\end{equation}
This is essentially the equation \eqref{toy} analysed before and
so all the phenomena found previously occur here also.
\section{Further remarks}
In this paper it has been possible to give a rigorous analysis of
the asymptotics of solutions of the system of Lev Bar-Or
\cite{levbaror} for some values of the parameters. Unfortunately
there are large ranges of the parameters for which no conclusions
were obtained or for which those conclusions are incomplete. The
latter statement even includes the two cases for which numerical
plots were included in \cite{levbaror}. The following questions
have not been answered for the system \eqref{basic} with general
parameters.
\begin{itemize}
\item Are there periodic solutions?
\item Are there homoclinic orbits?
\item Are there heteroclinic cycles?
\item Are there strange attractors?
\end{itemize}
It should be emphasized that there is no evidence, analytical or
numerical, that the answer to any of these questions is yes in the
general case. We are left with a picture of a dynamical system
where the long-time behaviour seems to be simple but it is quite
unclear how to prove it except for a restricted set of values of
the parameters.
It is consistent with everything which has been found here that
the system of four equations with eighteen parameters produces no
qualitatively new phenomena in comparison with a reduced
two-dimensional system with three parameters. No effects were
found which are specifically dependent on the inclusion of the
presentation of antigen by macrophages. The interactions between
Th1 and Th2 cells appear to be sufficient. It is seen that under
some circumstances the model with T cells and cytokines alone
predicts a situation of bistability where either a Th1 or Th2
dominated state can be approached, depending on the initial data.
\begin{thebibliography}{00}
\bibitem{deleenheer03} De Leenheer, P. and Smith, H. L.;
\emph{Virus dynamics: a global analysis}. SIAM J. Appl. Math.
63 (2003), 1313--1327.
\bibitem{egami09} Egami, C.;
\emph{Bifurcation analysis of the Nowak-Bangham
model in CTL dynamics}. Math. Biosci. 221(2009), 33--42.
\bibitem{hale} Hale, J. K.;
\emph{Ordinary differential equations}. Wiley
Interscience, New York, 1969.
\bibitem{hirsch} Hirsch, M. W.;
\emph{Systems of differential equations which
are competitive or cooperative. I: limit sets}. SIAM J. Math.
Anal. 13(1982), 167--179.
\bibitem{korobeinikov04} Korobeinikov, A.;
\emph{Global properties of basic
virus dynamic models}. Bull. Math. Biol. 66 (2004), 879--883.
\bibitem{kuznetsov} Kuznetsov, Y. A.;
\emph{Elements of applied bifurcation theory}.
Springer, Berlin, 1995.
\bibitem{lassmann} Lassmann, H. and Ransohoff, R. M.;
\emph{The CD4-Th1 model for multiple sclerosis: a crucial
re-appraisal}. Trends in Immunology 25 (2004), 132--137.
\bibitem{levbaror} Lev Bar-Or, R.;
\emph{Feedback mechanisms between T helper
cells and macrophages in the determination of the immune response}.
Math. Biosciences 163 (2000), 35--58.
\bibitem{marino} Marino, S. and Kirschner, D. E.;
\emph{The human immune reponse to Mycobacterium tuberculosis
in lung and lymph node}. J. Theor. Biol. 227 (2004),
463--486.
\bibitem{murphy} Murphy, K. M., Travers, P. and Walport, M.;
\emph{Janeway's immunobiology}. Garland Science, New York, 2007.
\bibitem{nowak} Nowak, M. A. and May, R.;
\emph{Virus dynamics}. Oxford University
Press, Oxford, 2000.
\bibitem{oldstone} Oldstone, M. B. A.;
\emph{Molecular mimcry and autoimmune
disease}. Cell 50 (1987), 819--820.
\bibitem{perelson99} Perelson, A. S. and Nelson, P. W.;
\emph{Mathematical
analysis of HIV-1 dynamics in vivo}. SIAM Rev. 41 (1999), 3-44.
\bibitem{smith86} Smith, H. L.;
\emph{Competing subcommunities of mutualists and
a generalized Kamke theorem}. SIAM J. Appl. Math. 46 (1986), 856--874.
\bibitem{wigginton} Wigginton, J. E. and Kirschner, D.;
\emph{A model to predict cell-mediated immune regulatory
mechanisms during human infection with Mycobacterium tuberculosis}
J. Immunol. 166 (2001), 1951--1967.
\end{thebibliography}
\end{document}