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   . 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  . 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  . 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     ). The model by Nobel Prize winner Ronald Ross  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  . In fact, there have been several extensions to the Ross-Macdonald model. For instance, Lou and Zhou  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  and Zika  . 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.  , 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   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.  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 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 susceptible deer, infected deer, susceptible midges, and infected midges. The total deer population at time t is , and the total midge population is . 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  . 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  . 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
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 and migration out of the patch is denoted by and . See Table 1 for a summary of the parameters and variables.
formulated with logistic growth models. These are the terms that include , , , and in model (1). Similar to  , the carrying capacity for the deer population exists and must be positive. Hence, it is required that
Also, the carrying capacity for midges exists and is positive. Thus,
Individual deer immigrate from the patch at a constant per capita rate ( and ) 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  . The function 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 is the probability that a successfully dispersing deer departing at time t completes the trip between time and . As is a probability density function, it is normalized so that . The functions and in the integrals are the probabilities of a deer surviving a trip of duration z given constant probabilities per unit of time and 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., and ) 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
we rewrite the first equation as
Similarly, we rewrite the second equation as
As the bottom two equations of model (1) have no integral term, we let and equal the right-hand side of the third and fourth equations in model (1), respectively. Let a solution of model (1) nearby an equilibrium solution be in the form of
for some , , , and . 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.
We know that equilibrium E satisfies the first equation of model (1), hence
Substituting (12) into (10) yields
Applying the same procedure to equation (7), we get that the second equation of model (1) is linearized by
Using Equations (9)-(14), model (1) is linearized about equilibrium E and takes the form
where 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  and look for exponential solutions of the form
We also let be the (one-sided) Laplace transform of the travel-time distribution . That is,
We have the following Lemma.
Lemma 1 The Laplace transform is a positive, decreasing function that is bounded above by 1 for all non-negative values of .
Proof. Let be a probability density function as described above. Because the function is positive for all real x and fixed z, when , and decreases for all . Therefore, it must be the case that and decreases for all non-negative x. Thus, 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
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 and are satisfied.
Proof. Noting that , , and at the DFE, the first equation in model (1) gives us . Similarly,
and , and the third equation of model (1) gives rise to
. As and 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
By linearizing model (1) about the DFE, we get the characteristic equation
where is the matrix in (18) evaluated at , and it simplifies to
Hence, the characteristic equation (20) is rewritten
Since and are not polynomials, the Routh-Hurwitz criteria  is not applicable for determining stability. However, with a specific form of , 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 is defined as the expected number of secondary infections produced by a single case of an infection introduced to a completely susceptible population  . When , the infection will spread as the number of infected individuals increases. When , the infection will die out in the long run. Thus, we seek conditions and parameter values so that .
The magnitude of determines the severity of infection. Larger values of lead to faster disease spread, whereas smaller values of lead to the disease dying out more rapidly. Using the Next Generation Matrix (NGM) approach   , the expression for can be derived. Specifically, the next generation matrix is given by , and the spectral radius of K is equal to . The elements of matrix F, using the extended definition of the matrix F  , represent new infections, where the entry 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 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
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
. Using (27), this simplifies to
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
As mentioned earlier, the basic reproduction number is the spectral radius of . Since is a positive definite matrix, is equal to the largest eigenvalue of . After simplifying, the expression for can be written as
representing the contribution of deer migration to disease outbreaks, and
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 ( ). This occurs when either or 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), implies . In reality, this can effectively occur when the midge population in the patch is negligible. It can be seen that is a concave down increasing function of . Thus, the flux rate of infected deer may increase . From Equation (32), we get that . Using Lemma 1, . Therefore, alone cannot cause an outbreak even though it increases the value. In fact, using Equations (32) and (35), it can be easily shown that 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 ( ). This occurs when , which means that infected deer have limited mobility and cannot leave or enter the patch due to illness. In this case, implies . In this case, an epidemic may be prevented if . This, in fact, may be possible as the harvest rate, , is a part of the expression of . On the other hand, small values of (i.e., low mortality of midges) may result in an outbreak.
The following proposition indicates the effects of parameter values on in general.
Proposition 2 The basic reproduction number is defined in Equation (34) and it has the following properties:
1) is an increasing function of and .
2) is a decreasing function of .
3) is an increasing function of if or the product is sufficiently small.
4) is a decreasing function of if or the product is sufficiently large.
Proof. Part (i): As shown below, the partial derivative of with respect to is positive.
Note that because is a decreasing function (See Lemma 1). Similarly, the partial derivative of with respect to is positive.
Part (ii): The partial derivative of with respect to is negative.
To prove statements (iii) and (iv), note that the partial derivative of with respect to is given by
Also note that . The expression
is equivalent to
Recall that . When is sufficiently small, will be sufficiently large and the inequality
holds. When is sufficiently small, will be sufficiently
small and the inequality holds. Thus . Similarly, when either or the product is sufficiently large, .
Remark 3 Proposition 2 implies that the flux rate of infected deer can have two opposing effects based on the value of or the product . Because the directional behavior of changes due to the value of these, there must be critical values ( and ) such that is an increasing function of when or are below either of the critical values and is a decreasing function of when or are above either of them.
The following Lemma is associated with the structure of the expression in equation (34).
Lemma 2 For , if and only if .
Proof. (Þ) If , then . Also, as , . Thus,
Remark 4 Let and . Using Lemma 2, we get that is equivalent to . As indicated in   , the expression is known as a Type-Reduction number which can be more accurate than 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 or, equivalently, .
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:
where , , , and . The characteristic equation of this matrix, using for the eigenvalues, is
For the first eigenvalue , we note that since the DFE must satisfy , we can show that . Therefore,
Similarly, for the second eigenvalue, given that the DFE must satisfy , we can show , and thus .
For the remaining two eigenvalues, we rewrite the part of the characteristic equation in brackets as
This is a quadratic of the form . According to the Routh-Hurwitz criteria  , the roots of a quadratic will have negative real parts if the linear coefficient and the constant term are positive. The linear coefficient is and is positive as shown below.
If , then
Hence, the constant term of the characteristic equation
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 , the real parts of all eigenvalues of the Jacobian matrix A are
negative. By (50), this occurs when which is the same as .
We must now prove the existence of an endemic equilibrium solution in the proposed model. However, this is difficult as two of the variables, and , are contained within the integral dispersion terms. Therefore, we utilize a technique called the chain trick  to reduce model (1) to an ODE model.
3.3. Reduction to ODE Model
Using the chain trick method  , we can rewrite the first two equations as
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 given by
The second is the gamma distribution given by
where are parameters which determine the shape of the distribution and the mean of the distribution is . In the case when , the result is the exponential distribution, . Using (56) with , the
expression for is computed to be:
Thus, by the product rule for differentiation,
The simplification is the same for , and so the delayed model in (1) is reduced to the ODE model formulated by
The disease free equilibrium (DFE) is computed to be
In the next section, we provide the numerical simulations of the ODE model (59) and the expression (34).
4. Numerical Simulations
Using Matlab 9.1, we generated the surface plots of values based on the model parameters and (See Figure 2). As proven in Proposition 2,
Figure 2. Numerical simulations of as a function of the selected model parameters. (a) values increase with provided values are small. When values are large, decreases with ; (b) increases both with and ; (c) increases with and decreases with ; (d) increases linearly with and increases parabolically with .
Figure 2(a) shows that is an increasing function with respect to . 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 is an increasing function with respect to and a decreasing function with respect to . Figure 2(a) and Figure 2(b) demonstrate the behavior of with respect to the influx of infected deer, . For smaller values of or , is an increasing function with respect to ; for larger values of or , it is a decreasing function with respect to . Thus, there must be a critical value ( or ) where the behavior changes.
If we consider as a function of the deer-midge interactions, then is essentially a linear function of and a function of the square root of . The graph of would be increasing and concave down with respect to an increase in (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 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 and achieves a disease free equilibrium if . To do this, a MATLAB code was written utilizing the ODE45 solver, and the results were verified against the computed value for a given set of parameters. At time , we have the following initial values: , , , , , and . 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 values.
Note: The 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 , 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 , the system stabilizes to its endemic equilibrium. See Table 2 for the specific values used and the corresponding values of .
Figure 3(a) and Figure 3(b) indicate that when , the system will stabilize to its disease free equilibrium. Figure 3(c) and Figure 3(d) show that when , the system will stabilize to an endemic equilibrium. These outcomes are robust for large sets of initial values and parameter values.
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., and ). The results show that there are critical values for the interaction parameters and rates of susceptible deer dispersion . 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  and Weisser et al.  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  . 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  . We also note that weather has an effect on both the midge population and the life cycle of the HD virus   . 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  . In addition, record warm temperatures in that spring and summer may cause midges to become more active sooner than normal  . 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  . 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.
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 expression provides insights into the effects of deer movement on the spread of disease.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.