\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}