Modeling Gender-Structured Wildlife Diseases with Harvesting: Chronic Wasting Disease as an Example

Chronic wasting disease (CWD) is a prion infectious disease that affects members of the deer family in North America. Concerns about the economic consequences of the presence of CWD have led management agencies to seek effective strategies to control CWD distribution and prevalence. Current mathematical models are either based on complex simulations or overly simplified compartmental models. We develop a mathematical model that includes gender structure to describe CWD in a logistically growing population. The model includes harvesting as a management strategy for the disease. We determine the stability conditions of the disease-free equilibrium for the model and calculate the basic reproduction number. We find an optimum interval of harvesting: with too little harvesting, the disease persists, whereas too much harvesting results in extinction of the population. A sensitivity analysis shows that the disease threshold is more sensitive to female than male harvesting and that harvesting has the greatest effect on the basic reproduction number. However, while harvesting may be a way to control CWD, the range of admissible harvesting rates may be very narrow, depending on other parameters.


Introduction
Many wild animals come into contact with humans via harvesting, for meat, for pelts, or for population control.Transmissible diseases in wild populations pose a threat to ecology, the environment, and humans.Harvesting may also be a way that such disease transmission can be controlled; however, accurate harvesting of only sick animals is usually impractical or impossible.For example, many diseases have an extended latency period during which the animal is asymptomatic; as a result, targeted harvesting is not feasible.When targeted harvesting is impossible, harvesting for the sake of disease control has to be tightly balanced to ensure population survival.Mathematical models are important tools to determine optimal harvesting rates.
A mathematical model should be as simple as possible to allow for easy parameterization and robust results; on the other hand, models need to include all important mechanisms of disease transmission and population behaviour.For example, if the latent period is comparable in duration to the life cycle of the animal, then incorporating population dynamics into a model becomes critical.Also, vertical transmission can sustain many diseases; in such cases, infection is passed from mother to offspring.Consequently, gender structure of the herd may play a crucial role in the life cycle of the disease.Conversely, many animals display distinct gender behaviour such as herding or fighting, which may affect horizontal transmission.
Chronic wasting disease (CWD) is a prion infectious disease [1] that affects members of the deer family, including mule deer, white-tailed deer, elk, and moose in North America [2,3].It is associated with protein changes in the brain, which lead to changes in behaviour, mood, and metabolism.These changes include excessive salivation, grinding of teeth, increased drinking and urination, dramatic loss of weight and body condition, poor hair coat, and staggering.Individuals cannot maintain weight and slowly waste away.CWD belongs to the same family of diseases as bovine spongiform encephalopathy (BSE), scrapie and Creutzfeldt-Jacob disease.
There is no known cure or treatment for CWD and all infected individuals eventually die [4].CWD can be passed vertically from mother to offspring and horizontally upon direct contact as well as indirectly through contamination of water and feed by saliva and/or feces [5][6][7].CWD seems more likely to occur where elk or deer are crowded or where they congregate at human-made feed and water stations [8].Concerns about the economic consequences of the presence of CWD have led management agencies to seek effective strategies to control CWD distribution and prevalence [9,10].
The relationship between the disease and the ecology of wild cervids is complex.Horizontal transmission depends on the social structure of the species; the incubation period is relatively long (one to three years [11]) compared to gestation, maturation, and even survival; and lifespan and harvesting rates depend on gender.
To the best of our knowledge, previous modeling approaches to CWD consist of either highly complex simulation models with tens if not hundreds of different classes [12][13][14] or simple differential equation models of four compartments or fewer [7,15].In particular, there is no genderstructured CWD model that is analytically tractable for general insights and observations.In the absence of disease, a simple age-and-gender model for cervids was successfully analyzed in [16].
Here, we develop an analytically tractable model to describe the spread of an infectious disease in a wild animal population with harvesting.We use data from the spread of CWD in a cervid population to illustrate the effects of horizontal and vertical transmission, gender, latency, and asymptomatic infection.We address the following research questions: (1) what is the role of gender structure in the spread of the disease?(2) is it possible to use harvesting to control the disease without eradicating the population?and (3) how does the spread of the disease depend on parameters?
This paper is organized as follows.In Section 2, we derive our models.In Section 3, we analyze a model of latency with constant birth and linearly increasing death rate.In Section 4, we extend the model to incorporate an asymptomatic class.In Section 5, we incorporate gender structure.In Section 6, we apply our results to CWD.We conclude with a discussion.

Model Derivations
A standard population dynamics model that accounts for reproduction, death, and harvesting is the logistic equation with a linear harvesting term (e.g., [17]).The logistic equation does not specify births and deaths separately, but rather their net effect.When the population is structured into different compartments such as gender or disease state, we need to specify births into and deaths from each class separately, since their values may depend on the state of the individual.Hence we start with a general description of the dynamics of the total population N: where B and D are the per capita birth and death rates, respectively, and h is the harvesting rate.We can recover the logistic equation by setting the birth rate constant and the death rate proportional to population density.In some species, birth rates may decrease with density whereas death rates are relatively unaffected [18].In general, as long as the difference B(N) − D(N) is a monotone decreasing function, (1) will exhibit logistic-like behaviour.Next, we split the population into compartments according to disease state as susceptible (S), exposed (E), and infectious (I).We assume that all individuals are born susceptible; that is, there is no vertical transmission.Infectious individuals experience an additional disease-induced mortality μ d .Disease transmission is by mass action with rate β; individuals progress from the exposed to the infectious state with rate σ.In particular, 1/σ is the average duration of the exposed period.Then we obtain the system of equations where N = S + E + I.If infectious individuals are unable to reproduce, then the birth term in the S-equation reads B(N)(S + E).We did not include an explicit harvesting term in the infectious class.If infectious individuals can be distinguished easily from noninfectious ones, a simple control strategy would be to cull only these individuals.In that case, we interpret μ d as the combined effect of disease-induced death and disease-specific culling.When the disease is not easily detected, then we can replace μ d with μ d + h in the model and the analysis remains the same.In CWD, the disease becomes only apparent in the final stage where the disease-induced mortality is much higher than any normal harvesting rate, in which case μ d + h ≈ μ d .For this reason, we do not explicitly include harvesting in the infectious class.
A crucial aspect of many diseases, and certainly of CWD, is that an exposed individual can transmit the disease long before clinical symptoms become apparent.In fact, the final, clinical stage of CWD is very short compared to the entire duration of the disease.We include such an infectious but nonclinical stage into our model as asymptomatic (A).We allow for differential transmission rates from the asymptomatic and infectious class (β A , β I ) and denote the sojourn times in the exposed and asymptomatic class as 1/σ E and 1/σ A , respectively.Finally, we assume that individuals in the asymptomatic class can have offspring, and that vertical transmission occurs with probability 1 − η.Then the equations read ( Finally, we introduce gender structure into our equations.This inclusion does not only increase the number of The full model.Animals can be susceptible, exposed, asymptomatic, or infected.We assume offspring of infected animals do not survive and that some proportion of offspring of asymptomatic animals will be susceptible, the remainder exposed.
compartments but it also requires a more detailed examination of the birth rate.We begin again in the absence of disease and split the total population into two gender classes where N f is the female population and N m the male population.The question of how the per capita birth rate depends on the size of the male and female population is studied in pair-formation models [19].In general, we can write the equations for N f and N m as where, for simplicity, we assume a 1 : 1 offspring ratio.One possible choice of how the birth function depends on the female-to-male ratio is In this case, b is the fecundity, which is half the maximum per-female number of offspring at low densities, and k −1 is the average harem size [20].We consider a different choice of birth function in Section 7.
When we include disease classes in a gender-structured model, contact rates γ q (q ∈ {m, f , f m, m f }) and transmissibility β i j ( j ∈ {A, I}) can depend on gender.For simplicity, we assume that infectious individuals cannot reproduce and denote the effective population size for reproduction as N f = S f +E f +A f and N m , accordingly.This assumption is justified for CWD.We have the following system of equations: Here, the force of infection is given by The model is illustrated in Figure 1, with variables and parameters listed in Table 1.

The Basic Model
We consider the simplest model (2) with constant birth rate and linearly increasing death rate.We determine the equilibria, their stability, and the nature of the arising bifurcations.We give several equivalent characterizations of the disease threshold R 0 , also known as the basic reproduction number [21,22].The model reads Average harem size η The probability that the newborn of an asymptomatic mother is susceptible where α measures the strength of density dependence in the death rate.

The Equilibria.
We set the right-hand side of ( 8)- (10) equal to zero and reduce the nonlinear system to a single polynomial of degree four in the total population density N.
Adding the steady state equations for S and E, we obtain Here, we use S e , E e , I e , and N e to represent the equilibrium values.
Substituting S + E = N − I and solving for E in terms of I and N from (10) with I = 0, we obtain We can express the steady-state values of S e , E e , and I e in terms of N e as under the condition b − h + μ d / = 0. Using ( 9) with E = 0, we obtain where and There are solutions N = 0, corresponding to the absence of population, and S = N = K, E = I = 0, corresponding to a disease-free equilibrium (DFE) at carrying capacity.For this steady state to be biologically meaningful and for the population to survive (K / = 0), we require Finally, we obtain another solution 0 < N e < K, S e , E e , I e > 0, where This solution is positive under the condition b − h + μ d > 0.

Stability.
The zero equilibrium of system ( 8)-( 10) is unstable if condition ( 16) holds; that is, b > h + μ.We start by analyzing the stability of the DFE (K, 0, 0).The corresponding Jacobian matrix has one negative eigenvalue λ = −b + μ + h (see (16)) and the two eigenvalues of the matrix Since all parameters are positive, the trace of J 1 is negative, due to (16).If the determinant of J 1 is positive, then the eigenvalues of J 1 are negative.Equivalently, we require Therefore, we have proved the following proposition.
Note that ( 20) is equivalent to where K C is the critical carrying capacity.As a result, we have the following corollary.
Corollary 2. The equilibrium (K, 0, 0) is stable for K < K C and unstable for K > K C .
This observation that the disease will die out for small populations is the basis for using harvesting (of healthy individuals) to control the disease.We explore this management option in Section 6.
Next, we would like to show that the endemic equilibrium (EE) is positive and locally asymptotically stable exactly when the DFE is unstable; that is, a transcritical bifurcation occurs.We first focus on the transmission rate β.
Solving N e = K, we obtain the equation This implies from which we conclude Thus Note also that, at this critical rate of transmission, we have E e = I e = 0 and S e = K; see (13).
As a result, we have the following two corollaries.
(1) The EE exists if and only if R 0 > 1.
(2) Note that β > β C implies N e < K and (K, 0, 0) is unstable, which means the population either goes to extinction or the disease persists.So if ( 16) holds, then we can avoid extinction because (0, 0, 0) is unstable in this case.
We now employ centre manifold theory (see, e.g., [23]) to analyze the local stability of the EE.
Theorem 5.For R 0 > 1, but close to 1, the unique EE is locally asymptotically stable.
Proof.We utilize Theorem 4 in [24].First, note that the Jacobian matrix for the DFE has eigenvalues λ = −b + μ + h, −b − σ and 0, when β = β C .Also, the zero eigenvalue has right eigenvector (u 1 , u 2 , u 3 ) and left eigenvector (v 1 , v 2 , v 3 ), where and u 3 > 0 is free.Also, and v 3 > 0 is free.We have where E 0 is the DFE.This implies that (i) or (iv) in Theorem 4 [24] is applicable.However, for β < β C , E 0 is locally asymptotically stable and we have no other positive equilibria.As a result, (iv) is the only applicable case.This means that, when β changes from β < β C to β > β C , E 0 changes from stable to unstable, and the EE changes from negative to positive and locally asymptotically stable.
To summarize, we have analyzed the basic SEI model with constant birth and linearly increasing death rate.We determined the equilibria and found the conditions that ensured the existence of the EE.We also calculated the basic reproduction number, a critical carrying capacity, and a critical rate of transmission.Finally, we used each of these three values to describe the stability of the DFE and the EE.

The Basic Model with Asymptomatic Individuals
As in the previous section, we choose a constant birth rate and a linearly increasing death rate to analyze this model.Hence the equations are 4.1.The Equilibria.The equilibria of ( 28)- (31) satisfy Equation ( 29) with with ISRN Biomathematics 7 Therefore, the equilibrium solutions are N = 0, corresponding to the absence of the population, and N = K = S, E = A = I = 0, corresponding to the DFE when ( 16) holds.Finally, when we have 0 < N e < K, S e , E e , A e , I e > 0.
We define the polynomial (1) P(0) < 0 for η close to 1, and (2) Therefore, (39) has at least one positive root for η close to 1 (i.e., when there is not too much vertical transmission).

Stability.
The zero equilibrium is unstable when condition ( 16) holds; that is, b > h + μ.For the DFE (K, 0, 0, 0), the Jacobian matrix has the negative eigenvalue λ = −b + μ + h (see (16)) and three eigenvalues of the matrix which are the zeros of the characteristic polynomial λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0, with coefficients Note that if a 0 = 0, then a 1 > 0, which means that there are two negative eigenvalues and one zero eigenvalue.Also, a 0 = 0 implies In order to have three negative eigenvalues, we require a 0 > 0, which is equivalent to R 0 < 1.Therefore, we have proved the following proposition.
Condition (45) is equivalent to the critical carrying capacity.As a result, we have the following corollary.
Corollary 7. The equilibrium (K, 0, 0, 0) is stable for K < K C and unstable for K > K C .
Remark 8.Note that if β A → 0 and σ A → ∞ (the probability of the disease transmission from asymptomatic animals and the duration of the asymptomatic period both tend to zero), then (45) and ( 20) are equivalent.
Substituting N e = K in ( 29) with E = 0 and solving for

ISRN Biomathematics
We have thus proved that Moreover, at this critical rate of transmission, we have E e = A e = I e = 0 and S e = K.
As a result, we have the following two corollaries.
(1) For R 0 > 1 and η → 1, the EE exists.Since P(K) > 0 and P(0) < 0, the result follows from applying the intermediate value theorem on the polynomial P in [0, K].However, for R 0 < 1, the existence of the EE is not clear.
(2) The condition β > β C implies N e < K, which means the population is driven to extinction or the disease persists.If (16) holds, then we can avoid the extinction case because (0, 0, 0, 0) is unstable.The stability of the EE is an open problem.
(3) Note that R 0 is a decreasing function of η while β C is an increasing function of η.Thus a vaccine that reduces vertical transmission may help control CWD.
In this section, we extended the basic model to include asymptomatic individuals.We calculated the equilibria and found conditions that ensure the existence of the EE.We also calculated the basic reproduction number, a critical carrying capacity, and a critical rate of transmission; these three values are used to describe the stability of both the DFE and the EE.

The Full Model
We now analyze stability in the full, gender-structured model. Let The equilibria thus satisfy We have two cases.
Case 1.If S f = S m = 0, we have the zero equilibrium.
Case 2. The nonzero disease-free equilibrium is calculated as follows.First, we multiply (49) by μ m /S f and (50) by μ f /S m .Then we solve the resulting system for z ≡ S m /S f and K = S f + S m .We obtain Equality of the two expressions for K results from the solution for z.For the equilibrium densities for male and females, we obtain A positive equilibrium requires K > 0, which implies To investigate the stability of the disease-free equilibrium, we use the method of the next-generation operator [25,26].Set The Jacobian matrix for the equilibrium (E f , E m , I f , I m , A f , A m , S f , S m ) = (0, 0, 0, 0, 0, 0, K/(1 + z), zK/(1 + z)) has the form where The 2 × 6 matrix J 1 turns out to be irrelevant to our calculations.We have The threshold parameter R 0 is the spectral radius of the matrix FV −1 , provided that all eigenvalues of J 2 have negative real part [26].Thus where (60)

Note that if b
In order to prove that J 2 has only eigenvalues with negative real part, we prove that tr(J 2 ) < 0 and det(J 2 ) > 0. First, we note the inequalities From (53), we obtain tr(J 2 ) < 0. Now, using (53), (61), as well as the fact that z and all parameters are positive, we have As a result, we have proved the following theorem.
Theorem 11.The DFE is stable for R 0 < 1 and unstable for R 0 > 1.
Note that the value of R 0 derived here is a bifurcation parameter, not necessarily the average number of secondary infections [22].Furthermore, we have illustrated local stability of the DFE, but cannot rule out the possibility of backward bifurcations.
Remark 12.Note that, for the full model, we have the following.Let where 3 , Moreover, if γ f = γ f m and γ m = γ m f , then R 0 = 1/2a 1 .If k = 0 and all parameters are identical for males and females, then z = 1 and which is the value of R 0 for the basic model (20).In the same way, if k = 0 and all parameters for females and males are identical, then we will have z = 1 and R 0 for the full model (58) matches R 0 for the basic model with asymptomatic individuals (45).
In this section, we have analyzed the full model.We calculated the DFE and the basic reproduction number and determined the stability of the DFE.We also related the basic reproduction number for the full model to the basic reproduction number for the two special cases.

Application to CWD
For the purpose of simulations, we used Wildlife Management Unit 150-Alberta (WMU 150), which has an area equal to 1841 km 2 [27].This area is not isolated and contains different kind of deer, so the simulation is a lower estimate for the number of infected animals in this area.The values of parameters used in the simulations are given in Table 2.
We used the following information.
(1) We consider harvesting at the same rate as in Alberta in the year 2000.Out of 120,000 mule deer, 8,256 does and 8,586 bucks were harvested.Also, the ratio of number of bucks to does is 34 : 100 (range 5-62 : 100) [29].As a result, the average numbers are 30,448 bucks and 89,552 does.The harvesting rates are estimated by h f = 0.0922, h m = 0.282, and h = 0.1512 (estimated using P(t) = P 0 e −ht , where P(t) is the population after t years and P 0 is the initial population).
(3) Following a widespread culling program, WMU 150 had 650 mule deer in 2006 [27].Among 245 adult animals tested during this program, there were 9 infected cases in early stages (mostly around the perimeter of the area).We estimate the number of cases in the early stages in the remaining 650 as E 0 = 24.We assume that there are no infected animals with clinical signs (I 0 f = I 0 m = 0).All other animals are considered susceptible (S 0 f + S 0 m = 626).Also, we assume that the disease prevalence in males is initially twice that in females [32], so we have E 0 f = 8 and E 0 m = 16.
(4) It is extremely difficult to obtain reasonable estimates or even for the disease transmission rates.Male fighting could increase their contact rates significantly, but injuries during such fights are rare.Contact rates could be higher for females as they live in large herds and have grooming contacts.Contact rates may also depend on season; for example, when animals gather around feeding areas in the winter.Since no clear data are available, we choose a model situation with minimal gender differences.
We assume that contact rates between the different compartments are identical (i.e., there is only one value of γ), but the transmission rates differ between females and males.We set β m , and β f = β − 0.01, β m = β + 0.01.See Table 2.
For the basic model, the endemic equilibrium (S e , E e , I e ) = (100, 56, 17) is stable and the basic reproduction number is R 0 = 5.407828.For the basic model with asymptomatic individuals, the endemic equilibrium (S e , E e , A e , I e ) = (47, 27,12,11) is stable and the basic reproduction number is R 0 = 7.026539.Finally, for the two-gender nonasymptomatic model, with η = 1, σ A → ∞ and β A → 0, the endemic equilibrium (S f , S m , E f , E m , I f , I m ) = (103, 73, 15, 8, 4, 2) is stable and the basic reproduction number is R 0 = 1.56363.These numbers indicate that the disease persists in the absence of control measures.Note how prevalence and R 0 are significantly reduced in the full model, compared to the simpler models.We return to this point in Section 7.
We examine harvesting as a method of controlling the disease.Figure 2 shows the graph of R 0 and K as functions of harvesting for the basic model.Any harvesting rate between h 1 = 0.4130 and h 2 = 0.4647 implies that the disease goes to extinction, but the population will survive.In Figure 3, we take three values for harvesting (h = 0.1512, 0.45, 0.55) and track the deer population in WMU 150 for 50 years.
Figure 4 shows the graph of R 0 and K as functions of harvesting for the basic model with asymptomatic individuals.Any harvesting rate between h 1 = 0.4262 and h 2 = 0.4647 implies that the disease goes extinct, but the population survives.In Figure 5, we take three values for harvesting (h = 0.1512, 0.45, 0.55) and track the deer population in WMU 150 for 50 years.
Figure 6 shows the level curves R 0 = 1 and K = 0 as functions of gender-specific harvesting parameters h f and h m for the full model.Any point (h f , h m ) between R 0 = 1 and K = 0 implies that the disease is controlled and the population survives.In Figure 7, we take three different points ((h f , h m ) = (0.0922, 0.282), (0.4, 0), (0.6, 0)) and track the deer population in WMU 150 for 50 years.Note also that the equilibria for the female population are different from the equilibria for the male population.6.1.Sensitivity Analysis.Due to the degree of uncertainty in the parameters values, we considered a range of parameters to examine the dependence of R 0 on parameter variation for the basic model.We used Latin hypercube sampling and partial rank correlation coefficients (PRCCs) to identify which parameters R 0 is most sensitive to [33].Latin hypercube sampling is a statistical sampling method that evaluates sensitivity of an outcome variable to all input variables.PRCCs measure the relative degree of sensitivity to each parameter, regardless of whether the parameter has a positive or negative influence on the outcome variable.Figure 8 plots PRCCs for each input parameter.This demonstrates that R 0 is most sensitive to variations in harvesting rates.Figure 9 shows that the disease is reliably controlled only for high harvesting rates or if transmissibility is sufficiently small.
It should be noted that variations in h will change R 0 from values much greater than 1 to very low values, resulting in significant dependence of R 0 on h in Figure 8.However, sufficiently high values of h result in population extinction, so we constrained h to avoid such harvesting.Hence the effect of realistic values of h in Figure 9 is smaller than those predicted by Figure 8.It follows that β may have a larger effect on the outcome in practice, but we have little control over the transmissibility.

Discussion
We developed a model to evaluate the effects of birth and death dynamics, gender structure, and disease transmission on the fate of the disease and the population.Despite the model complexity, we found explicit expressions for the basic  reproductive number and thereby the stability of the diseasefree equilibrium.For some model simplifications, we also found explicit stability conditions for the endemic equilibrium.By building models of increasing complexity, we gained a crucial understanding of important parameters and mechanisms from comparing the simpler models.
We modelled logistic growth as a density-dependent death rate (with parameter α i ), whereas some sources in the literature suggest that birth, rather than death, is density dependent in large herbivores [34].We chose a densitydependent death term since we were able to obtain more detailed analytical results this way.We also considered density-independent death combined with density-dependent birth; that is, we set α i = 0 and multiplied the birth terms by a factor (1 + κN) −1 .In the simplest case of a basic model, this approach gives a disease-free equilibrium with and the basic reproduction number As before, we find an interval of harvesting rates that guarantees persistence of the population and extinction of the disease; see Figure 10.We expect that all the results presented in this work hold qualitatively also for the case of density-dependent birth functions; however, the algebraic equations are much harder to solve.The existence of a population-level threshold for the outbreak or control of a disease depends on whether one chooses to model disease transmission via density or frequency dependence [14,35].In particular, such a threshold typically does not exist in simple models with frequencydependent transmission.While density-dependent transmission, as chosen here, is often the first choice for human disease dynamics, evidence for either modeling assumption is hard to get in wildlife populations.Conner et al. state that the question is open [36].Wasserberg et al. found a slight preference for frequency-dependent transmission in a CWD dataset from South-Central Wisconsin [14], but an acceptable fit to data with density-dependent transmission as well.Frequency-dependent transmission has the advantage that it avoids the assumption of mass-action transmission and describes populations with complex structure [37].An increase in the transmission rate with population density is to be expected in wild populations during the winter when animals gather around (human-created) feed stations.Density-dependent transmission is also typical in small populations.Hence when harvesting reduces population size, then density dependence should arise.Future work could explore transmission rates that vary between density and frequency dependence as population size varies [38].
The model is motivated by the study of CWD dynamics in wild cervid populations.Comparing the simpler models with the full model can give significant insights into the relative importance of model assumptions.Parameterizing the models with data from WMU 150 in Alberta, we highlighted the importance of considering the asymptomatic class and, even more so, gender structure.While the value of R 0 calculated for the models without gender structure is relatively high (in fact, most likely too high), the value for R 0 in the gender-structured model is quite reasonable.Accordingly, the prevalence of CWD in the simpler models is very high, but the corresponding numbers for the genderstructured case are reasonably low.In a study in 2000, Miller et al. found that prevalence of CWD in Colorado and Wyoming ranged between 1% and 5%, depending on species [13].for a variant of the model that considers density-dependent births and density-independent deaths.The disease persists when R 0 > 1, and the population persists when the harvesting rate is less than the larger of the two vertical lines.In between the two vertical lines, the population persists but the disease dies out.
We investigated the possibility of disease control through harvesting.We found that there are two harvesting thresholds.When harvesting is too low, then the disease remains endemic; when harvesting is too high, then the population declines to extinction.Alternatively, these thresholds can be expressed in terms of the critical carrying capacity or the critical transmission rates.Furthermore, we note that a vaccine that controlled vertical transmission would reduce R 0 and increase the critical transmission rate.This reduction in R 0 may be enough to control the disease, depending on other parameters.
Note that the harvesting window is small for the parameters considered here.Furthermore, due to the long latency and asymptomatic classes, harvesting healthy animals is unavoidable.While that may be desirable in some cases, for disease control the result is a high burden on the population.It is thus crucial that the ecological impact is carefully assessed before harvesting is used to control the spread of disease.
In our numerical simulations, we saw that R 0 is most sensitive to the harvesting rate.Using data from Alberta, we found an admissible interval for harvesting rates of male and female deer to ensure disease control and population survival.Too much harvesting and the population will die out; too little and the disease persists.This illustrates the importance of strict harvesting guidelines.A more detailed analysis of the full model also demonstrates that the outcome is significantly more dependent on harvesting of females than males (Figure 7).This result is to be expected, since females contribute more strongly to population growth.
Finally, we note that applying the ODE framework to populations of extremely small populations is limited; although our numerical simulations illustrate small populations, it should be noted that the assumptions of such models may break down at this scale.In particular, stochasticity in small populations may overwhelm the deterministic component; such fluctuations have the potential to drive small populations to extinction.
One open mathematical question is the stability of the endemic equilibrium when R 0 is greater than one.Our simulations suggest that the situation is the same as in the simpler models, but we cannot exclude the possibility of oscillating solutions or a backward bifurcation at R 0 = 1 [22].Two major challenges for future research are to include environmental disease transmission and to consider seasonality in the model.The importance of environmental transmission results from extremely long persistence times of prions in the environment.Abandoned cervid farms could act as reservoirs for prions.Migrating herds could transmit the disease even without physical contact with infected resident herds.Since the social behaviour of wild cervids shows strong seasonal variation, disease transmission is likely to vary within a year.In addition, reproduction in wild cervids is strongly seasonal.A combination of discrete year-to-year dynamics with a continuous seasonal transmission and survival model could be the next step in that direction.

1 Figure 2 :
Figure 2: Critical harvesting rates in the basic model.Too little harvesting results in disease persistence, while too much eradicates the population.

Figure 3 : 1 Figure 4 :Figure 5 :
Figure 3: Deer population in the WMU 150 area for the basic model when (a) h < h 1 , resulting in disease persistence; (b) h 1 < h < h 2 , resulting in disease control; (c) h > h 2 , resulting in eradication of the population (solid curve; susceptible; diamonds: exposed; dashed curve: infected).

Figure 6 :
Figure 6: Critical harvesting rates in the full model.Too little harvesting results in disease persistence.Too much eradicates the population.Note that the outcome varies with harvesting of females much more than with male harvesting.

Figure 7 :Figure 10 :
Figure 7: Deer population in the WMU 150 area for the full model.((a) and (b)): If h < h 1 , then the disease persists, although it is stronger in males.((c) and (d)): If h 1 < h < h 2 ,then the disease can be controlled and the population survives.In this case, only females are harvested, resulting in greater survival of males.((e) and (f)): if h > h 2 , then the population is eradicated.In this case, only females are harvested, so their numbers decline quicker (solid curve; susceptible; diamonds: exposed; dashed curve: infected).
m Density-dependent reduction rate for males k−1

Table 2 :
Sample values and sensitivity ranges for mule deer.