A Mathematical Model to Analyze Spread of Hemorrhagic Disease in White-Tailed Deer Population

Show more

1. Introduction

Hemorrhagic disease (HD) is a fatal disease of white-tailed deer (Odocoileus virginianus). It is the collective term used for epizootic hemorrhagic disease and bluetongue disease (genus Orbivirus). These diseases have similar symptoms and are frequently grouped together and referred to as HD. Symptoms include hemorrhaging, swelling due to fluid accumulation, sores, ulcers, sloughing of hooves, high fever, and loss of fear of humans [1] [2] . There are three different forms of HD (peracute, acute, and chronic) which dictate how long a deer will survive. Death can occur in as little as one to two weeks [3] . It is possible for a deer to survive, but it is rare. In addition to white-tailed deer, HD can be transmitted to other wild ruminants and domestic animals, most commonly hoof stock, but it rarely causes disease. The infection does not affect humans or non-ruminant animals [1] . The vector that spreads HD is small biting midge (Culicoides Ceratopogondiae). These midges are tiny, blood-sucking flies that are merely pests to humans, but they are the vectors in the spread of the disease in deer and livestock.

In the present work, we build a mathematical model to investigate the dynamics of HD. The amount of literature dedicated to the mathematical modeling of vector-borne diseases is extensive (See for example [4] [5] [6] [7] ). The model by Nobel Prize winner Ronald Ross [4] is at the cornerstone of such models, and he used his model to investigate the spread of malaria. Over four decades later, George Macdonald developed it further [5] . In fact, there have been several extensions to the Ross-Macdonald model. For instance, Lou and Zhou [6] included advection and diffusion terms to take the spatial movements of individuals into account. Reaction-diffusion models have also been used for investigating dynamics of vector-borne diseases such as dengue fever [7] and Zika [8] . Using a deterministic modeling approach, the main objective of the present study is to have a better understanding of the possible effects of deer-midge interactions and deer migrations on HD dynamics in a deer population.

In recent years, more realistic models have been constructed which take into account dispersion time and host movements. A key article is the work by Neubert et al. [9] , which argues that dispersion in Lotka-Volterra predator-prey models is unrealistic as individuals leaving an area (i.e., a patch) immediately appear in another. In nature, an individual requires a finite amount of time to complete a trip from one patch to another or to complete a round trip leaving and returning to the same patch. During this time, the migrating individuals are not interacting with other predators or prey in this patch. Thus, Neubert and his co-authors [9] [10] demonstrate that models that incorporate explicit travel-time are often more stable.

Few models have been constructed to analyze the dynamics of HD in white-tailed deer populations and dairy farms. Park et al. [11] studied these dynamics by first fitting a statistical model to predict HD incidents as a function of seroprevalence (i.e., the number of individuals in a population who tested positive for HD). Then, using ordinary differential equations (ODE), they formulated a mechanistic model to support the theory that there is a correlation between the number of HD cases and the number of deer in a population with the virus. Their study suggests that the maximum number of cases occurs at intermediate levels of this seroprevalence. By constructing a realistic model, we will be able to analyze and simulate the dynamics of HD. A better understanding of HD dynamics gives epidemiologists and biologists the capacity to control and predict future epidemics in white-tailed deer populations. The present work is the first step toward realistic modeling of HD dynamics with a focus on migrating effects of white-tailed deer population.

The rest of this paper is organized as follows. In Section 2 we propose the vector-borne model of HD spread. This model takes into account the migration and immigration of deer from and into a single patch. In Section 3, we study the model through linearization, chain-trick method, and equilibrium analysis. We also calculate the ${R}_{0}$ expression and use it in Section 4 to numerically investigate the effects of the model parameters on outbreaks. Finally, in Section 5, we provide a discussion of results and outline the limitations of this study.

2. The Single-Patch Model

In the attempt to create a mathematical model of HD outbreak in a population of white-tailed deer, we make certain assumptions based on the ecology of deer and midge populations and the characteristics of HD. The deer (host) and midge (vector) populations are divided into susceptible and infected classes. At time t, there are ${D}_{S}\left(t\right)$ susceptible deer, ${D}_{I}\left(t\right)$ infected deer, ${M}_{S}\left(t\right)$ susceptible midges, and ${M}_{I}\left(t\right)$ infected midges. The total deer population at time t is ${D}_{N}\left(t\right)={D}_{S}\left(t\right)+{D}_{I}\left(t\right)$ , and the total midge population is ${M}_{N}\left(t\right)={M}_{S}\left(t\right)+{M}_{I}\left(t\right)$ . Susceptible deer become infected through bites of infected midges; susceptible midges become infected when they feed on the blood of an infected deer. As observed in the wild, deer will migrate (disperse) out of and back into a region (i.e., a patch) due to seasonal variations, availability of food, or predators; midges, however, will not. They are weak fliers and typically disperse no more than about a mile from the site of larval development, with females flying farther than males [12] . Moreover, their flying activity is greatly reduced in windy conditions. They may fly as far as six miles or more, but this is very rare [13] . We therefore consider the following assumptions in the model construction:

1) All newborns are susceptible in both populations of deer and midges (i.e., no inherited infection or vertical transmission is considered).

2) Susceptible deer become infected only by adequate contact with infected midges and cannot become infected via contact with an infected deer.

3) Once infected, a deer will die from the disease. (Note, in actuality, there are cases where a deer survives the infection, but it is rare.)

4) Individuals in both populations will die naturally by both density independent and density dependent factors.

5) By the law of mass action, we assume that infection transmission is proportional to the population densities of deer and midges.

6) Deer will frequently travel out of and into a geographic area (a patch), but midges do not (as the amount of dispersal in midge populations is negligible).

A compartmental diagram of the proposed HD model is seen in Figure 1, and a summary of parameters and variables is given in Table 1. All parameters are assumed to be non-negative. Given the above-mentioned assumptions and the model diagram, the set of delayed differential equations representing the model is given by

$\begin{array}{l}\frac{\text{d}{D}_{S}}{\text{d}t}={\lambda}_{D}{D}_{N}-\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left({\mu}_{D}+\rho +{d}_{S}+{\mu}_{2D}{D}_{N}\right){D}_{S}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z\\ \frac{\text{d}{D}_{I}}{\text{d}t}=\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left({\mu}_{D}+\rho +{\gamma}_{D}+{d}_{I}+{\mu}_{2D}{D}_{N}\right){D}_{I}+{d}_{I}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{I}z}{D}_{I}\left(t-z\right)\text{d}z\\ \frac{\text{d}{M}_{S}}{\text{d}t}={\lambda}_{M}{M}_{N}-\frac{{\beta}_{M}{D}_{I}{M}_{S}}{{D}_{N}}-\left({\mu}_{M}+\sigma +{\mu}_{2M}{M}_{N}\right){M}_{S}\\ \frac{\text{d}{M}_{I}}{\text{d}t}=\frac{{\beta}_{M}{D}_{I}{M}_{S}}{{D}_{N}}-\left({\mu}_{M}+\sigma +{\mu}_{2M}{M}_{N}\right){M}_{I}\end{array}$ (1)

In absence of the disease, population growths of deer and midges are

Table 1. Summary of the variables and parameters used in the delayed HD model (1).

Note: All variables and parameters are non-negative. The specific parameter values used in the analysis will be indicated in Table 2, Section 4.

Figure 1. A compartmental diagram of the HD model (1) with population dispersal. Dashed lines represent the HD transmission between the vector and host. Deer migration into the patch is denoted by $\int}g$ and migration out of the patch is denoted by ${d}_{S}$ and ${d}_{I}$ . See Table 1 for a summary of the parameters and variables.

formulated with logistic growth models. These are the terms that include ${\lambda}_{D}$ , ${\lambda}_{M}$ , ${\mu}_{2D}$ , and ${\mu}_{2M}$ in model (1). Similar to [14] , the carrying capacity for the deer population exists and must be positive. Hence, it is required that

${H}_{1}:{\lambda}_{D}>{\mu}_{D}+\rho +{d}_{S}$ (2)

and

${H}_{2}:{\lambda}_{D}>{\mu}_{D}+\rho +{\gamma}_{D}+{d}_{I}.$ (3)

Also, the carrying capacity for midges exists and is positive. Thus,

${H}_{3}:{\lambda}_{M}>{\mu}_{M}+\sigma .$ (4)

Individual deer immigrate from the patch at a constant per capita rate ( ${d}_{S}$ and ${d}_{I}$ ) and return z units of time after their departure. The integrals in the first two equations of model (1) are distributed delay terms representing the influx of susceptible and infected deer, respectively, from all points in time in the past up to and including the present time [15] . The function $g\left(z\right)$ in the integrals is a probability density function for the time it takes for a deer to disperse given that the deer survives the trip, and $g\left(z\right)\text{d}z$ is the probability that a successfully dispersing deer departing at time t completes the trip between time $t+z$ and $t+z+\text{d}z$ . As $g\left(z\right)$ is a probability density function, it is normalized so that ${\int}_{0}^{\infty}}g\left(z\right)\text{d}z=1$ . The functions ${\text{e}}^{-{\delta}_{S}z}$ and ${\text{e}}^{-{\delta}_{I}z}$ in the integrals are the probabilities of a deer surviving a trip of duration z given constant probabilities per unit of time ${\delta}_{S}$ and ${\delta}_{I}$ for the mortality during travel of susceptible and infected deer, respectively. All deer migrating back into this single patch originated in the patch; in other words, there are no new deer entering the patching that originated from somewhere else. Hence, we are studying a herd of deer concentrated within a patch with the ability of migrating in and out of it.

3. Analysis of the Single-Patch Model

3.1. Linear Stability Analysis

In this section, we provide a formal procedure of linear stability analysis which leads to the characteristic equation and the stability conditions for the equilibrium solutions. Specifically, Disease Free Equilibrium (DFE) (i.e., ${D}_{I}^{*}=0$ and ${M}_{I}^{*}=0$ ) and Endemic Equilibrium (EE) are the constant solutions of model (1). In epidemiology, a stable DFE is always desired whereas a stable EE can be of great concern. The first two equations of model (1) have an integral influx term that may be simplified by the following method. Letting

${f}^{\left[1\right]}\left({D}_{S},{D}_{I},{M}_{S},{M}_{I}\right)={\lambda}_{D}{D}_{N}-\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left({\mu}_{D}+\rho +{d}_{S}+{\mu}_{2D}{D}_{N}\right){D}_{S},$ (5)

we rewrite the first equation as

$\frac{\text{d}{D}_{S}}{\text{d}t}={f}^{\left[1\right]}\left({D}_{S},{D}_{I},{M}_{S},{M}_{I}\right)+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z.$ (6)

Similarly, we rewrite the second equation as

$\frac{\text{d}{D}_{I}}{\text{d}t}={f}^{\left[2\right]}\left({D}_{S},{D}_{I},{M}_{S},{M}_{I}\right)+{d}_{I}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{I}z}{D}_{I}\left(t-z\right)\text{d}z,$ (7)

where

${f}^{\left[2\right]}\left({D}_{S}\mathrm{,}{D}_{I}\mathrm{,}{M}_{S}\mathrm{,}{M}_{I}\right)=\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left({\mu}_{D}+\rho +{\gamma}_{D}+{d}_{I}+{\mu}_{2D}{D}_{N}\right){D}_{I}\mathrm{.}$ (8)

As the bottom two equations of model (1) have no integral term, we let ${f}^{\left[3\right]}$ and ${f}^{\left[4\right]}$ equal the right-hand side of the third and fourth equations in model (1), respectively. Let a solution $\left({D}_{S}\left(t\right)\mathrm{,}{D}_{I}\left(t\right)\mathrm{,}{M}_{S}\left(t\right)\mathrm{,}{M}_{I}\left(t\right)\right)$ of model (1) nearby an equilibrium solution $E=\left({D}_{S}^{*},{D}_{I}^{*},{M}_{S}^{*},{M}_{I}^{*}\right)$ be in the form of

$\begin{array}{l}{D}_{S}\left(t\right)={D}_{S}^{*}+{\stackrel{\u02dc}{D}}_{S}\left(t\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{D}_{I}\left(t\right)={D}_{I}^{*}+{\stackrel{\u02dc}{D}}_{I}\left(t\right),\\ {M}_{S}\left(t\right)={M}_{S}^{*}+{\stackrel{\u02dc}{M}}_{S}\left(t\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{M}_{I}\left(t\right)={M}_{I}^{*}+{\stackrel{\u02dc}{M}}_{I}\left(t\right)\end{array}$ (9)

for some ${\stackrel{\u02dc}{D}}_{S}\left(t\right)$ , ${\stackrel{\u02dc}{D}}_{I}\left(t\right)$ , ${\stackrel{\u02dc}{M}}_{S}\left(t\right)$ , and ${\stackrel{\u02dc}{M}}_{I}\left(t\right)$ . Using the Taylor expansion, we linearize the first equation in model (1) about equilibrium E by substituting (9) into (6) and dropping the nonlinear terms. Thus, the first equation of model (1) is linearized as follows.

$\begin{array}{c}\frac{\text{d}{D}_{S}}{\text{d}t}=\frac{\text{d}{\stackrel{\u02dc}{D}}_{S}}{\text{d}t}={f}^{\left[1\right]}\left({D}_{S}\mathrm{,}{D}_{I}\mathrm{,}{M}_{S}\mathrm{,}{M}_{I}\right)+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z\\ ={f}^{\left[1\right]}\left({D}_{S}^{\mathrm{*}}\mathrm{,}{D}_{I}^{\mathrm{*}}\mathrm{,}{M}_{S}^{\mathrm{*}}\mathrm{,}{M}_{I}^{\mathrm{*}}\right)+\frac{\partial {f}^{\left[1\right]}}{\partial {D}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{S}+\frac{\partial {f}^{\left[1\right]}}{\partial {D}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{I}+\frac{\partial {f}^{\left[1\right]}}{\partial {M}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{S}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\partial {f}^{\left[1\right]}}{\partial {M}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{I}+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}\left({D}_{S}^{\mathrm{*}}+{\stackrel{\u02dc}{D}}_{S}\left(t-z\right)\right)\text{d}z\\ ={f}^{\left[1\right]}\left({D}_{S}^{\mathrm{*}}\mathrm{,}{D}_{I}^{\mathrm{*}}\mathrm{,}{M}_{S}^{\mathrm{*}}\mathrm{,}{M}_{I}^{\mathrm{*}}\right)+\frac{\partial {f}^{\left[1\right]}}{\partial {D}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{S}+\frac{\partial {f}^{\left[1\right]}}{\partial {D}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{I}+\frac{\partial {f}^{\left[1\right]}}{\partial {M}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{S}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\partial {f}^{\left[1\right]}}{\partial {M}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{I}+{d}_{S}{D}_{S}^{\mathrm{*}}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}\text{d}z+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{\stackrel{\u02dc}{D}}_{S}\left(t-z\right)\text{d}z\mathrm{.}\end{array}$ (10)

We know that equilibrium E satisfies the first equation of model (1), hence

${f}^{\left[1\right]}\left({D}_{S}^{*},{D}_{I}^{*},{M}_{S}^{*},{M}_{I}^{*}\right)+{d}_{S}{D}_{S}^{*}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}\text{d}z=0,$ (11)

and thus

${d}_{S}{D}_{S}^{*}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}\text{d}z=-{f}_{1}\left({D}_{S}^{*},{D}_{I}^{*},{M}_{S}^{*},{M}_{I}^{*}\right).$ (12)

Substituting (12) into (10) yields

$\begin{array}{c}\frac{\text{d}{\stackrel{\u02dc}{D}}_{S}}{\text{d}t}=\frac{\partial {f}^{\left[1\right]}}{\partial {D}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{S}+\frac{\partial {f}^{\left[1\right]}}{\partial {D}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{I}+\frac{\partial {f}^{\left[1\right]}}{\partial {M}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{S}+\frac{\partial {f}^{\left[1\right]}}{\partial {M}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{I}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{\stackrel{\u02dc}{D}}_{S}\left(t-z\right)\text{d}z.\end{array}$ (13)

Applying the same procedure to equation (7), we get that the second equation of model (1) is linearized by

$\begin{array}{c}\frac{\text{d}{\stackrel{\u02dc}{D}}_{I}}{\text{d}t}=\frac{\partial {f}^{\left[2\right]}}{\partial {D}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{S}+\frac{\partial {f}^{\left[2\right]}}{\partial {D}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{D}}_{I}+\frac{\partial {f}^{\left[2\right]}}{\partial {M}_{S}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{S}+\frac{\partial {f}^{\left[2\right]}}{\partial {M}_{I}}\left(E\right)\cdot {\stackrel{\u02dc}{M}}_{I}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{d}_{I}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{I}z}{\stackrel{\u02dc}{D}}_{I}\left(t-z\right)\text{d}z.\end{array}$ (14)

Using Equations (9)-(14), model (1) is linearized about equilibrium E and takes the form

${Y}^{\prime}\left(t\right)=AY\left(t\right),$ (15)

where $Y\left(t\right)={\left[{\stackrel{\u02dc}{D}}_{S}\left(t\right),{\stackrel{\u02dc}{D}}_{I}\left(t\right),{\stackrel{\u02dc}{M}}_{S}\left(t\right),{\stackrel{\u02dc}{M}}_{I}\left(t\right)\right]}^{\text{T}}$ and A is the Jacobian matrix evaluated at E. However, the specific form of matrix A cannot be extracted due to the presence of the integral terms in (13) and (14). To bypass this issue, we use the Fundamental Theorem of linear systems of differential equations [16] and look for exponential solutions of the form

$\left[\begin{array}{c}{\stackrel{\u02dc}{D}}_{S}\left(t\right)\\ {\stackrel{\u02dc}{D}}_{I}\left(t\right)\\ {\stackrel{\u02dc}{M}}_{S}\left(t\right)\\ {\stackrel{\u02dc}{M}}_{I}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{r}_{1}\\ {r}_{2}\\ {r}_{3}\\ {r}_{4}\end{array}\right]{\text{e}}^{\lambda t}=R{e}^{\lambda t}.$ (16)

We also let $\stackrel{\u02dc}{g}$ be the (one-sided) Laplace transform of the travel-time distribution $g\left(z\right)$ . That is,

$\stackrel{\u02dc}{g}\left(x\right)\equiv {\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-xz}\text{d}z\mathrm{.}$ (17)

We have the following Lemma.

Lemma 1 The Laplace transform $\stackrel{\u02dc}{g}$ is a positive, decreasing function that is bounded above by 1 for all non-negative values of $x$ .

Proof. Let $g\left(z\right)$ be a probability density function as described above. Because the function ${\text{e}}^{-xz}$ is positive for all real x and fixed z, ${\text{e}}^{-xz}=1$ when $x=0$ , and ${\text{e}}^{-xz}$ decreases for all $x>0$ . Therefore, it must be the case that $0<g\left(z\right){\text{e}}^{-xz}\le 1$ and $g\left(z\right){\text{e}}^{-xz}$ decreases for all non-negative x. Thus, $\stackrel{\u02dc}{g}\left(x\right)\equiv {\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-xz}\text{d}z$ is a positive decreasing function bounded above by 1.

By substituting (16) into (15) and simplifying the terms, we get the specific form of matrix A, and 15 is rewritten as

$\left[\begin{array}{cccc}\frac{\partial {f}^{\left[1\right]}\left(E\right)}{\partial {D}_{S}}+{d}_{S}\stackrel{\u02dc}{g}\left(\lambda +{\delta}_{S}\right)-\lambda & \frac{\partial {f}^{\left[1\right]}\left(E\right)}{\partial {D}_{I}}& \frac{\partial {f}^{\left[1\right]}\left(E\right)}{\partial {M}_{S}}& \frac{\partial {f}^{\left[1\right]}\left(E\right)}{\partial {M}_{I}}\\ \frac{\partial {f}^{\left[2\right]}\left(E\right)}{\partial {D}_{S}}& \frac{\partial {f}^{\left[2\right]}\left(E\right)}{\partial {D}_{I}}+{d}_{I}\stackrel{\u02dc}{g}\left(\lambda +{\delta}_{I}\right)-\lambda & \frac{\partial {f}^{\left[2\right]}\left(E\right)}{\partial {M}_{S}}& \frac{\partial {f}^{\left[2\right]}\left(E\right)}{\partial {M}_{I}}\\ \frac{\partial {f}^{\left[3\right]}\left(E\right)}{\partial {D}_{S}}& \frac{\partial {f}^{\left[3\right]}\left(E\right)}{\partial {D}_{I}}& \frac{\partial {f}^{\left[3\right]}\left(E\right)}{\partial {M}_{S}}-\lambda & \frac{\partial {f}^{\left[3\right]}\left(E\right)}{\partial {M}_{I}}\\ \frac{\partial {f}^{\left[4\right]}\left(E\right)}{\partial {D}_{S}}& \frac{\partial {f}^{\left[4\right]}\left(E\right)}{\partial {D}_{I}}& \frac{\partial {f}^{\left[4\right]}\left(E\right)}{\partial {M}_{S}}& \frac{\partial {f}^{\left[4\right]}\left(E\right)}{\partial {M}_{I}}-\lambda \end{array}\right]\left[\begin{array}{c}{r}_{1}\\ {r}_{2}\\ {r}_{3}\\ {r}_{4}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right].$

(18)
The linear system in (18) has a nontrivial solution if and only if the determinant of the matrix is zero. This leads to the characteristic equation corresponding to model (1) linearized about E. Before deriving the characteristic equation, we prove the existence of DFE.

Proposition 1 The disease free equilibrium of model (1) exists if and only if ${\lambda}_{D}>{\mu}_{D}+\rho +{d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)$ and ${\lambda}_{M}>{\mu}_{M}+\sigma $ are satisfied.

Proof. Noting that ${D}_{I}^{*}=0$ , ${D}_{N}={D}_{S}^{*}$ , and $\frac{\text{d}{D}_{S}}{\text{d}t}=0$ at the DFE, the first equation in model (1) gives us ${D}_{S}^{*}=\frac{{\lambda}_{D}-\left({\mu}_{D}+\rho \right)-{d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)}{{\mu}_{2D}}$ . Similarly,

${M}_{I}^{*}=0$ and ${M}_{N}={M}_{S}^{*}$ , and the third equation of model (1) gives rise to

${M}_{S}^{*}=\frac{{\lambda}_{M}-\left({\mu}_{M}+\sigma \right)}{{\mu}_{2M}}$ . As ${D}_{S}^{*}>0$ and ${M}_{S}^{*}>0$ by parameter assumptions, the

disease free equilibrium exists.

Remark 1 The inequalities (2) and (4) and Lemma 1 imply that the conditions of Proposition 1 are always satisfied. Hence, the DFE always exists and it is given by

$\begin{array}{l}{D}_{S}^{*}=\frac{{\lambda}_{D}-\left({\mu}_{D}+\rho \right)-{d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)}{{\mu}_{2D}}\\ {D}_{I}^{*}=0\\ {M}_{S}^{*}=\frac{{\lambda}_{M}-\left({\mu}_{M}+\sigma \right)}{{\mu}_{2M}}\\ {M}_{I}^{*}=0\end{array}$ (19)

By linearizing model (1) about the DFE, we get the characteristic equation

$det\left(J\left(\lambda \right)\right)=\mathrm{0,}$ (20)

where $J\left(\lambda \right)$ is the matrix in (18) evaluated at $E=DFE$ , and it simplifies to

$J\left(\lambda \right)=\left[\begin{array}{cccc}{J}_{1}\left(\lambda \right)& {\lambda}_{D}-{\mu}_{2D}{D}_{S}^{*}& 0& -{\beta}_{D}\\ 0& {J}_{2}\left(\lambda \right)& 0& {\beta}_{D}\\ 0& -\frac{{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}& {J}_{3}\left(\lambda \right)& {\lambda}_{M}-{\mu}_{2M}{M}_{S}^{*}\\ 0& \frac{{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}& 0& {J}_{4}\left(\lambda \right)\end{array}\right]$ (21)

such that

${J}_{1}\left(\lambda \right)={\lambda}_{D}-{\mu}_{D}-\rho -{d}_{S}-2{\mu}_{2D}{D}_{S}^{*}+{d}_{S}\stackrel{\u02dc}{g}\left(\lambda +{\delta}_{S}\right)-\lambda ,$ (22)

${J}_{2}\left(\lambda \right)=-{\mu}_{D}-\rho -{\gamma}_{D}-{d}_{I}-{\mu}_{2D}{D}_{S}^{*}+{d}_{I}\stackrel{\u02dc}{g}\left(\lambda +{\delta}_{I}\right)-\lambda ,$ (23)

${J}_{3}\left(\lambda \right)={\lambda}_{M}-{\mu}_{M}-\sigma -2{\mu}_{2M}{M}_{S}^{*}-\lambda ,$ (24)

and

${J}_{4}\left(\lambda \right)=-{\mu}_{M}-\sigma -{\mu}_{2M}{M}_{S}^{*}-\lambda .$ (25)

Hence, the characteristic equation (20) is rewritten

${J}_{1}\left(\lambda \right){J}_{3}\left(\lambda \right)\left[{J}_{2}\left(\lambda \right){J}_{4}\left(\lambda \right)-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}\right]=0.$ (26)

Since ${J}_{1}\left(\lambda \right)$ and ${J}_{2}\left(\lambda \right)$ are not polynomials, the Routh-Hurwitz criteria [17] is not applicable for determining stability. However, with a specific form of $g\left(z\right)$ , we may compute the roots of the characteristic equation and determine the necessary and sufficient conditions for the stability of the DFE.

3.2. Basic Reproduction Number

The basic reproduction number ${R}_{0}$ is defined as the expected number of secondary infections produced by a single case of an infection introduced to a completely susceptible population [18] . When ${R}_{0}>1$ , the infection will spread as the number of infected individuals increases. When ${R}_{0}<1$ , the infection will die out in the long run. Thus, we seek conditions and parameter values so that ${R}_{0}<1$ .

The magnitude of ${R}_{0}$ determines the severity of infection. Larger values of ${R}_{0}>1$ lead to faster disease spread, whereas smaller values of ${R}_{0}<1$ lead to the disease dying out more rapidly. Using the Next Generation Matrix (NGM) approach [19] [20] , the expression for ${R}_{0}$ can be derived. Specifically, the next generation matrix is given by $K=F{V}^{-1}$ , and the spectral radius of K is equal to ${R}_{0}$ . The elements of matrix F, using the extended definition of the matrix F [21] , represent new infections, where the entry $\left(i\mathrm{,}j\right)$ of F represents the rate at which secondary individuals appear in class i per individual of type j. The elements of matrix V are the transition of infections.

In order to calculate the ${R}_{0}$ expression, we make some simplifying assumptions in our model. In particular, we assume the integral terms in the first and second equations of model (1) are simplified to

${\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z=\stackrel{\u02dc}{g}\left({\delta}_{S}\right){D}_{S}\left(t\right)$ (27)

and

${\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{I}z}{D}_{I}\left(t-z\right)\text{d}z=\stackrel{\u02dc}{g}\left({\delta}_{I}\right){D}_{I}\left(t\right)$ (28)

respectively.

Remark 2 The assumptions in (27) and (28) result in a positive outflow of deer out of the patch. The first equation of model (1) contains the expression

$-{d}_{S}{D}_{S}\left(t\right)+{d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z$ . Using (27), this simplifies to

${d}_{S}\left(\stackrel{\u02dc}{g}\left({\delta}_{S}\right)-1\right){D}_{S}\left(t\right)$ which is negative by the above Lemma. In other words,

there are more susceptible deer leaving the patch than entering it. The same is true for the infected deer as concluded from the second equation of model (1) and assumption (28).

Using the assumptions in (27) and (28), we get that

$F=\left[\begin{array}{cc}{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)& {\beta}_{D}\\ \frac{{\beta}_{M}{M}_{S}^{\mathrm{*}}}{{D}_{S}^{\mathrm{*}}}& 0\end{array}\right]\mathrm{,}$ (29)

$V=\left[\begin{array}{cc}{V}_{1}& 0\\ 0& {V}_{2}\end{array}\right],$ (30)

and

$F{V}^{-1}=\left[\begin{array}{cc}\frac{{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)}{{V}_{1}}& \frac{{\beta}_{D}}{{V}_{2}}\\ \frac{{\beta}_{M}{M}_{S}^{*}}{{V}_{1}{D}_{S}^{*}}& 0\end{array}\right],$ (31)

where

${V}_{1}={\mu}_{D}+\rho +{\gamma}_{D}+{d}_{I}+{\mu}_{2D}{D}_{S}^{*}$ (32)

and

${V}_{2}={\mu}_{M}+\sigma +{\mu}_{2M}{M}_{S}^{*}.$ (33)

As mentioned earlier, the basic reproduction number ${R}_{0}$ is the spectral radius of $F{V}^{-1}$ . Since $F{V}^{-1}$ is a positive definite matrix, ${R}_{0}$ is equal to the largest eigenvalue of $F{V}^{-1}$ . After simplifying, the expression for ${R}_{0}$ can be written as

${R}_{0}=\frac{1}{2}\left({R}_{0}^{\left[1\right]}+\sqrt{{\left({R}_{0}^{\left[1\right]}\right)}^{2}+4{R}_{0}^{\left[2\right]}}\right)\mathrm{,}$ (34)

where

${R}_{0}^{\left[1\right]}=\frac{{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)}{{V}_{1}}\mathrm{,}$ (35)

representing the contribution of deer migration to disease outbreaks, and

${R}_{0}^{\left[2\right]}=\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{\mathrm{*}}}{{V}_{1}{V}_{2}{D}_{S}^{\mathrm{*}}}\mathrm{,}$ (36)

representing the effects of the deer-midge interactions on disease outbreaks. Therefore, the migration effects of infected deer and the effects of deer-midge interactions within the patch on HD outbreaks can be studied separately.

1) Pure migration effects ( ${R}_{0}^{\left[2\right]}=0$ ). This occurs when either ${\beta}_{D}$ or ${\beta}_{M}$ is zero, and thus there is no transmission of the disease between the midges and the deer (or vice-versa) within the patch. Using Equation (34), ${R}_{0}^{\left[2\right]}=0$ implies ${R}_{0}={R}_{0}^{\left[1\right]}$ . In reality, this can effectively occur when the midge population in the patch is negligible. It can be seen that ${R}_{0}^{\left[1\right]}$ is a concave down increasing function of ${d}_{I}$ . Thus, the flux rate of infected deer ${d}_{I}$ may increase ${R}_{0}^{\left[1\right]}$ . From Equation (32), we get that ${\mathrm{lim}}_{{d}_{I}\to \infty}{R}_{0}^{\left[1\right]}=\stackrel{\u02dc}{g}\left({\delta}_{I}\right)$ . Using Lemma 1, $\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\le 1$ . Therefore, ${d}_{I}$ alone cannot cause an outbreak even though it increases the ${R}_{0}^{\left[1\right]}$ value. In fact, using Equations (32) and (35), it can be easily shown that ${R}_{0}^{\left[1\right]}<1$ for all parameter values of the model. Hence, assumptions (27) and (28) are underestimating the migration effects of deer population on disease outbreak.

2) Residential effects ( ${R}_{0}^{\left[1\right]}=0$ ). This occurs when ${d}_{I}=0$ , which means that infected deer have limited mobility and cannot leave or enter the patch due to illness. In this case, ${R}_{0}^{\left[1\right]}=0$ implies ${R}_{0}=\sqrt{{R}_{0}^{\left[2\right]}}$ . In this case, an epidemic may be prevented if ${R}_{0}^{\left[2\right]}<1$ . This, in fact, may be possible as the harvest rate, $\rho $ , is a part of the expression of ${R}_{0}^{\left[2\right]}$ . On the other hand, small values of ${V}_{2}$ (i.e., low mortality of midges) may result in an outbreak.

The following proposition indicates the effects of parameter values on ${R}_{0}$ in general.

Proposition 2 The basic reproduction number ${R}_{0}$ is defined in Equation (34) and it has the following properties:

1) ${R}_{0}$ is an increasing function of ${\delta}_{S}$ and ${d}_{S}$ .

2) ${R}_{0}$ is a decreasing function of ${\delta}_{I}$ .

3) ${R}_{0}$ is an increasing function of ${d}_{I}$ if ${d}_{S}$ or the product ${\beta}_{D}{\beta}_{M}$ is sufficiently small.

4) ${R}_{0}$ is a decreasing function of ${d}_{I}$ if ${d}_{S}$ or the product ${\beta}_{D}{\beta}_{M}$ is sufficiently large.

Proof. Part (i): As shown below, the partial derivative of ${R}_{0}$ with respect to ${\delta}_{S}$ is positive.

$\frac{\partial {R}_{0}}{\partial {\delta}_{S}}=\frac{-{\beta}_{D}{\beta}_{M}{d}_{S}{M}_{S}^{*}{\stackrel{\u02dc}{g}}^{\prime}\left({\delta}_{S}\right)}{{\mu}_{2D}{V}_{1}{V}_{2}\left({D}_{S}^{*}\right)\sqrt{{\left({R}_{0}^{\left[1\right]}\right)}^{2}+4{R}_{0}^{\left[2\right]}}}>0.$ (37)

Note that ${\stackrel{\u02dc}{g}}^{\prime}\left({\delta}_{S}\right)<0$ because $\stackrel{\u02dc}{g}\left({\delta}_{S}\right)$ is a decreasing function (See Lemma 1). Similarly, the partial derivative of ${R}_{0}$ with respect to ${d}_{S}$ is positive.

$\begin{array}{l}\frac{\partial {R}_{0}}{\partial {d}_{S}}=\frac{{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)}{2{V}_{1}^{2}}+\frac{1}{\sqrt{{\left({R}_{0}^{\left[1\right]}\right)}^{2}+4{R}_{0}^{\left[2\right]}}}[\frac{{\left({d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\right)}^{2}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)}{2{V}_{1}^{3}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{\mu}_{2D}{V}_{2}{\left({V}_{1}{D}_{S}^{*}\right)}^{2}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)}\left({\mu}_{2D}{D}_{S}^{*}+{V}_{1}\right)]>0.\end{array}$ (38)

Part (ii): The partial derivative of ${R}_{0}$ with respect to ${\delta}_{I}$ is negative.

$\frac{\partial {R}_{0}}{\partial {\delta}_{I}}=\frac{{d}_{I}{\stackrel{\u02dc}{g}}^{\prime}\left({\delta}_{I}\right)}{2}\left(\frac{1}{{V}_{1}}+\frac{{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)}{{V}_{1}^{2}\sqrt{{\left({R}_{0}^{\left[1\right]}\right)}^{2}+4{R}_{0}^{\left[2\right]}}}\right)<0.$ (39)

To prove statements (iii) and (iv), note that the partial derivative of ${R}_{0}$ with respect to ${d}_{I}$ is given by

$\frac{\partial {R}_{0}}{\partial {d}_{I}}=\frac{\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\left({V}_{1}-{d}_{I}\right)}{2{V}_{1}^{2}}+\frac{1}{{V}_{1}^{2}\sqrt{{\left({R}_{0}^{\left[1\right]}\right)}^{2}+4{R}_{0}^{\left[2\right]}}}\left[\frac{{d}_{I}\left({V}_{1}-{d}_{I}\right){\left(\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\right)}^{2}}{2{V}_{1}}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{V}_{2}{D}_{S}^{*}}\right]$ (40)

Also note that ${V}_{1}-{d}_{I}={\mu}_{D}+\rho +{\gamma}_{D}+{\mu}_{2D}{D}_{S}^{*}>0$ . The expression

$\frac{{d}_{I}\left({V}_{1}-{d}_{I}\right){\left(\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\right)}^{2}}{2{V}_{1}}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{V}_{2}{D}_{S}^{*}}>0$ (41)

is equivalent to

${d}_{I}{V}_{2}{D}_{S}^{*}\left({V}_{1}-{d}_{I}\right){\left(\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\right)}^{2}-2{V}_{1}{\beta}_{D}{\beta}_{M}{M}_{S}^{*}>0.$ (42)

Recall that ${D}_{S}^{*}=\frac{{\lambda}_{D}-\left({\mu}_{D}+\rho \right)-{d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)}{{\mu}_{2D}}$ . When ${d}_{S}$ is sufficiently small, ${d}_{I}{V}_{2}{D}_{S}^{\mathrm{*}}\left({V}_{1}-{d}_{I}\right){\left(\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\right)}^{2}$ will be sufficiently large and the inequality

holds. When ${\beta}_{D}{\beta}_{M}$ is sufficiently small, $2{V}_{1}{\beta}_{D}{\beta}_{M}{M}_{S}^{\mathrm{*}}$ will be sufficiently

small and the inequality holds. Thus $\frac{\partial {R}_{0}}{\partial {d}_{I}}>0$ . Similarly, when either ${d}_{S}$ or the product ${\beta}_{D}{\beta}_{M}$ is sufficiently large, $\frac{\partial {R}_{0}}{\partial {d}_{I}}<0$ .

Remark 3 Proposition 2 implies that the flux rate ${d}_{I}$ of infected deer can have two opposing effects based on the value of ${d}_{S}$ or the product ${\beta}_{D}{\beta}_{M}$ . Because the directional behavior of ${R}_{0}$ changes due to the value of these, there must be critical values ( ${d}_{S}^{\left[c\right]}$ and ${\left({\beta}_{D}{\beta}_{M}\right)}^{\left[c\right]}$ ) such that ${R}_{0}$ is an increasing function of ${d}_{I}$ when ${d}_{S}$ or ${\beta}_{D}{\beta}_{M}$ are below either of the critical values and ${R}_{0}$ is a decreasing function of ${d}_{I}$ when ${d}_{S}$ or ${\beta}_{D}{\beta}_{M}$ are above either of them.

The following Lemma is associated with the structure of the ${R}_{0}$ expression in equation (34).

Lemma 2 For $a\mathrm{,}b\ge 0$ , $a+b<1$ if and only if $\frac{1}{2}\left(a+\sqrt{{a}^{2}+4b}\right)<1$ .

Proof. (Þ) If $a+b<1$ , then $b<1-a$ . Also, as $0\le a<1$ , $\left|a-2\right|=2-a$ . Thus,

$\begin{array}{c}\frac{1}{2}\left(a+\sqrt{{a}^{2}+4b}\right)<\frac{1}{2}\left(a+\sqrt{{a}^{2}+4\left(1-a\right)}\right)=\frac{1}{2}\left(a+\sqrt{{a}^{2}-4a+4}\right)\\ =\frac{1}{2}\left(a+\left|a-2\right|\right)=\frac{1}{2}\left(a+2-a\right)=1\end{array}$ (43))

(Ü)

$\begin{array}{l}\frac{1}{2}\left(a+\sqrt{{a}^{2}+4b}\right)<1\\ a+\sqrt{{a}^{2}+4b}<2\\ \sqrt{{a}^{2}+4b}<2-a\end{array}$ (44)

$\begin{array}{l}{a}^{2}+4b<4-4a+{a}^{2}\\ 4a+4b<4\\ a+b<1\end{array}$

Remark 4 Let $a={R}_{0}^{\left[1\right]}$ and $b={R}_{0}^{\left[2\right]}$ . Using Lemma 2, we get that ${R}_{0}<1$ is equivalent to ${R}_{0}^{\left[1\right]}+{R}_{0}^{\left[2\right]}<1$ . As indicated in [22] [23] , the expression ${R}_{0}^{\left[1\right]}+{R}_{0}^{\left[2\right]}$ is known as a Type-Reduction number which can be more accurate than ${R}_{0}$ to calculate the minimum disease eradication efforts.

Proposition 3 Under the assumptions (27) and (28), the DFE of model (1) is locally asymptotically stable if and only if ${R}_{0}^{\left[1\right]}+{R}_{0}^{\left[2\right]}<1$ or, equivalently, ${R}_{0}<1$ .

Proof. (Ü) We determine stability conditions at the DFE by using the Jacobian of the system of equations. The DFE is locally asymptotically stable if the real parts of all eigenvalues of the Jacobian matrix are negative as explained in Section 3.1. Using assumptions (27) and (28), the Jacobian matrix evaluated at the DFE is given by:

$A=\left[\begin{array}{cccc}{A}_{1}& {\lambda}_{D}-{\mu}_{2D}{D}_{S}^{*}& 0& -{\beta}_{D}\\ 0& {A}_{2}& 0& {\beta}_{D}\\ 0& -\frac{{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}& {A}_{3}& {\lambda}_{M}-{\mu}_{2M}{M}_{S}^{*}\\ 0& \frac{{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}& 0& {A}_{4}\end{array}\right],$ (45)

where ${A}_{1}={\lambda}_{D}-{\mu}_{D}-\rho -\left({d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)+2{\mu}_{2D}{D}_{S}^{*}\right)$ , ${A}_{2}={d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)-{V}_{1}$ , ${A}_{3}={\lambda}_{M}-{\mu}_{M}-\sigma -s{\mu}_{2M}{M}_{S}^{*}$ , and ${A}_{4}=-{V}_{2}$ . The characteristic equation of this matrix, using $\Lambda $ for the eigenvalues, is

$f\left(\Lambda \right)=\left({A}_{1}-\Lambda \right)\left({A}_{3}-\Lambda \right)\left[\left({A}_{2}-\Lambda \right)\left({A}_{4}-\Lambda \right)-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}\right].$ (46)

For the first eigenvalue ${A}_{1}$ , we note that since the DFE must satisfy ${{D}^{\prime}}_{S}=0$ , we can show that ${\lambda}_{D}={\mu}_{D}+\rho +{d}_{S}+{\mu}_{2D}{D}_{S}^{*}-{d}_{S}\stackrel{\u02dc}{g}\left({\delta}_{S}\right)={V}_{1}-{d}_{S}\stackrel{\u02dc}{g}\left({\delta}_{S}\right)$ . Therefore,

$\begin{array}{l}{A}_{1}{\lambda}_{D}-{\mu}_{D}-\rho -\left({d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)+2{\mu}_{2D}{D}_{S}^{*}\right)\\ =\left({\mu}_{D}+\rho +{d}_{S}+{\mu}_{2D}{D}_{S}^{*}-{d}_{S}\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)-{\mu}_{D}-\rho -\left({d}_{S}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{S}\right)\right)+2{\mu}_{2D}{D}_{S}^{*}\right)\\ =-{\mu}_{2D}{D}_{S}^{*}<0\end{array}$ (47)

Similarly, for the second eigenvalue, given that the DFE must satisfy ${{M}^{\prime}}_{S}=0$ , we can show ${\lambda}_{M}={\mu}_{M}+\sigma +{\mu}_{2M}{M}_{S}^{*}$ , and thus ${A}_{3}={\lambda}_{M}-{\mu}_{M}-\sigma -2{\mu}_{2M}{M}_{S}^{*}=-{\mu}_{2M}{M}_{S}^{*}<0$ .

For the remaining two eigenvalues, we rewrite the part of the characteristic equation in brackets as

${\Lambda}^{2}-\left({A}_{2}+{A}_{4}\right)\Lambda +{A}_{2}{A}_{4}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}=0.$ (48)

This is a quadratic of the form ${\Lambda}^{2}+b\Lambda +c$ . According to the Routh-Hurwitz criteria [17] , the roots of a quadratic will have negative real parts if the linear coefficient and the constant term are positive. The linear coefficient is $-\left({A}_{2}+{A}_{4}\right)$ and is positive as shown below.

$\begin{array}{c}{A}_{2}+{A}_{4}={d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)-{V}_{1}-{V}_{2}\\ ={d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)-{\mu}_{D}-\rho -{\gamma}_{D}-{d}_{I}-{\mu}_{2D}{D}_{S}^{*}-{V}_{2}\\ =-{d}_{I}\left(1-\stackrel{\u02dc}{g}\left({\delta}_{I}\right)\right)-{\mu}_{D}-\rho -{\gamma}_{D}-{d}_{I}-{\mu}_{2D}{D}_{S}^{*}-{V}_{2}<0\end{array}$ (49)

If ${R}_{0}^{\left[1\right]}+{R}_{0}^{\left[2\right]}<1$ , then

$\begin{array}{l}\frac{{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)}{{V}_{1}}+\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{V}_{1}{V}_{2}{D}_{S}^{*}}<1\\ {d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right){V}_{2}+\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}<{V}_{1}{V}_{2}\\ {V}_{1}{V}_{2}-{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right){V}_{2}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}>0\end{array}$ (50)

Hence, the constant term of the characteristic equation

$\begin{array}{c}{A}_{2}{A}_{4}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{\mathrm{*}}}{{D}_{S}^{\mathrm{*}}}=-\left({d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right)-{V}_{1}\right){V}_{2}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{\mathrm{*}}}{{D}_{S}^{\mathrm{*}}}\\ ={V}_{1}{V}_{2}-{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right){V}_{2}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{\mathrm{*}}}{{D}_{S}^{\mathrm{*}}}>0\end{array}$ (51)

Therefore, both roots of the quadratic (i.e. the two eigenvalues) must have negative real parts. Thus, under the given conditions, the system is stable at DFE.

(Þ) If the DFE of model (1) is locally asymptotically stable, then by Theorem 8.12. iii of [24], the real parts of all eigenvalues of the Jacobian matrix A are

negative. By (50), this occurs when ${V}_{1}{V}_{2}-{d}_{I}\stackrel{\u02dc}{g}\left({\delta}_{I}\right){V}_{2}-\frac{{\beta}_{D}{\beta}_{M}{M}_{S}^{*}}{{D}_{S}^{*}}>0$ which is the same as ${R}_{0}^{\left[1\right]}+{R}_{0}^{\left[2\right]}<1$ .

We must now prove the existence of an endemic equilibrium solution in the proposed model. However, this is difficult as two of the variables, ${D}_{S}$ and ${D}_{I}$ , are contained within the integral dispersion terms. Therefore, we utilize a technique called the chain trick [15] to reduce model (1) to an ODE model.

3.3. Reduction to ODE Model

Using the chain trick method [15] , we can rewrite the first two equations as

$\begin{array}{l}\frac{\text{d}{D}_{S}}{\text{d}t}={\lambda}_{D}{D}_{N}-\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left(\left({\mu}_{D}+\rho \right)+{d}_{S}+{\mu}_{2D}{D}_{N}\right){D}_{S}+{\stackrel{\xaf}{D}}_{S}\\ \frac{\text{d}{D}_{I}}{\text{d}t}=\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left(\left({\mu}_{D}+\rho \right)+{\gamma}_{D}+{d}_{I}+{\mu}_{2D}{D}_{N}\right){D}_{I}+{\stackrel{\xaf}{D}}_{I}\end{array}$ (52)

where

${\stackrel{\xaf}{D}}_{S}={d}_{S}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z$ (53)

and

${\stackrel{\xaf}{D}}_{I}={d}_{I}{\displaystyle {\int}_{0}^{\infty}}g\left(z\right){\text{e}}^{-{\delta}_{I}z}{D}_{I}\left(t-z\right)\text{d}z.$ (54)

These quantities are treated as new model variables, so we may now differentiate both of them and amend the existing set of equations.

In time delay models, there are two distributions that are commonly used. The first is a uniform distribution with mean $\tau $ given by

$g\left(u\right)=(\begin{array}{l}\frac{1}{\tau \rho}\mathrm{,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\tau \left(1-\frac{\rho}{2}\right)\le u\le \tau \left(1+\frac{\rho}{2}\right)\\ \mathrm{0,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{elsewhere}\text{.}\end{array}$ (55)

The second is the gamma distribution given by

$g\left(u\right)=\frac{{u}^{p-1}{\alpha}^{p}{\text{e}}^{-\alpha u}}{\Gamma \left(p\right)}\mathrm{,}$ (56)

where $\alpha \mathrm{,}p\ge 0$ are parameters which determine the shape of the distribution and the mean of the distribution is $p/\alpha $ . In the case when $p=1$ , the result is the exponential distribution, $g\left(z\right)=\alpha {\text{e}}^{-\alpha z}$ . Using (56) with $p=1$ , the

expression for $\frac{\text{d}{\stackrel{\xaf}{D}}_{S}}{\text{d}t}$ is computed to be:

$\begin{array}{c}{\stackrel{\xaf}{D}}_{S}={\displaystyle {\int}_{0}^{\infty}}\text{\hspace{0.05em}}g\left(z\right){\text{e}}^{-{\delta}_{S}z}{D}_{S}\left(t-z\right)\text{d}z\\ ={\displaystyle {\int}_{-\infty}^{t}}\text{\hspace{0.05em}}g\left(t-u\right){\text{e}}^{-{\delta}_{S}\left(t-u\right)}{D}_{S}\left(u\right)\text{d}u\\ ={\text{e}}^{-{\delta}_{S}t}{\displaystyle {\int}_{-\infty}^{t}}\text{\hspace{0.05em}}g\left(t-u\right){\text{e}}^{-{\delta}_{S}u}{D}_{S}\left(u\right)\text{d}u\\ ={\text{e}}^{-{\delta}_{S}t}{\displaystyle {\int}_{-\infty}^{t}}\text{\hspace{0.05em}}\alpha {\text{e}}^{-\alpha \left(t-u\right)}{\text{e}}^{-{\delta}_{S}u}{D}_{S}\left(u\right)\text{d}u\\ =\alpha {\text{e}}^{-\left({\delta}_{S}+\alpha \right)t}{\displaystyle {\int}_{-\infty}^{t}}\text{\hspace{0.05em}}{\text{e}}^{\left({\delta}_{S}+\alpha \right)u}{D}_{S}\left(u\right)\text{d}u\end{array}$ (57)

Thus, by the product rule for differentiation,

$\begin{array}{c}\frac{\text{d}{\stackrel{\xaf}{D}}_{S}}{\text{d}t}=\alpha \left(-\left({\delta}_{S}+\alpha \right)\right){\text{e}}^{-\left({\delta}_{S}+\alpha \right)t}{\displaystyle {\int}_{-\infty}^{t}}\text{\hspace{0.05em}}{\text{e}}^{\left({\delta}_{S}+\alpha \right)u}{D}_{S}\left(u\right)\text{d}u+{d}_{S}\alpha {\text{e}}^{-\left({\delta}_{S}+\alpha \right)t}{\text{e}}^{\left({\delta}_{S}+\alpha \right)u}{D}_{S}\left(t\right)\\ =-\left({\delta}_{S}+\alpha \right){\stackrel{\xaf}{D}}_{S}+\alpha {D}_{S}\end{array}$ (58)

The simplification is the same for ${\stackrel{\xaf}{D}}_{I}$ , and so the delayed model in (1) is reduced to the ODE model formulated by

$\begin{array}{l}\frac{\text{d}{D}_{S}}{\text{d}t}={\lambda}_{D}{D}_{N}-\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left({\mu}_{D}+\rho +{d}_{S}+{\mu}_{2D}{D}_{N}\right){D}_{S}+{d}_{S}{\stackrel{\xaf}{D}}_{S}\\ \frac{\text{d}{D}_{I}}{\text{d}t}=\frac{{\beta}_{D}{M}_{I}{D}_{S}}{{D}_{N}}-\left({\mu}_{D}+\rho +{\gamma}_{D}+{d}_{I}+{\mu}_{2D}{D}_{N}\right){D}_{I}+{d}_{I}{\stackrel{\xaf}{D}}_{I}\\ \frac{\text{d}{M}_{S}}{\text{d}t}={\lambda}_{M}{M}_{N}-\frac{{\beta}_{M}{D}_{I}{M}_{S}}{{D}_{N}}-\left({\mu}_{M}+\sigma +{\mu}_{2M}{M}_{N}\right){M}_{S}\\ \frac{\text{d}{M}_{I}}{\text{d}t}=\frac{{\beta}_{M}{D}_{I}{M}_{S}}{{D}_{N}}-\left({\mu}_{M}+\sigma +{\mu}_{2M}{M}_{N}\right){M}_{I}\\ \frac{\text{d}{\stackrel{\xaf}{D}}_{S}}{\text{d}t}=-\left({\delta}_{S}+\alpha \right){\stackrel{\xaf}{D}}_{S}+\alpha {D}_{S},\text{\hspace{0.17em}}\frac{\text{d}{\stackrel{\xaf}{D}}_{I}}{\text{d}t}=-\left({\delta}_{I}+\alpha \right){\stackrel{\xaf}{D}}_{I}+\alpha {D}_{I}\end{array}$ (59)

The disease free equilibrium (DFE) is computed to be

$\begin{array}{l}{D}_{S}^{*}\frac{\alpha +\left({\delta}_{S}+\alpha \right)\left({\lambda}_{D}-\left({\mu}_{D}+\rho \right)-{d}_{S}\right)}{\left({\delta}_{S}+\alpha \right){\mu}_{2D}}\\ {D}_{I}^{*}=0\\ {M}_{S}^{*}=\frac{{\lambda}_{M}-\left({\mu}_{M}+\sigma \right)}{{\mu}_{2M}}\\ {M}_{I}^{*}=0\\ {\stackrel{\xaf}{D}}_{S}^{*}=\frac{\alpha \left[{d}_{S}\alpha +\left({\delta}_{D}+\alpha \right)\left({\lambda}_{D}-\left({\mu}_{D}+\rho \right)-{d}_{S}\right)\right]}{{\left({\delta}_{S}+\alpha \right)}^{2}{\mu}_{2D}}\\ {\stackrel{\xaf}{D}}_{I}^{*}=0\end{array}$ (60)

In the next section, we provide the numerical simulations of the ODE model (59) and the ${R}_{0}$ expression (34).

4. Numerical Simulations

Using Matlab 9.1, we generated the surface plots of ${R}_{0}$ values based on the model parameters ${R}_{0}^{\left[1\right]}$ and ${R}_{0}^{\left[2\right]}$ (See Figure 2). As proven in Proposition 2,

Figure 2. Numerical simulations of ${R}_{0}$ as a function of the selected model parameters. (a) ${R}_{0}$ values increase with ${d}_{I}$ provided ${d}_{S}$ values are small. When ${d}_{S}$ values are large, ${R}_{0}$ decreases with ${d}_{I}$ ; (b) ${R}_{0}$ increases both with ${\beta}_{D}{\beta}_{M}$ and ${d}_{I}$ ; (c) ${R}_{0}$ increases with ${\delta}_{S}$ and decreases with ${\delta}_{I}$ ; (d) ${R}_{0}$ increases linearly with ${R}_{0}^{\left[1\right]}$ and increases parabolically with ${R}_{0}^{\left[2\right]}$ .

Figure 2(a) shows that ${R}_{0}$ is an increasing function with respect to ${d}_{S}$ . The influx of additional, susceptible deer into a patch leads to an increased number of potential interactions with infected midges and thus an increase in the number of infections overall. Figure 2(c) shows that ${R}_{0}$ is an increasing function with respect to ${\delta}_{S}$ and a decreasing function with respect to ${\delta}_{I}$ . Figure 2(a) and Figure 2(b) demonstrate the behavior of ${R}_{0}$ with respect to the influx of infected deer, ${d}_{I}$ . For smaller values of ${d}_{S}$ or ${\beta}_{D}{\beta}_{M}$ , ${R}_{0}$ is an increasing function with respect to ${d}_{I}$ ; for larger values of ${d}_{S}$ or ${\beta}_{D}{\beta}_{M}$ , it is a decreasing function with respect to ${d}_{I}$ . Thus, there must be a critical value ( ${d}_{S}^{\left[c\right]}$ or ${\left({\beta}_{D}{\beta}_{M}\right)}^{\left[c\right]}$ ) where the behavior changes.

If we consider ${R}_{0}$ as a function of the deer-midge interactions, then ${R}_{0}$ is essentially a linear function of ${R}_{0}^{\left[1\right]}$ and a function of the square root of ${R}_{0}^{\left[2\right]}$ . The graph of ${R}_{0}$ would be increasing and concave down with respect to an increase in ${R}_{0}^{\left[2\right]}$ (See Figure 2(d)). This is consistent with what we would expect to happen. As the amount of interaction increases, so does the number of potential new infections with a greater chance of an outbreak occurring. Plus, as a greater proportion of the deer population becomes infected, the rate of increase of ${R}_{0}$ must decrease as the number of uninfected deer will consequently drop.

We also demonstrate numerically that the solutions of model (1) converge to the endemic equilibrium if ${R}_{0}>1$ and achieves a disease free equilibrium if ${R}_{0}<1$ . To do this, a MATLAB code was written utilizing the ODE45 solver, and the results were verified against the computed ${R}_{0}$ value for a given set of parameters. At time $t=0$ , we have the following initial values: ${D}_{S}\left(0\right)=30$ , ${D}_{I}\left(0\right)=10$ , ${M}_{S}\left(0\right)=20$ , ${M}_{I}\left(0\right)=5$ , ${\stackrel{\xaf}{D}}_{S}\left(0\right)=10$ , and ${\stackrel{\xaf}{D}}_{I}\left(0\right)=1$ . See Table 2 for the specific parameter values used for the numerical simulations.

Figure 3(a) and Figure 3(c) show the long-term behavior of the four classes of deer populations-total susceptible, total infected, susceptible influx, and infected influx-plotted on the same graph, while Figure 3(b) and Figure 3(d) show the long-term behavior of the susceptible and infected midge populations.

Table 2. Parameter values used in model simulation and the calculated ${R}_{0}$ values.

Note: The ${R}_{0}$ values are consistent with the numerical simulations shown in Figure 3. Similar results were obtained using different sets of parameter values.

Figure 3. (a) (b) When the basic reproduction number ${R}_{0}<1$ , the system stabilizes to its disease free equilibrium and the number of infected deer, the number of dispersing infected deer, and the number of infected midges tends to zero as t increases; (c) (d) When the basic reproduction number ${R}_{0}>1$ , the system stabilizes to its endemic equilibrium. See Table 2 for the specific values used and the corresponding values of ${R}_{0}$ .

Figure 3(a) and Figure 3(b) indicate that when ${R}_{0}<1$ , the system will stabilize to its disease free equilibrium. Figure 3(c) and Figure 3(d) show that when ${R}_{0}>1$ , the system will stabilize to an endemic equilibrium. These outcomes are robust for large sets of initial values and parameter values.

5. Discussion

In this paper, we have developed a distributed delay model for transmission dynamics of HD in a deer population. Though mathematical models for disease and HD specifically are established, we chose to focus on how the dynamics are affected by the dispersion (migration) of deer specifically and how the basic reproduction number is affected by these dispersion rates (i.e., ${d}_{S}$ and ${d}_{I}$ ). The results show that there are critical values for the interaction parameters ${\left({\beta}_{D}{\beta}_{M}\right)}^{\left[c\right]}$ and rates of susceptible deer dispersion ${d}_{S}^{\left[c\right]}$ . Hence, possible outbreaks could be avoided by controlling how and where these deer move.

One of the primary limitations of this study is the lack of actual parameter values. Although the qualitative behavior of model (1) remains fairly distinctive, (i.e., convergence to DFE or EE) for large sets of parameter values, many of the values were chosen randomly. It is our goal to estimate some of the parameter values using data from the Missouri Department of Conservation concerning the prevalence of HD in Missouri’s white-tailed deer. Nevertheless, the graphs presented in Figure 2 and Figure 3 show consistent tendencies in the behavior in the model. We also have not considered behavior in a multi-patch system, where migrating individuals leave one patch and eventually enter a neighboring patch, nor did we consider a delay in the traveling time. Holt [25] and Weisser et al. [26] extended their results to a system of multiple patches joined through a pool of dispersing individuals. Moreover, the proposed model (1) does not include the effect of predators on the population of white-tailed deer. As a prey species, deer are linked with local predators. In Missouri, the coyote is one such predator. Some coyote predator studies have been done, but these are admittedly outdated. However, deer make up a portion of a coyote’s diet and that large increases or decreases in predator populations may influence deer mortality rates [28] . Finally, our model assumed only one vector for the transmission of HD. With the species richness of the Culicoides genus, we may reasonably expect more and different interaction rates and different levels of control efficacy [27] . We also note that weather has an effect on both the midge population and the life cycle of the HD virus [2] [29] . Midge populations thrive in damper areas, and in 2012, there was an above average amount of rain in the late winter/early spring, filling ponds and other water bodies in Missouri [28] . In addition, record warm temperatures in that spring and summer may cause midges to become more active sooner than normal [28] . Next, the high temperatures caused water sources to dry up, and not only did the resulting mud flats become ideal breeding areas for subsequent generations of midges, but also caused deer to visit water sources more frequently due to lower water content in the plants they ate as part of their diet. These same high temperatures also cause female midges to lay more eggs, and Wittmann et al. also revealed that higher temperatures decrease the extrinsic incubation period of the HD virus within the midges [30] . Thus, the virus develops faster and allows a midge to infect more deer during its life span. None of these factors have been considered in the model (1). Instead, the main focus has been on migration effects of deer population on overall HD dynamics within a patch.

6. Conclusion

The above mentioned limitations demand model extensions to study the effectiveness of control and preventive strategies. Deer species are important members of the ecosystem as they feed on brush and grass in a given area and keep them in check. In conclusion, the present work is the first step towards inclusion of migration effects of deer population modeling of HD dynamics. The ${R}_{0}$ expression provides insights into the effects of deer movement on the spread of disease.

References

[1] Flinn, E. and Sumners, J. (2013a) Breaking Down the Hemorrhagic Disease Outbreak. Missouri Conservationist, 74, 24-29.

[2] Sleeman, J.M., Howell, J.E., Knox, W.M. and Stenger, P.J. (2009) Incidence of Hemorrhagic Disease in White-Tailed Deer Is Associated with Winter and Summer Climactic Conditions. EcoHealth, 6, 11-15.

https://doi.org/10.1007/s10393-009-0220-6

[3] Foster, N.M., Breckon, R.D., Luedke, A.J., Jones, R.H. and Metcalf, H.E. (1977) Transmission of Two Strains of Epizootic Hemorrhagic Disease Virus in Deer by Culicoides variipennis. Journal of Wildlife Diseases, 13, 9-16.

https://doi.org/10.7589/0090-3558-13.1.9

[4] Ross, R. (1911) The Prevention of Malaria. John Murray, London.

[5] Macdonald, G. (1957) The Epidemiology and Control of Malaria. Oxford University Press, London.

[6] Lou, Y. and Zhao, X.-Q. (2009) The Periodic Ross-Macdonald Model with Diffusion and Advection. Applicable Analysis, 89, 1067-1089.

https://doi.org/10.1080/00036810903437804

[7] Wang, W. and Zhao, X.-Q. (2011) A Nonlocal and Time-Delayed Reaction Diffusion Model of Dengue Transmission. SIAM Journal of Applied Mathematics, 71, 147-168. https://doi.org/10.1137/090775890

[8] Fitzgibbon, W.E., Morgan, J.J. and Webb, G.B. (2017) An Outbreak Vector-Host Epidemic Model with Spatial Structure: The 2015-2016 Zika Outbreak in Rio De Janeiro. Theoretical Biology and Medical Modeling, 14, 7.

https://doi.org/10.1186/s12976-017-0051-z

[9] Neubert, M.G., Klepac, P. and Van Den Driessche, P. (2001) Stabilizing Dispersal Delays in Predator-Prey Metapopulations Models. Theoretical Population Biology, 61, 339-347.

https://doi.org/10.1006/tpbi.2002.1578

[10] Klepac, P., Neuber, M.G. and Van Den Driessche, P. (2007) Dispersal Delays, Predator-Prey Stability, and the Paradox of Enrichment. Theoretical Population Biology, 7, 436-444.

https://doi.org/10.1016/j.tpb.2007.02.002

[11] Park, A.W., Magori, K., White, B.A. and Stallknecht, D.E. (2013) When More Transmission Equals Less Disease: Reconciling the Disconnect between Disease Hotspots and Parasite Transmission. PLoS ONE, 8, e61501.

https://doi.org/10.1371/journal.pone.0061501

[12] Purdue University, Extension E-250-W. Biting Midges: Biology and Public Health Risk.

https://extension.entm.purdue.edu/publichealth/insects/bitingmidge.html

[13] Sedda, L., Brown, H.E., Purse, B.V., Burgin, L., Gloster, J. and Rogers, D.J. (2012) A New Algorithm Quantifies the Roles of Wind and Midge Flight Activity in the Bluetongue Epizootic in Northwest Europe. Proceedings: Biological Sciences, 279, 235462.

http://www.jstor.org.proxy.library.umkc.edu/stable/41549546

https://doi.org/10.1098/rspb.2011.2555

[14] Ngwa, G. and Shu, W. (1999) A Mathematical Model for Endemic Malaria with Variable Human and Mosquito Populations. United Nations Educational Scientific and Cultural Organization and International Atomic Energy Agency, The Abdus Salam International Centre for Theoretical Physics, Miramare-Trieste.

[15] Kuang, Y. (1993) Delay Differential Equations with Applications in Population Dynamics. Academic Press, Inc., San Diego.

[16] Perko, L. (2001) Differential Equations and Dynamical Systems. 3rd Edition, Springer, New York.

https://doi.org/10.1007/978-1-4613-0003-8

[17] Routh, E.J. (1877) A Treatise on the Stability of a Given State of Motion: Particularly Steady Motion. Macmillan, London.

[18] Dietz, K. (1993) The Estimation of the Basic Reproduction Number for Infectious Diseases. Statistical Methods in Medical Research, 2, 23-41.

https://doi.org/10.1177/096228029300200103

[19] Bani-Yaghoub, M., Gautam, R., Ivanek, R., van den Driessche, P. and Shuai, Z. (2012) Reproduction Numbers for Infections with Free-Living Pathogens Growing in the Environment. Journal of Biological Dynamics, 6, 923-940.

https://doi.org/10.1080/17513758.2012.693206

[20] Van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48.

https://doi.org/10.1016/S0025-5564(02)00108-6

[21] Hurford, A., Cownden, D. and Day, T. (2010) Next-Generation Tools for Evolutionary Invasion Analyses. Journal of the Royal Society Interface, 7, 561-571.

https://doi.org/10.1098/rsif.2009.0448

[22] Heesterbeek, J.A.P. and Dietz, K. (1996) The Concept of R0 in Epidemic Theory. Statistica Neerlandica, 50, 89-110.

https://doi.org/10.1111/j.1467-9574.1996.tb01482.x

[23] Heesterbeek, J.A.P. and Roberts, M.G. (2007) The Type-Reproduction Number T in Models for Infectious Disease Control. Mathematical Biosciences, 206, 3-10.

https://doi.org/10.1016/j.mbs.2004.10.013

[24] Jordan, D.W. and Smith, P. (1999) Nonlinear Ordinary Differential Equations: An Introduction to Dynamical Systems. Oxford University Press, Oxford.

[25] Holt, R.D. (1984) Spatial Heterogeneity, Indirect Interaction, and the Coexistence of Prey Species. American Naturalist, 124, 377-406.

https://doi.org/10.1086/284280

[26] Weisser, W.W., Jansen, V.A.A. and Hassell, M.P. (1997) The Effects of a Pool of Dispersers on Host-Parasitoid Systems. Journal of Theoretical Biology, 189, 413-425.

https://doi.org/10.1006/jtbi.1997.0529

[27] Park, A.W., Cleveland, C.A., Dallas, T.A. and Corn, J.L. (2016) Vector Species Richness Increases Haemorrhagic Disease Prevalence through Functional Diversity Modulating the Duration of Seasonal Transmission. Parasitology, 143, 874-879.

https://doi.org/10.1017/S0031182015000578

[28] Flinn, E. and Sumners, J. (2013b) State of the States Deer Herd. Missouri Conservationist, 74, 24-29.

[29] Mellor, P.S., Boorman, J. and Baylis, M. (2000) Culicoides Biting Midges: Their Role as Arbovirus Vectors. Annual Review of Entomology, 45, 307-340.

https://doi.org/10.1146/annurev.ento.45.1.307

[30] Wittmann, E.J., Mellor, P.S. and Baylis, M. (2002) Effect of Temperature on the Transmission of Orbiviruses by the Biting Midge, Culicoides sonorensis. Medical and Veterinary Entomology, 16, 147-156.

https://doi.org/10.1046/j.1365-2915.2002.00357.x