\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2011 (2011), No. 142, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2011 Texas State University - San Marcos.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2011/142\hfil Recurrent epidemics]
{Recurrent epidemics resulting from secondary transmission
routes in SIR models}
\author[M. R. Adams, S. Obara\hfil EJDE-2011/142\hfilneg]
{Malcolm R. Adams, Samuel Obara} % in alphabetical order
\address{Malcolm R. Adams \newline
Department of Mathematics, University of Georgia,
Athens, GA 30602, USA}
\email{mradams@uga.edu, Fax 706-542-2573}
\address{Samuel Obara \newline
Department of Mathematics, Texas State University,
San Marcos, TX 78666, USA}
\email{so16@txstate.edu}
\thanks{Submitted March 9, 2011. Published October 28, 2011.}
\subjclass[2000]{34C15, 34C23, 92D30}
\keywords{Epidemiology; recurrence; Hopf bifurcation; SIR model}
\begin{abstract}
In this article, we analyze the behavior of solutions to
a variant of the SIR (susceptible, infected, recovered)
model from epidemiology. The model studied includes a secondary
route for susceptible individuals to be exposed to the infectious
agent. This secondary route provides a feedback mechanism that,
within certain parameter regimes, allows for a limit cycle; i.e.,
sustained periodic behavior in the solutions.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks
\section{Introduction}
The SIR model in epidemiology gives a simple dynamic description
of three interacting populations, the Susceptibles, the Infected,
and the Recovered. In spite of its simplicity, the SIR model,
exhibits the basic structure generally associated to the spread
of a disease in a population: after a possible initial epidemic,
the infected population either tapers to zero or to a stable
endemic level. Many variations of the SIR model have been studied
in recent years to more accurately model more complex diseases and
infection mechanisms. For instance the susceptible population may
be divided into subgroups with different infection rates
\cite{AM1, Esteva, Macdonald, Zhien}, the disease may affect
reproductive rates in the infected or recovered population
\cite{Anderson}, or there may be multiple levels of infections,
some lethal, some sublethal \cite{Boots, Sait}. These more
complicated models can lead to instabilities, and even cyclic
behavior in the infected population.
Here we study a simple variation of the SIR model
introduced in \cite{Fiorello} to model the dynamics of
distemper and parvo virus in jaguars of the Bolivian jungle.
In this model, the disease can infect a susceptible individual
not only through contact with an infected individual, but also
through a secondary method, in this case through contact with
infected feces. Models with a similar feedback mechanism were
studied in \cite{Ewald} and later in \cite{Thieme} in the context
of the evolution of virulence in waterborne diseases.
Numerical evidence in \cite{Fiorello} showed that this model
can support enduring cyclic epidemics.
Heuristically, this can be explained by the idea that the
delayed feedback of the fecal infections can cause a recurrence
of the epidemic. Of course, infection rates and fecal decay
rates must be right for such recurrence to endure.
There are many diseases that bring great economic and human cost
to society which accommodate multiple infection routes. Other
than through direct contact, diseases are often transmitted
through animal or insect vectors (e.g. malaria or plague), through
fecal contamination of water or food supplies (cholera, typhus,
avian influenza), or through ingestion of infected tissues
(brucellosis, bovine spongiform encephalopathy). Detailed
modelling of these infection routes can be quite complicated,
involving such issues as the life cycle of the vector, or the
diffusion of infected materials through water sources. The model
we study does not attempt to understand the details of such
mechanisms, but, as with the SIR model, is meant as a broad
conceptual tool for giving rudimentary insight into the general
behavior of the dynamics of such diseases.
In this article we will provide an analytic proof that an SIR
type model with two transmission routes does support enduring
cyclic epidemics. We will outline the proof that a Hopf bifurcation
gives rise to this cyclic behavior. We will also provide a
numerical example showing that the cyclic behavior can persist
in the full model studied \cite{Fiorello, Obara}.
\section{Modelling}
We begin with a basic SIR model with logistic term:
\begin{equation} \label{eq1}
\begin{gathered}
\dot S = bN(1-\frac{N}{K}) - \beta IS \\
\dot I = \beta IS - (d+r+d_i)I \\
\dot R = rI - dR
\end{gathered}
\end{equation}
where $S$ denotes the susceptible population, $I$ the infected
population, $R$ the recovered population, and $N= S+I+R$ is the
total population. The parameters are $b$, the growth rate; $K$,
the carrying capacity; $\beta$, the infection rate; $d$, the
natural death rate; $r$, the recovery rate; and $d_i$, the
additional death rate due to the disease.
Most analytic discussions of the SIR model treat only the case of
infinite carrying capacity ($K = \infty$) in order to make the
algebra simpler. The case of finite carrying capacity is much more
realistic in that it restricts the total population to remain
bounded. Indeed, if the total population is ever larger than (or
equal to) $K$, then $dN/dt = dS/dt +dI/dt +dR/dt < 0$.
It should be noted that there are other ways to introduce a
logistic term in SIR models. Indeed the model we have chosen
might be rejected since if $S= 0$ and $N>K$ this model yields
$dS/dt <0$, (and thus $S(t)$ would become negative). However, the
region described by $S>0$, $I>0$, $R >0$, and $N (d+r+d_i)/\beta > K$, but for
all initial values, the infected population eventually tends to
zero.)
(2) When $\beta K > d +d_i +r$ the equilibrium $(K,0,0)$ becomes
unstable (a saddle) and there is a third equilibrium in the
positive octant, which we denote $(S_0,I_0,R_0)$. This
equilibrium is always stably attracting, representing an endemic
population of infected individuals. Thus, after an initial
epidemic, solution curves approache $(S_0,I_0,R_0)$. The limiting
behavior can be a damped oscillation (see Figure \ref{fig1}).
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1}
\end{center}
\caption{A damped oscillation limiting to an endemic infected
population for the SIR model with a logistic term.
Parameter values are
$\beta =.12$, $b = .20$, $K = 5$, $r = .2$, $d = .05$
$d1 = 0$. The blue curve is $S(t)$, the green one is $I(t)$,
and the red one is $R(t)$} \label{fig1}
\end{figure}
In order to include in this model the possibility of disease
transmission through a route other than direct contact we
introduce a new dependent variable $F(t)$ representing a
quantity which can infect susceptibles through mutual
contact at a rate given by the parameter $\gamma$.
This secondary transmission factor is produced by the infected
population $I$ at a relative rate given by $\delta$ and it
leaves the environment at a relative rate given by a parameter $\alpha$. Thus, we study the four dimensional system given by
\begin{gather*}
\dot S = bN(1-\frac{N}{K}) - \beta IS -\gamma FS \\
\dot I = \beta IS +\gamma FS - (d+r+d_i)I \\
\dot R = rI - dR \\
\dot F = \delta I - \alpha F
\end{gather*}
Although this model is fairly simplistic, it seems to be too
complicated to obtain complete analytic results. Instead,
we will study a lower dimensional reduction of this system for
which we are able to obtain a complete analysis. At the end of
this section we will show numerically that some of the salient
features of our reduced system are preserved in the full
four - dimensional system.
The simplifying assumption that we use to reduce this system to
an analytically amenable case is to decouple the $R$ dependence
from the rest of the system. This could happen by a variety of
mechanisms. For instance, if the disease is fatal, so $r=0$, then,
with $R(0) = 0$, we have that $R(t) = 0$ for all $t$.
Another mechanism for decoupling $R$ would be replacing $N$ in the
logistic term by $S+I$. This would be appropriate if the recovered
population does not contribute to the reproductive role represented
by the logistic term (e.g. the recovereds are sterile) and nor do
they contribute to the limiting nature of the carrying capacity
(this might be reasonable if $R$ remains relatively small).
After assuming this decoupling we end up with the three dimensional
system given by
\begin{gather*}
\dot S = b(S+I)(1-\frac{S+I}{K}) - \beta IS -\gamma FS \\
\dot I = \beta IS +\gamma FS - (d+r+d_i)I \\
\dot F = \delta I - \alpha F.
\end{gather*}
In the next section we will provide an analytic proof of the
following results.
(1) Solutions remain bounded for all time.
(2) If $K(\beta \alpha +\gamma \delta) < d\alpha$ then there
are only two biologically relevant equilibria (in the positive octant)
located at $(0,0,0)$ and $(K,0,0)$. The origin is a saddle point
and $(K,0,0)$ is an attractor. So, under these conditions,
the long term behavior predicts that the infected population
eventually approaches zero (after a possible epidemic) and
the susceptible population will approach the carrying capacity, $K$.
(3) On the other hand, when
$K(\beta \alpha +\gamma \delta) > d\alpha$ a third equilibrium moves
into the positive quadrant and so becomes biologically relevant.
At the same time the equilibrium at $(K,0,0)$ becomes unstable
(a saddle). The third equilibrium can be stable or unstable
according to the parameter values. To obtain a glimpse of the
nature of this equilibrium we make a change of variable in
the parameter space and then study as a bifurcation parameter
the quantity $c_4 = \frac{\delta \gamma K}{\alpha b}$.
When $c_4 = 0$ the third equilibrium is a sink. But it can be
shown that there is a region of parameter values for which there
is a value $c_4^*$ at which a linear Hopf bifurcation occurs.
That is, there is some $\epsilon >0$ such that if $c_4$ satisfies
$c_4^* - \epsilon < c_4 < c_4^*$ then the third equilibrium is a
stable spiral point, but when $c_4^* K$,
so $N$ (and hence each of $S$ and $I$) must remain
bounded. Finally, if $I$ is bounded, say by $M$, we have that
$\dot F$ becomes negative if ever $F$ grows beyond $\delta M/\alpha$.
Hence $F$ cannot grow without bound either. We also remark here
that since $\dot N <0$ when $N = K$, it follows that the region
given by $S\ge0, I\ge0$, and $N\le K$ is invariant, and so we need
not worry about negative population components as long as our
initial value $N(0)$ is less than $K$.
Now we move on to the equilibrium analysis of our equations.
We begin by non-dimensionalizing these equations by taking
$x=S/K, y = I/K, z= \gamma F/b$ and rescaling time by $\tau = bt$.
This yields
\begin{gather*}
\dot x =(x+y)(1-x-y) -\frac{\beta K}{b}xy -xz\\
\dot y = \frac{\beta K}{b} xy +xz -\frac{d}{b}y \\
\dot z = \frac{\delta \gamma K}{b^2} y - \frac{\alpha}{b} z.
\end{gather*}
Setting $c_1 = \frac{\beta K}{b}$, $c_2 = \frac{d}{b}$,
$c_3 = \frac{\alpha}{b}$ and $c_4 = \frac{\delta \gamma K}{\alpha b}$
we obtain
\begin{gather*}
\dot x =(x+y)(1-x-y) -c_1 xy -xz\\
\dot y =c_1 xy +xz -c_2 y \\
\dot z =c_3(c_4 y - z).
\end{gather*}
(Now, of course, the region of interest is given by $x\ge 0$,
$y\ge 0$, and $x+y \le 1$.)
At an equilibrium point, we must have
\begin{gather*}
(x+y)(1-x-y) -c_1xy -xz =0 \\
c_1xy +xz - c_2 y =0 \\
c_4y-z =0.
\end{gather*}
Substituting $z = c_4y$ into the first two equations yields:
\begin{gather*}
(x+y)(1-x-y)- (c_1+c_4)xy = 0 \\
y\big( (c_1+c_4)x - c_2\big) = 0.
\end{gather*}
There are two equilibria with $y=0$:
$ (0,0,0)$ and $(1,0,0)$.
When $y\ne 0$ there are two others, $(x_0, y_{0+}, c_4y_{0+})$
and $(x_0, y_{0-}, c_4y_{0-})$, where $x_0 = \frac{c_2}{c_1+c_4}$ and
$$
y_{0\pm} = \frac{1}{2}(1-2x_0 -c_2) \pm\frac{1}{2}
\sqrt{(2x_0 +c_2 -1)^2-4x_0(x_0-1)}.
$$
Since $(2x_0 +c_2-1)^2 -4x_0(x_0-1) = (c_2-1)^2 +c_2x_0$,
it follows that the two roots, $y_{0\pm}$, are always real.
Furthermore, if $x_0 <1$, then $y_{0+} >0$ and $y_{0-}<0$
so only $y_{0+}$ is biologically relevant.
Indeed, a straightforward computation shows that if $x_0 <1$,
then also $x_0 +y_0 <1$, so this equilibria is in the relevant
region for our model. For this reason we will henceforth refer
to $y_{0+}$ as simply $y_0$. (Finally, note that if $x_0 >1$,
then both $y_{0+}$ and $y_{0-}$ are negative, and so not
biologically relevant.)
We note that if $c_4 = 0$ then $z$ decreases to zero so we limit
to the 2 dimensional system
\begin{gather*}
\dot x =(x+y)(1-x-y) - c_1 xy \\
\dot y = c_1 xy -c_2 y
\end{gather*}
which is a standard SIR model with a logistic term
(and $R$ is removed). Thus we think of $c_4$ as a perturbation
parameter and analyze bifurcations of the system as $c_4$ varies,
holding the other parameters constant.
The Jacobi matrix for system \eqref{eq2} is
$$
J = \begin{pmatrix} 1-2x-2y-c_1y-z & 1-2x-2y -c_1x & -x \\
c_1y +z & c_1x-c_2 & x \\ 0 & c_3c_4 & -c_3 \end{pmatrix}.
$$
At the point $(0,0,0)$ this gives
$$
J_0 = \begin{pmatrix} 1 & 1 & 0 \\ 0 & - c_2 &0 \\
0 & c_3c_4 & -c_3 \end{pmatrix}.
$$
Since the eigenvalues are $1, -c_2$, and $-c_3$, the point $(0,0,0)$
is a saddle.
At the point $(1,0,0)$, we obtain
$$
J_1 = \begin{pmatrix}- 1 & -1 -c_1 & -1 \\ 0 & c_1-c_2 & 1 \\
0 & c_3c_4 & -c_3 \end{pmatrix}.
$$
The eigenvalues of this matrix are $ -1$, and
\begin{align*}
&\frac{1}{2}\Big(c_1-c_2-c_3 \pm \sqrt{(c_1-c_2-c_3)^2
+4c_3( c_1-c_2+c_4)}\Big)\\
&= \frac{1}{2}\Big(c_1-c_2-c_3 \pm \sqrt{(c_1-c_2+c_3)^2 +4c_3c_4}\Big).
\end{align*}
The above expression shows that all eigenvalues are always real.
The signs of the last two eigenvalues depend on the sign of $x_0-1$.
When $x_0 >1$ we have $c_2 > c_1 +c_4$ and hence the terms
$c_1-c_2-c_3$ and $c_1-c_2 +c_4$ are both negative.
This implies that all three eigenvalues are negative; i.e.,
this equilibrium is a sink. When $x_0 <1$, we have $c_2< c_1+c_4$
so the term $c_1-c_2 +c_4$ is positive. This implies that one
of the eigenvalues is positive while the other two are negative;
the equilibrium is a saddle.
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.49\textwidth]{fig2a}
\includegraphics[width=0.49\textwidth]{fig2b}
\end{center}
\caption{In the left figure, the parameter values are
$c_1 = .2$, $c_2 = .5$, $c_3 = 1$, and $c_4 = .2$, so $x_0 = 1.25 >1$.
Thus all solutions limit to $x=1$ (blue), $y=0$ (green),
and $z = 0$ (red). In the right hand figure, the parameter values
are $c_1 = 4, c_2 = 2, c_3 = .01$, and $c_4 = 5$, so $x_0 = 2/9 <1$.
Here we see the trend to a stable endemic diseased population}
\label{fig2}
\end{figure}
Summarizing our results so far, we have determined the following
(see Figure \ref{fig2}):
(1) If $x_0 >1$ then there are only two biologically relevant
equilibria: A saddle at $(0,0,0)$ and a sink at $(1,0,0)$.
(2) If $x_0 <1$ then there are three biologically relevant
equilibria: Saddles at $(0,0,0)$ and $(1,0,0)$.
We study the nature of the third equilibrium $(x_0,y_0,c_4y_0)$ next.
In the following discussion we will assume that $x_0<1$ so that
the equilibrium $(x_0,y_0,c_4y_0)$ is biologically relevant.
Here the Jacobian becomes
$$
J_2 = \begin{pmatrix} 1- 2x_0 -2y_0 - (c_1 +c_4)y_0
& 1-2x_0-2y_0-c_1x_0 & -x_0 \\ (c_1 +c_4)y_0 & c_1x_0 - c_2 & x_0
\\ 0 & c_3c_4 & -c_3 \end{pmatrix},
$$
where $x_0 = c_2/(c_1 +c_4) <1$.
Here it is useful to make a change of variable in the parameter space:
we replace the parameter $c_1$ by $x_0$ using the relation
$c_1 = (c_2/x_0) -c_4$. Substituting this into $J_2$ yields
$$
J_2 = \begin{pmatrix}1- 2(x_0+y_0) - \frac{c_2y_0}{x_0}
&1- 2(x_0+y_0)+ c_4x_0 -c_2 & -x_0 \\
c_2y_0/x_0& -c_4 x_0 & x_0 \\ 0 & c_3c_4 & -c_3
\end{pmatrix}.
$$
The term $1- 2(x_0+y_0)$ can be simplified to
$c_2-\sqrt{(c_2-1)^2 +4c_2x_0}$, which, for the moment,
we will simply denote by $\omega$. Note that $\omega -c_2 <0$.
Also, from $(x_0 +y_0)(1-x_0-y_0) = (c_1 +c_4)x_0y_0 = c_2y_0$,
so we have $1-x_0 - y_0 = \frac{c_2y_0}{x_0 + y_0}$ and so
$\omega = \frac{c_2y_0}{x_0 + y_0} - (x_0 + y_0)$.
Then, a computation shows that the (negative of the) characteristic
polynomial of $J_2$ is
\begin{align*}
P(\lambda) &= \lambda ^3 + \big(c_3 +\frac{c_2y_0}{x_0}-\omega \big)
\lambda^2 + \big(\frac {c_2y_0}{x_0}(c_3+c_2 -\omega) -c_3 \omega
\big)\lambda + c_3(c_2-\omega)\frac{c_2y_0}{x_0} \\
&\quad + c_4 x_0\big( \lambda (\lambda - \omega) \big).
\end{align*}
When $c_4 = 0$, the roots are $-c_3$ and
$$
\frac {1}{2} \left(\omega - \frac{c_2y_0}{x_0}
\pm \sqrt{(\omega -\frac{c_2y_0}{x_0})^2 +4(\frac{c_2y_0}{x_0}
(\omega -c_2)})\right).
$$
Now, as noted above, $\omega - c_2 <0$, so the real part of the
last two roots of the characteristic polynomial have the same
sign as $\omega - \frac{c_2y_0}{x_0}$. But, again
from above, $\omega = \frac{c_2y_0}{y_0+x_0} -(x_0+y_0)$,
so
$$
w-\frac{c_2y_0}{x_0} = c_2y_0(\frac{1}{x_0+y_0}
- \frac{1}{x_0} ) - (x_0 +y_0) <0.
$$
Thus, this equilibrium is a sink when $c_4 = 0$.
To analyze the general case we use the following basic results
about a general cubic polynomial of the form
$P(\lambda) = \lambda^3 + A\lambda^2 + B\lambda +C$:
\begin{lemma} \label{lem1}
Let $\Delta = A^2B^2 + 18 ABC -27C^2 - 4B^3 - 4A^3C$, then
$P(\lambda) = 0$ has two complex roots if and only if
$\Delta < 0$. Moreover, if this is the case, then the sign of
the real parts of these complex roots is the same as $C - AB$.
\end{lemma}
Here we set
\begin{gather*}
A_0 = \Big( c_3+\frac{c_2y_0}{x_0} -\omega\Big), \\
B_0 = \Big(c_3(\frac{c_2y_0}{x_0} -\omega),
+\frac{c_2y_0}{x_0}(c_2 - \omega)\Big), \\
C_0 = c_3 (c_2- \omega)\frac{c_2y_0}{x_0}, \\
A_1 = x_0, \quad
B_1 = -x_0\omega .
\end{gather*}
(Note that $A_0$, $B_0$, and $C_0$ are all positive.)
Also, $A = A_0 + c_4A_1$,
$B= B_0 + c_4 B_1$, and $C = C_0$. So
\begin{equation}
AB -C = A_0B_0 - C_0 +c_4(A_0 B_1 +A_1 B_0) + c_4^2A_1B_1.
\label{quadrat}
\end{equation}
Note that
$$
A_0B_0 - C_0 = (\frac{c_2y_0}{x_0}-\omega)
\Big(c_3^2+c_3(\frac{c_2y_0}{x_0} -\omega)
+\frac{c_2y_0}{x_0}(c_2-\omega)\Big)
$$
which is greater than zero since $\frac{c_2y_0}{s_0} -\omega >0$
and $c_2-\omega >0$. Since (\ref{quadrat}) is quadratic in $c_4$
and positive at $c_4 = 0$ we see that $AB-C$ becomes negative
after some critical value, $c_4^*$, of $c_4$ if and only if the
coefficient, $A_1B_1$, of $c_4^2$ is negative.
Since $A_1B_1 =- x_0^2\omega$ we see that this is negative only
if $\omega >0$; i.e.,
\begin{equation}
c_2> \sqrt{(c_2-1)^2 +4c_2x_0}. \label{negcond}
\end{equation}
Squaring both sides, and simplifying, yields that inequality
\eqref{negcond} holds exactly when
$$
x_0< \frac{1}{2}\big(1-\frac{1}{2c_2}\big).
$$
Thus we have proven the following theorem.
\begin{theorem} \label{thm1}
Assume that $J_2$ has two complex eigenvalues and let
$A_0$, $B_0$, $C_0$, $A_1$, $B_1$ be defined as above.
Let $p(c_4)= A_0B_0 - C_0 +c_4(A_0 B_1 +A_1 B_0) +c_4^2A_1B_1$,
and assume that $x_0 < \frac{1}{2}(1-\frac{1}{2c_2})$.
Then $p$ has two real roots, one positive and one negative.
Letting $c_4^*$ denote the positive root, we have that for
$c_4 > c_4^*$ the sign of the real part of the complex eigenvalues
of $J_2$ will be positive.
\end{theorem}
Next we claim that $J_2$ has complex eigenvalues when $c_4$ is
in a neighborhood of $c_4^*$. When $c_4= c_4^*$ we have that
\begin{equation}
C(c_4^*) = A(c_4^*)B(c_4^*) \label{ABC}
\end{equation}
so
\begin{align*}
-\Delta(c_4^*)
& = 8A(c_4^*)^2B(c_4^*)^2 +4B(c_4^*)^3 +4A(c_4^*)^4B(c_4^*) \\
&= 8A(c_4^*)^2B(c_4^*)^2 +4B(c_4^*)(B(c_4^*)^2 +A(c_4^*)^4).
\end{align*}
Notice that $B(c_4^*) >0$ from \eqref{ABC} since $A(c_4)$
and $C(c_4)$ are positive for all values of $c_4$.
Thus we have proven the following result.
\begin{theorem} \label{thm2}
Assume that $x_0 <\frac{1}{2}(1-\frac{1}{2c_2})$, and let
$c_4^*$ be defined as above. Then the system has a linear
Hopf bifurcation at $c_4^*$; i.e., there is a $\delta >0$ such
that for $c_4$ satisfying $c_4^*-\delta c_4^*$, we have an outward spiral from $(x_0, y_0,z_0)$
which is being compressed in the complementary direction and
remains bounded. This certainly indicate a likely limit cycle.
(See Figure \ref{fig3}.)
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3}
\end{center}
\caption{Oscillatory behavior for the system studied in this section.
The parameter values are $x_0=.05$, $c_2 = 5$, $c_3 =1.2$, and
$ c_4 = 99$ (so $c_1 = 1)$. The critical value is $c_4^* = 95.4$.
In the graph, the values of $z(t)$ (red) are divided by 10 to
make the scaling compatible} \label{fig3}
\end{figure}
\section{The full four dimensional model}
Finally, we return to the full 4 dimensional model.
If the reduced model is robust, we would expect that for small
values of $r$, we will see the same behavior that we see when $r= 0$.
To verify that this is the case, we look at two sets of parameter
values, the first satisfying $c_4 < c_4*$ and the second with
$c_4 > c_4*$, and numerically compute the eigenvalues of the
linearized system near the equilibrium. Before proceeding,
we should first repeat the dimensional reduction for this system.
We begin with the system
\begin{gather*}
\dot S = bN(1-\frac{N}{K}) - \beta IS -\gamma FS \\
\dot I = \beta IS +\gamma FS - (d+r+d_i)I \\
\dot R = rI - dR \\
\dot F = \delta I - \alpha F.
\end{gather*}
As above, we rescale time by $\tau = bt$ and set
$x= S/K$, $y = I/K$, $z=\gamma F/b$, as well as $w= R/K$.
This yields the equations
\begin{gather*}
\dot x = (x+y+w)(1-x-y-w) -c_1 xy -xz\\
\dot y = c_1 xy +xz -(c_2 +c_5 +c+6)y \\
\dot z = c_3(c_4 y - z)\\
\dot w = c_5y - c_2w
\end{gather*}
with $c_1 = \frac{\beta K}{b}$, $c_2 = \frac{d}{b}$,
$c_3 = \frac{\alpha}{b}$, $c_4 = \frac{\delta \gamma K}{\alpha b}$,
and $c_5 = r/b$, $c_6 = d_i/b$.
This system has an equilibrium $(x_0, y_0, z_0, w_0)$, where
$x_0 = (c_2+c_5+c_6)/(c_1 +c_4)$,
$z_0 = c_4y_0, w_0 = \frac{c_5}{c_2}y_0$, and
$ y_0$ is the positive root of the polynomial equation
$$
\Big(\frac{c_2+c_5}{c_2}\Big)^2y^2
+ \Big(c_2+c_5+c_6 + \big(\frac{c_2+c_5}{c_2}\big)(2x_0 -1)\Big)y
+x_0(x_0-1) = 0.
$$
We study this system near the parameter values given in
Figure \ref{fig3} above. We will keep the recovery rate, $r$,
small so that we are near to the 3 dimensional system studied above.
Also, to keep $x_0$ fixed as in the above analysis, we will need
to keep $c_1 +c_4$ constant. In Figure \ref{fig4} we show graphs of $x(t)$
and $y(t)$ with two nearby sets of parameter values.
The graph on the left shows persistent periodic behavior as
in Figure \ref{fig3}, while the graph on the right shows damped
periodic behavior. For the second graph, we have lowered the
value of $c_4$ below $c_4^*$, while raising $c_1$ to keep $x_0$
constant.
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.49\textwidth]{fig4a}
\includegraphics[width=0.49\textwidth]{fig4b}
\end{center}
\caption{Graphs of $x(t)$ and $y(t)$ for the full four dimensional
model. In the left, the parameter values are
$c_1 = 1$, $c_2 = 1$, $c_3 = 1.2$, $c_4 = 99$, $c_5 = 0.01$, and
$c_6= 4$. In the right, the parameter values are
$c_1 = 10$, $c_2 = 1$, $c_3 = 1.2$, $c_4 = 90$, $c_5 = 0.01$, and
$c_6 = 4$. Here we see the trend to a stable endemic diseased
population} \label{fig4}
\end{figure}
The Jacobian for the four dimensional system is
$$
J = \left(\begin{smallmatrix} 1-2x-2y-2w-c_1y-z & 1-2x-2y-2w -c_1x & -x &
1-2x-2y-2w\\ c_1y +z & c_1x-(c_2+c_5+c_6) & x & 0\\
0 & c_3c_4 & -c_3 & 0\\ 0& c_5 &0 & -c_2
\end{smallmatrix}\right).
$$
Using the parameter values $c_1 = 1$, $c_2 = 1$, $c_3 = 1.2$,
$c_4 = 99$, $c_5 = 0.01$, and $c_6= 4$, one can calculate that
the equilibrium is at the point
$(x_0, y_0, z_0, w_0) = (0.0501, 0.0116, 1.1455, 0.0001)$.
Evaluating the Jacobian at this point yields
$$
J1 = \begin{pmatrix} -0.2806 & 0.8263 & -0.0501& 0.8764\\
1.1571& -4.9599& 0.0501& 0\\ 0 & 118.8& -1.2 & 0\\ 0& 0.01 &0 & -1
\end{pmatrix},
$$
whose eigenvalues are
$\lambda_1 = 0.00878 + 0.9418\sqrt{-1}$,
$\lambda_2 = 0.00878 - 0.9418\sqrt{-1}$,
$\lambda_3 = -6.4583$, and $\lambda_4 = -0.9998$.
The positive real parts of the complex eigenvalues explains
the persistent periodic behavior at these parameter values.
Similarly, if we compute the Jacobian at the values
$c_1 = 10$, $c_2 = 1$, $c_3 = 1.2$, $c_4 = 90$, $c_5 = 0.01$, and
$c_6= 4$, one can calculate that the equilibrium is at the point
$(x_0, y_0, z_0, w_0) = (0.0501, 0.0116, 1.1455, 0.0001)$.
Evaluating the Jacobian at this point yields
$$
J2 = \begin{pmatrix} -0.2806 & 0.3754 & -0.0501& 0.8764\\
1.1571& -4.509& 0.0501& 0\\ 0 & 108.8& -1.2 & 0\\ 0& 0.01 &0 & -1
\end{pmatrix},
$$
whose eigenvalues are
$\lambda_1 = -0.0174 + 0.9806\sqrt{-1}$,
$\lambda_2 = -0.0174 - 0.9806\sqrt{-1}, \lambda_3 = -5.955$, and
$\lambda_4 = -0.9998$. The negative real parts of the complex
eigenvalues show that solutions for these parameter values damp
to the equilibrium point.
\subsection*{Concluding Remarks}
The analytic and numerical results of the previous sections indicate
that a secondary transmission route for an infectious disease
provides a robust feedback mechanism that can give rise to
sustained periodic epidemics. Of course this is a very crude
model which does not take into account detailed modelling of the
secondary transmission route, but it is simple enough that
dependence of solutions on the parameters can be studied
analytically, and so lends insight into how the feedback affects
the general behavior of solutions.
\subsection*{Acknowledgements}
We would like to thank Christine
Fiorello and Andrew Sornborger for several helpful discussions
during the preparation of this paper.
\begin{thebibliography}{00}
\bibitem{AM1} Anderson, R. M.; May, R. M.;
Infectious Diseases of Humans,
Oxford University Press, Oxford, 1991.
\bibitem{Anderson} Anderson, R. M., May, R. M.;
The Population Dynamics of Microparasites and their Invertibrate Hosts.
\emph{Philosophical Transactions of the Royal Society of London,
Series B, Biological Sciences}, \textbf{291} (1981) 451-524.
\bibitem{Boots} Boots, M.; Norman, R.;
Sublethal infection and the population dynamics of host-micropar\-asite
interactions. \emph{J. Animal Ecology} \textbf{69} (2000) 517-524.
\bibitem{Esteva} Esteva, L.; Vargas, C.;
Analysis of a dengue disease transmission model.
\emph{Math. Biosci.} \textbf{150} (1998) 131-151.
\bibitem{Ewald} Ewald, P. W.; De Leo G.;
Alternative transmission modes and the evolution of virulence,
\emph{Adaptive Dynamics of Infectious Diseases: in Pursuit of Virulence
Management}, (Dieckmann, U.; Metz, J.A.J., Sabelis, M.W., Sigmund, K., eds.), 10-25, International Institute for Applied Systems Analysis,
Cambridge University Press, Cambridge 2002
\bibitem{Fiorello} Fiorello, C. V.;
Disease Ecology of Wild and Domestic Carnivores in Bolivia,
\emph{Doctoral Dissertation} Columbia University, 2004.
\bibitem{Macdonald} MacDonald, G.;
The Epidemiology and Control of Malaria, Oxford University Press,
London, 1957.
\bibitem{Obara} Obara, Mathematical Models in Epidemiology.
\emph{Master's Thesis} University of Georgia, 2005.
\bibitem{Sait} Sait, S. M.; Begon, M.; Thompson, D. J.;
The effects of a sublethal baculovirus infection in the Indian
meal moth, Plodia interpunctella.
\emph{J. Animal Ecology} \textbf{63} (1994) 541-550.
\bibitem{Thieme} Thieme, H. R.;
Pathogen competition and coexistence and the evolution of virulence.
\emph{Mathematics for Life Sciences and Medicine.}
(Y. Takeuchi, Y. Iwasa, K. Sato, eds.), 123-153, Springer,
Berlin Heidelberg 2007
\bibitem{Zhien} Zhien, M.; Liu, J.; Li, J.;
Stability analysis for differential infectivity epidemic models.
\emph{Nonlinear Analysis: Real World Applications} \textbf{4}
(2003) 841-856.
\end{thebibliography}
\end{document}