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

Authors: Mo'tassem Al-Arydah; Robert J. Smith?; Frithjof Lutscher
Journal: ISRN Biomathematics
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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–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–14] or simple differential equation models of four compartments or fewer [7, 15]. In particular, there is no gender-structured 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 Section2, 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.


## 2. 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:
(1)N′=B(N)N-D(N)N-hN,
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
(2)S′=B(N)N-D(N)S-hS-βIS,E′=-D(N)E+βIS-(σ+h)E,I′=-D(N)I+σE-μdI,
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
(3)S′=B(N)(S+E+ηA)-D(N)S-hS-(βAA+βII)S,E′=B(N)(1-η)A-D(N)E+(βAA+βII)S-(σE+h)E,A′=σEE-D(N)A-(σA+h)A,I′=-D(N)I+σAA-μdI.Finally, we introduce gender structure into our equations. This inclusion does not only increase the number of 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 classesN=Nf+Nm, where Nf is the female population and Nm 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 Nf and Nm as
(4)Nf′(t)=B(Nf,Nm)Nf-D(N)Nf-hfNf,Nm′(t)=B(Nf,Nm)Nm-D(N)Nm-hmNm,
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
(5)B(Nf,Nm)=bNmkNf+Nm.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,fm,mf}) and transmissibility βji(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=Sf+Ef+Af and N^m, accordingly. This assumption is justified for CWD. We have the following system of equations:
(6)Sf′(t)=B(N^f,N^m)(Sf+Ef+ηAf)-D(N)Sf-hfSf-FfSf,Sm′(t)=B(N^f,N^m)(Sf+Ef+ηAf)-D(N)Sm-hmSm-FmSm,Ef′(t)=B(N^f,N^m)(1-η)Af-σEEf-D(N)Ef-hfEf+FfSf,Em′(t)=B(N^f,N^m)(1-η)Af-σEEm-D(N)Em-hmEm+FmSm,Af′(t)=σEEf-σAAf-D(N)Af-hfAf,Am′(t)=σEEm-σAAm-D(N)Am-hmAm,If′(t)=σAAf-D(N)If-μdIf,Im′(t)=σAAm-D(N)Im-μdIm.
Here, the force of infection is given by
(7)Ff=βAf(γfAf+γmfAm)+βIf(γfIf+γmfIm),Fm=βAm(γfmAf+γmAm)+βIm(γfmIf+γmIm).
The model is illustrated in Figure 1, with variables and parameters listed in Table 1.Table 1
Parameter definitions.
Variable
Definition
S
f
(
t
)
Number of susceptible female deer at timet
S
m
(
t
)
Number of susceptible male deer at timet
E
f
(
t
)
Number of exposed female deer at timet
E
m
(
t
)
Number of exposed male deer at timet
I
f
(
t
)
Number of infected female deer at timet
I
m
(
t
)
Number of infected male deer at timet
b
Fecundity
β
Transmission rate
h
f
Harvesting rate of females
h
f
Harvesting rate of males
β
A
f
Transmissibility to females from an asymptomatic deer
β
I
f
Transmissibility to females from infected deer
β
A
m
Transmissibility to males from asymptomatic deer
β
I
m
Transmissibility to males from infected deer
β
A
Rate of transmission from asymptomatic deer
β
I
Rate of transmission from infected deer
γ
m
f
Rate of contact between infected males and susceptible females
γ
f
m
Rate of contact between infected females and susceptible males
γ
m
Rate of contact between infected males and susceptible males
γ
f
Rate of contact between infected females and susceptible females
L
Incubation period
L
E
Lifespan for exposed deer
L
A
Lifespan for asymptomatic deer
σ
Rate of progression from exposed class to infected class
σ
E
Rate of progression from exposed class to asymptomatic class
σ
A
Rate of progression from asymptomatic class to infected class
μ
f
Natural death rate for females
μ
m
Natural death rate for males
μ
d
Death rate due to disease
α
f
-
1
Density-dependent reduction rate for females
α
m
-
1
Density-dependent reduction rate for males
k
-
1
Average harem size
η
The probability that the newborn of an asymptomatic mother is susceptibleFigure 1
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.


## 3. 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 R0, also known as the basic reproduction number [21, 22]. The model reads
(8)S′(t)=b(S+E)-μ(1+αN)S-hS-βIS,(9)E′(t)=-μ(1+αN)E+βIS-(σ+h)E,(10)I′(t)=-μ(1+αN)I+σE-μdI,
where α measures the strength of density dependence in the death rate.


### 3.1. 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
(11)(b-h)(Se+Ee)-μ(1+αNe)(Se+Ee)-σEe=0.
Here, we use Se, Ee, Ie, and Ne to represent the equilibrium values.SubstitutingS+E=N-I and solving for E in terms of I and N from (10) with I′=0, we obtain
(12)(b-h)Ne-μ(1+αNe)Ne-(b-h+μd)Ie=0.
We can express the steady-state values of Se, Ee,  and  Ie in terms of Ne as
(13)Ie=(b-h-μ(1+αNe))Neb-h+μd,Ee=(b-h-μ(1+αNe))(μd+μ(1+αNe))Neσ(b-h+μd),Se=Ne-Ie-Ee=Ne(μd+μ(1+αNe))(σ-b+h+μ(1+αNe))σ(b-h+μd),
under the condition b-h+μd≠0. Using (9) with E′=0, we obtain
(14)Ne(Ne-K)(n2Ne2+n1Ne-n0)=0,
where K=(1/μα)(b-h-μ) and
(15)n2=βμα,n1=((σ-b+h+μ)β-μα(b-h+μd)),n0=(b-h+μd)(σ+h+μ)>0.
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
(16)b-h-μ>0.
Finally, we obtain another solution 0<Ne<K, Se,Ee,Ie>0, where
(17)Ne=12n2(-n1+n12+4n2n0).
This solution is positive under the condition b-h+μd>0.


### 3.2. 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
(18)J(K,0,0)=[-b+μ+hh+μ-b+μ+h-βμα(b-μ-h)0-b-σβμα(b-μ-h)0σ-b+h-μd]
has one negative eigenvalue λ=-b+μ+h (see (16)) and the two eigenvalues of the matrix
(19)J1=[-b-σβμα(b-μ-h)σ-b+h-μd].
Since all parameters are positive, the trace of J1 is negative, due to (16). If the determinant of J1 is positive, then the eigenvalues of J1 are negative. Equivalently, we require
(20)R0≡βσ(b-μ-h)μα(σ+b)(b-h+μd)<1.Therefore, we have proved the following proposition.Proposition 1.
The equilibrium(K,0,0) is locally asymptotically stable for R0<1 and unstable for R0>1.Note that (20) is equivalent to
(21)K<(σ+b)(b-h+μd)βσ≡KC,
where KC is the critical carrying capacity. As a result, we have the following corollary.Corollary 2.
The equilibrium(K,0,0) is stable for K<KC and unstable for K>KC.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 Section6.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β.SolvingNe=K, we obtain the equation (22)1μα(b-h-μ)=12βCμα{[(μα(b-h+μd)-βC(σ-b+h+μ)]2+4βCμα(b-h+μd)(σ+h+μ)μα(b-h+μd)-βC(σ-b+h+μ)+[(μα(b-h+μd)-βC(σ-b+h+μ)]2+4βCμα(b-h+μd)(σ+h+μ)}.This implies
(23)σ(b-h-μ)=μαβC(b-h+μd)(b+σ),
from which we conclude
(24)β=μα(b-h+μd)(b+σ)σ(b-h-μ)≡βC.
Thus β=βC if and only if Ne=K. Also, it is clear that β=βC if and only if R0=1. Note also that, at this critical rate of transmission, we have Ee=Ie=0 and Se=K; see (13).As a result, we have the following two corollaries.Corollary 3.
Suppose0<Ne≤K. Then the following are equivalent: (i) β=βC, (ii) Ne=K, (iii) K=KC, (iv) R0=1 and (v) (Se,Ee,Ie)=(K,0,0).Corollary 4.
Suppose0<Ne≤K holds. Then(1)
forβ<βC, we have  R0<1 so the disease-free equilibrium  (K,0,0) is locally stable;(2)
forβ>βC, we have  R0>1 so the disease-free equilibrium  (K,0,0) is unstable.Remarks.
Suppose (16) holds.(1)
The EE exists if and only ifR0>1.(2)
Note thatβ>βC implies Ne<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.
ForR0>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 (u1,u2,u3) and left eigenvector (v1,v2,v3), where
(25)u1=(1-b-h-μdσ-b-h-μdb-μ-h)u3,u2=(b-h+μd)σu3
and u3>0 is free. Also,
(26)v1=0,v2=σb+σv3
and v3>0 is free. We have
(27)∑k,i=1n‍vkui∂2fk∂xi∂β(E0,βC)=Kv2u3>0,where E0 is the DFE. This implies that (i) or (iv) in Theorem 4 [24] is applicable. However, for β<βC, E0 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, E0 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.


## 3.1. 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
(11)(b-h)(Se+Ee)-μ(1+αNe)(Se+Ee)-σEe=0.
Here, we use Se, Ee, Ie, and Ne to represent the equilibrium values.SubstitutingS+E=N-I and solving for E in terms of I and N from (10) with I′=0, we obtain
(12)(b-h)Ne-μ(1+αNe)Ne-(b-h+μd)Ie=0.
We can express the steady-state values of Se, Ee,  and  Ie in terms of Ne as
(13)Ie=(b-h-μ(1+αNe))Neb-h+μd,Ee=(b-h-μ(1+αNe))(μd+μ(1+αNe))Neσ(b-h+μd),Se=Ne-Ie-Ee=Ne(μd+μ(1+αNe))(σ-b+h+μ(1+αNe))σ(b-h+μd),
under the condition b-h+μd≠0. Using (9) with E′=0, we obtain
(14)Ne(Ne-K)(n2Ne2+n1Ne-n0)=0,
where K=(1/μα)(b-h-μ) and
(15)n2=βμα,n1=((σ-b+h+μ)β-μα(b-h+μd)),n0=(b-h+μd)(σ+h+μ)>0.
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
(16)b-h-μ>0.
Finally, we obtain another solution 0<Ne<K, Se,Ee,Ie>0, where
(17)Ne=12n2(-n1+n12+4n2n0).
This solution is positive under the condition b-h+μd>0.


## 3.2. 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
(18)J(K,0,0)=[-b+μ+hh+μ-b+μ+h-βμα(b-μ-h)0-b-σβμα(b-μ-h)0σ-b+h-μd]
has one negative eigenvalue λ=-b+μ+h (see (16)) and the two eigenvalues of the matrix
(19)J1=[-b-σβμα(b-μ-h)σ-b+h-μd].
Since all parameters are positive, the trace of J1 is negative, due to (16). If the determinant of J1 is positive, then the eigenvalues of J1 are negative. Equivalently, we require
(20)R0≡βσ(b-μ-h)μα(σ+b)(b-h+μd)<1.Therefore, we have proved the following proposition.Proposition 1.
The equilibrium(K,0,0) is locally asymptotically stable for R0<1 and unstable for R0>1.Note that (20) is equivalent to
(21)K<(σ+b)(b-h+μd)βσ≡KC,
where KC is the critical carrying capacity. As a result, we have the following corollary.Corollary 2.
The equilibrium(K,0,0) is stable for K<KC and unstable for K>KC.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 Section6.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β.SolvingNe=K, we obtain the equation (22)1μα(b-h-μ)=12βCμα{[(μα(b-h+μd)-βC(σ-b+h+μ)]2+4βCμα(b-h+μd)(σ+h+μ)μα(b-h+μd)-βC(σ-b+h+μ)+[(μα(b-h+μd)-βC(σ-b+h+μ)]2+4βCμα(b-h+μd)(σ+h+μ)}.This implies
(23)σ(b-h-μ)=μαβC(b-h+μd)(b+σ),
from which we conclude
(24)β=μα(b-h+μd)(b+σ)σ(b-h-μ)≡βC.
Thus β=βC if and only if Ne=K. Also, it is clear that β=βC if and only if R0=1. Note also that, at this critical rate of transmission, we have Ee=Ie=0 and Se=K; see (13).As a result, we have the following two corollaries.Corollary 3.
Suppose0<Ne≤K. Then the following are equivalent: (i) β=βC, (ii) Ne=K, (iii) K=KC, (iv) R0=1 and (v) (Se,Ee,Ie)=(K,0,0).Corollary 4.
Suppose0<Ne≤K holds. Then(1)
forβ<βC, we have  R0<1 so the disease-free equilibrium  (K,0,0) is locally stable;(2)
forβ>βC, we have  R0>1 so the disease-free equilibrium  (K,0,0) is unstable.Remarks.
Suppose (16) holds.(1)
The EE exists if and only ifR0>1.(2)
Note thatβ>βC implies Ne<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.
ForR0>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 (u1,u2,u3) and left eigenvector (v1,v2,v3), where
(25)u1=(1-b-h-μdσ-b-h-μdb-μ-h)u3,u2=(b-h+μd)σu3
and u3>0 is free. Also,
(26)v1=0,v2=σb+σv3
and v3>0 is free. We have
(27)∑k,i=1n‍vkui∂2fk∂xi∂β(E0,βC)=Kv2u3>0,where E0 is the DFE. This implies that (i) or (iv) in Theorem 4 [24] is applicable. However, for β<βC, E0 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, E0 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.


## 4. 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(28)S′(t)=b(S+E+ηA)-μ(1+αN)S-hS-βAAS-βIIS,(29)E′(t)=b(1-η)A-μ(1+αN)E+βAAS+βIIS-(σE+h)E,(30)A′(t)=σEE-μ(1+αN)A-(σA+h)A,(31)I′(t)=-μ(1+αN)I+σAA-μdI.


### 4.1. The Equilibria

The equilibria of (28)–(31) satisfy
(32)Ie=(b-h-μ(1+αNe))Neb-h+μdAe=(μd+μ(1+αNe))σA(b-h+μd)×(b-h-μ(1+αNe))NeσA(b-h+μd)Ee=(σA+h+μ(1+αNe))σEσA(b-h+μd)×(μd+μ(1+αNe))(b-h-μ(1+αNe))NeσEσA(b-h+μd)×(b-h-μ(1+αNe))NeσEσA(b-h+μd)Se=(μd+μ(1+αNe))NeσEσA(b-h+μd)×(σEσA+(b-h-μ(1+αNe))×(-σE-σA-h-μ(1+αNe))).
Equation (29) with E′=0 implies
(33)Ne(Ne-K)(n4Ne4+n3Ne3+n2Ne2+n1Ne+n0)=0,
with
(34)n4=βA(μα)3,(35)n3=(μα)2(βAμ+μd+βIσA)+βA(μα)2(-b+2h+2μ+σE+σA),(36)n2=βAμα(σEσA+(b-h-μ)(-σE-σA-h-μ))+μα(-b+2h+2μ+σE+σA)(βAμ+μd+βI)-σA(b-h+μd)(μα)2,(37)n1=(σEσA+(b-h-μ)(-σE-σA-h-μ))×(βAμ+μd+βI)-σAμα(b-h+μd)(2μ+2h+σE+σA),(38)n0=σA(b-h+μd)×(b(1-η)σE-(μ+h+σE)(μ+h+σA)).
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
(39)n4Ne4+n3Ne3+n2Ne2+n1Ne+n0=0,
we have 0<Ne<K, Se,Ee,Ae,Ie>0.We define the polynomial(40)P(Ne)≡n4Ne4+n3Ne3+n2Ne2+n1Ne+n0,
which satisfies(1)
P(0)<0 for η close to 1, and(2)
P(Ne)→∞ as Ne→∞.Therefore, (39) has at least one positive root for η close to 1 (i.e., when there is not too much vertical transmission).


### 4.2. Stability

The zero equilibrium is unstable when condition (16) holds; that is, b>h+μ. For the DFE (K,0,0,0), the Jacobian matrix(41)J(K,0,0,0)=[-b+μ+hh+μh+μ-βAμα(b-μ-h)-b(1-η)+h+μ-βIμα(b-h-μ)0-b-σEb(1-η)+βAμα(b-μ-h)βIμα(b-h-μ)0σE-b-σA000σA-b+h-μd]has the negative eigenvalue λ=-b+μ+h (see (16)) and three eigenvalues of the matrix
(42)J1=[-b-σEb(1-η)+βAμα(b-μ-h)βIμα(b-h-μ)σE-b-σA00σA-b+h-μd],
which are the zeros of the characteristic polynomial λ3+a2λ2+a1λ+a0=0, with coefficients
(43)a2=σE+3b+σA-h+μd>0,(44)a1(σE+b)(2b+σA-h+μd)-σEβEμα(b-h-μ),a0=(σE+b)(σA+b)(b-h+μd)-σE(b-h+μd)×(b(1-η)+βAμα(b-h-μ))-σEσAβIμα(b-h-μ).
Note that if a0=0, then a1>0, which means that there are two negative eigenvalues and one zero eigenvalue. Also, a0=0 implies
(45)R0≡(b-μ-h)σE(βA(b-h+μd)+βIσA)μα(b-h+μd)(σE+b)(σA+b)+μασE(1-η)b(b-h+μd)μα(b-h+μd)(σE+b)(σA+b)=1.
In order to have three negative eigenvalues, we require a0>0, which is equivalent to R0<1.Therefore, we have proved the following proposition.Proposition 6.
The equilibrium(K,0,0,0) is stable for R0<1 and unstable for R0>1.Condition (45) is equivalent to
(46)K<(b-h+μd)((σE+b)(σA+b)-σE(1-η)b)σE(βA(b-h+μd)+βIσA)≡KC,
the critical carrying capacity. As a result, we have the following corollary.Corollary 7.
The equilibrium(K,0,0,0) is stable for K<KC and unstable for K>KC.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.
SubstitutingNe=K in (29) with E′=0 and solving for βA(b-h+μd)+βIσA, we obtain
(47)βA(b-h+μd)+βIσA=μα(b-h+μd)((σE+b)(σA+b)-σE(1-η)b)σE(b-h-μ)≡βC.
We have thus proved that βA(b-h+μd)+βIσA=βC if and only if Ne=K. Also, βA(b-h+μd)+βIσA=βC if and only if R0=1. Moreover, at this critical rate of transmission, we have Ee=Ae=Ie=0 and Se=K.As a result, we have the following two corollaries.Corollary 9.
Suppose0<Ne≤K. Then the following are equivalent: (i) βA(b-h+μd)+βIσA=βC, (ii) Ne=K, (iii) K=KC, (iv) R0=1, and (v) (Se,Ee,Ae,Ie)=(K,0,0,0).Corollary 10.
Suppose0<Ne≤K. Then(1)
forβA(b-h+μd)+βIσA<βC, we have  R0<1; so the disease-free equilibrium  (K,0,0,0) is locally asymptotically stable;(2)
forβA(b-h+μd)+βIσA>βC, we have  R0>1; so the disease-free equilibrium  (K,0,0,0) is unstable.Remarks.
Suppose (16) holds.(1)
ForR0>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 R0<1, the existence of the EE is not clear.(2)
The conditionβ>βC implies Ne<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 thatR0 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.


## 4.1. The Equilibria

The equilibria of (28)–(31) satisfy
(32)Ie=(b-h-μ(1+αNe))Neb-h+μdAe=(μd+μ(1+αNe))σA(b-h+μd)×(b-h-μ(1+αNe))NeσA(b-h+μd)Ee=(σA+h+μ(1+αNe))σEσA(b-h+μd)×(μd+μ(1+αNe))(b-h-μ(1+αNe))NeσEσA(b-h+μd)×(b-h-μ(1+αNe))NeσEσA(b-h+μd)Se=(μd+μ(1+αNe))NeσEσA(b-h+μd)×(σEσA+(b-h-μ(1+αNe))×(-σE-σA-h-μ(1+αNe))).
Equation (29) with E′=0 implies
(33)Ne(Ne-K)(n4Ne4+n3Ne3+n2Ne2+n1Ne+n0)=0,
with
(34)n4=βA(μα)3,(35)n3=(μα)2(βAμ+μd+βIσA)+βA(μα)2(-b+2h+2μ+σE+σA),(36)n2=βAμα(σEσA+(b-h-μ)(-σE-σA-h-μ))+μα(-b+2h+2μ+σE+σA)(βAμ+μd+βI)-σA(b-h+μd)(μα)2,(37)n1=(σEσA+(b-h-μ)(-σE-σA-h-μ))×(βAμ+μd+βI)-σAμα(b-h+μd)(2μ+2h+σE+σA),(38)n0=σA(b-h+μd)×(b(1-η)σE-(μ+h+σE)(μ+h+σA)).
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
(39)n4Ne4+n3Ne3+n2Ne2+n1Ne+n0=0,
we have 0<Ne<K, Se,Ee,Ae,Ie>0.We define the polynomial(40)P(Ne)≡n4Ne4+n3Ne3+n2Ne2+n1Ne+n0,
which satisfies(1)
P(0)<0 for η close to 1, and(2)
P(Ne)→∞ as Ne→∞.Therefore, (39) has at least one positive root for η close to 1 (i.e., when there is not too much vertical transmission).


## 4.2. Stability

The zero equilibrium is unstable when condition (16) holds; that is, b>h+μ. For the DFE (K,0,0,0), the Jacobian matrix(41)J(K,0,0,0)=[-b+μ+hh+μh+μ-βAμα(b-μ-h)-b(1-η)+h+μ-βIμα(b-h-μ)0-b-σEb(1-η)+βAμα(b-μ-h)βIμα(b-h-μ)0σE-b-σA000σA-b+h-μd]has the negative eigenvalue λ=-b+μ+h (see (16)) and three eigenvalues of the matrix
(42)J1=[-b-σEb(1-η)+βAμα(b-μ-h)βIμα(b-h-μ)σE-b-σA00σA-b+h-μd],
which are the zeros of the characteristic polynomial λ3+a2λ2+a1λ+a0=0, with coefficients
(43)a2=σE+3b+σA-h+μd>0,(44)a1(σE+b)(2b+σA-h+μd)-σEβEμα(b-h-μ),a0=(σE+b)(σA+b)(b-h+μd)-σE(b-h+μd)×(b(1-η)+βAμα(b-h-μ))-σEσAβIμα(b-h-μ).
Note that if a0=0, then a1>0, which means that there are two negative eigenvalues and one zero eigenvalue. Also, a0=0 implies
(45)R0≡(b-μ-h)σE(βA(b-h+μd)+βIσA)μα(b-h+μd)(σE+b)(σA+b)+μασE(1-η)b(b-h+μd)μα(b-h+μd)(σE+b)(σA+b)=1.
In order to have three negative eigenvalues, we require a0>0, which is equivalent to R0<1.Therefore, we have proved the following proposition.Proposition 6.
The equilibrium(K,0,0,0) is stable for R0<1 and unstable for R0>1.Condition (45) is equivalent to
(46)K<(b-h+μd)((σE+b)(σA+b)-σE(1-η)b)σE(βA(b-h+μd)+βIσA)≡KC,
the critical carrying capacity. As a result, we have the following corollary.Corollary 7.
The equilibrium(K,0,0,0) is stable for K<KC and unstable for K>KC.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.
SubstitutingNe=K in (29) with E′=0 and solving for βA(b-h+μd)+βIσA, we obtain
(47)βA(b-h+μd)+βIσA=μα(b-h+μd)((σE+b)(σA+b)-σE(1-η)b)σE(b-h-μ)≡βC.
We have thus proved that βA(b-h+μd)+βIσA=βC if and only if Ne=K. Also, βA(b-h+μd)+βIσA=βC if and only if R0=1. Moreover, at this critical rate of transmission, we have Ee=Ae=Ie=0 and Se=K.As a result, we have the following two corollaries.Corollary 9.
Suppose0<Ne≤K. Then the following are equivalent: (i) βA(b-h+μd)+βIσA=βC, (ii) Ne=K, (iii) K=KC, (iv) R0=1, and (v) (Se,Ee,Ae,Ie)=(K,0,0,0).Corollary 10.
Suppose0<Ne≤K. Then(1)
forβA(b-h+μd)+βIσA<βC, we have  R0<1; so the disease-free equilibrium  (K,0,0,0) is locally asymptotically stable;(2)
forβA(b-h+μd)+βIσA>βC, we have  R0>1; so the disease-free equilibrium  (K,0,0,0) is unstable.Remarks.
Suppose (16) holds.(1)
ForR0>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 R0<1, the existence of the EE is not clear.(2)
The conditionβ>βC implies Ne<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 thatR0 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.


## 5. The Full Model

We now analyze stability in the full, gender-structured model.LetEf=Em=If=Im=Af=Am=0. Then (6)-(7) imply
(48)Sf′=bSfSmkSf+Sm-μf(1+αfN)Sf-hfSf,Sm′=bSfSmkSf+Sm-μm(1+αmN)Sm-hmSm.
The equilibria thus satisfy
(49)bSfSmkSf+Sm-μf(1+αN)Sf-hfSf=0,(50)bSfSmkSf+Sm-μm(1+αN)Sm-hmSm=0.
We have two cases.Case 1.
IfSf=Sm=0, we have the zero equilibrium.Case 2.
The nonzero disease-free equilibrium is calculated as follows. First, we multiply (49) by μm/Sf and (50) by μf/Sm. Then we solve the resulting system for z≡Sm/Sf and K=Sf+Sm. We obtain
(51)z=bμf-k(μfhm-μmhf)bμm+(μfhm-μmhf),K=1αμf(bzk+z-hf-μf)=1αμm(bk+z-hm-μm).
Equality of the two expressions for K results from the solution for z. For the equilibrium densities for male and females, we obtain
(52)Sf=K1+z,Sm=zK1+z.
A positive equilibrium requires K>0, which implies
(53)bk+z-hm-μm>0,or,equivalently,bzk+z-hf-μf>0.
To investigate the stability of the disease-free equilibrium, we use the method of the next-generation operator [25, 26]. Set
(54)r1=σE+bzk+z,s1=σE+bk+z,r2=σA+bzk+z,s2=σA+bk+z,r3=bzk+z-hf+μd,s3=bk+z-hm+μd.
The Jacobian matrix for the equilibrium(Ef,Em,If,Im,Af,Am,Sf,Sm) = (0,0,0,0,0,0,K/(1+z),zK/(1+z)) has the form(55)M=[F-VOJ1J2],
where(56)F=[00z(1-η)bk+z+βAfγfK1+zβAfγmfK1+zβIfγfK1+zβIfγmfK1+z00z(1-η)bk+z+βAmγfmzK1+zβAmγmzK1+zβImγfmzK1+zβImγmzK1+z000000000000000000000000],V=[r1000000s10000-σE0r20000-σE0s20000-σA0r30000-σA0s3],O=[000000000000],J2=[bz2(k+z)2-bzk+z-11+zr3bk(k+z)2-11+zr3bz2(k+z)2-z1+zs3bk(k+z)2-bk+z-z1+zs3].The 2×6 matrix J1 turns out to be irrelevant to our calculations. We have
(57)V-1=[1r10000001s10000σEr1r201r20000σEs1s201s200σEσAr1r2r30σAr2r301r300σEσAs1s2s30σAs2s301s3].
The threshold parameter R0 is the spectral radius of the matrix FV-1, provided that all eigenvalues of J2 have negative real part [26]. Thus
(58)R0=ρ(FV-1)=max{|12b1+12b12-4b2|,|12b1-12b12-4b2|},
where
(59)b1=(z(1-η)bσE/(k+z))+σEβAfγf(K/(1+z))r1r2+σEσAβIfγf(K/(1+z))r1r2r3+σEβAmγm(Kz/(1+z))s1s2+σEσAβImγm(Kz/(1+z))s1s2s3,(60)b2=((zσE(1-η)b/(k+z))+σEβAfγf(K/(1+z))r1r2+σEσAβIfγf(K/(1+z))r1r2r3)×(σEβAmγm(Kz/(1+z))s1s2+σEσAβimγm(Kz/(1+z))s1s2s3)-(σE((z(1-η)b/(k+z))+βAmγfm(zK/(1+z)))r1r2+σEσAβImγfm(zK/(1+z))r1r2r3)×(σEβAfγmf(K/(1+z))s1s2+σEσAβIfγmf(K/(1+z))s1s2s3).
Note that if b12-4b2≥0, then R0=(1/2)b1+(1/2)b12-4b2. In order to prove that J2 has only eigenvalues with negative real part, we prove that tr(J2)<0 and det(J2)>0. First, we note the inequalities
(61)bz2(k+z)2-bzk+z=-bkz(k+z)2≤0,bk(k+z)2-bk+z=-bz(k+z)2≤0.
From (53), we obtain tr(J2)<0.
Now, using (53), (61), as well as the fact that z and all parameters are positive, we have
(62)det(J2)=(-bkz(k+z)2+11+z(bzk+z-hf-μf))×(-bzk+z+z1+z(bk+z-hm-μm))-(bk(k+z)2-11+z(bzk+z-hf-μf))×(bz2(k+z)2-z1+z(bk+z-hm-μm))=bkz(k+z)2z1+z(bk+z-hm-μm)+11+z(bzk+z-hf-μf)bzk+z+bk(k+z)2z1+z(bk+z-hm-μm)+11+z(bzk+z-hf-μf)bz2(k+z)2>0.
As a result, we have proved the following theorem.Theorem 11.
The DFE is stable forR0<1 and unstable for R0>1.Note that the value ofR0 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βAs→0, β=βIs for s∈{f,m},η→1, σ=σE and σA→∞. Then we have b1→a1 and b2→a2, where
(63)a1=σK1+z(βfγf(σ+(bz/(k+z)))r3+zβmγm(σ+(b/(k+z)))s3),a2=σK1+z·zβfβm(γfγm-γfmγmf)(σ+(bz/(k+z)))(σ+(b/(k+z)))r3s3.
Moreover, if γf=γfm and γm=γmf, then R0=1/2a1. If k=0 and all parameters are identical for males and females, then z=1 and
(64)R0=βγσ(b-h-μ)μα(σ+b)(b-h+μd),
which is the value of R0 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 R0 for the full model (58) matches R0 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.


## 6. Application to CWD

For the purpose of simulations, we used Wildlife Management Unit 150-Alberta (WMU 150), which has an area equal to 1841 km2 [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.Table 2
Sample values and sensitivity ranges for mule deer.
Parameter
Range
Sample value
Units
Reference
b
0.5–0.9
0.57
Years−1
[28]
β
0–0.05
0.0326
Individuals−1 Years−1
[6]
h
0-1
0.1512
Years−1
[29]
h
f
0.0922
Years−1
[29]
h
m
0.282
Years−1
[29]
L
1-2
1.5
Years
[30]
L
E
1.5
Years
[30]
L
A
0.5
Years
[28]
μ
0.045–0.15
0.1053
Years−1
[28]
μ
f
0.1
Years−1
[28]
μ
m
0.1111
Years−1
[28]
α
0.002–0.008
0.004
Individuals−1
Assumed
μ
d
1–4
2
Years−1
[6]
σ
0.045–0.15
0.05
Years−1
[28]
k
0.2
—
[31]
γ
f
=
γ
m
=
γ
m
f
=
γ
f
m
0.002
—
AssumedWe 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 hf=0.0922, hm=0.282, and h=0.1512 (estimated using P(t)=P0e-ht, where P(t) is the population after t years and P0 is the initial population).(2)
The average lifespan for bucks and does is 9 and 10 years, respectively [28], so the background death rates are estimated by μf=0.1, μm=0.1111, and μ=0.1053 (μ=1/L, where L is the lifespan).(3)
Following a widespread culling program, WMU 150 had650 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 E0=24. We assume that there are no infected animals with clinical signs (If0=Im0=0). All other animals are considered susceptible (Sf0+Sm0=626). Also, we assume that the disease prevalence in males is initially twice that in females [32], so we have Ef0=8 and Em0=16.(4)
It is extremely difficult to obtain reasonable estimates or even guesses 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 βIf=βAf=βf, βIm=βAm=βm, and βf=β-0.01, βm=β+0.01. See Table 2.For the basic model, the endemic equilibrium(Se,Ee,Ie)=(100,56,17) is stable and the basic reproduction number is R0=5.407828. For the basic model with asymptomatic individuals, the endemic equilibrium (Se,Ee,Ae,Ie)=(47,27,12,11) is stable and the basic reproduction number is R0=7.026539. Finally, for the two-gender nonasymptomatic model, with η=1, σA→∞ and βA→0, the endemic equilibrium (Sf,Sm,Ef,Em,If,Im)=(103,73,15,8,4,2) is stable and the basic reproduction number is R0=1.56363. These numbers indicate that the disease persists in the absence of control measures. Note how prevalence and R0 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. Figure2 shows the graph of R0 and K as functions of harvesting for the basic model. Any harvesting rate between h1=0.4130 and h2=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 2
Critical harvesting rates in the basic model. Too little harvesting results in disease persistence, while too much eradicates the population.Deer population in the WMU 150 area for the basic model when (a)h<h1, resulting in disease persistence; (b) h1<h<h2, resulting in disease control; (c) h>h2, resulting in eradication of the population (solid curve; susceptible; diamonds: exposed; dashed curve: infected).
(a)
(b)
(c)Figure4 shows the graph of R0 and K as functions of harvesting for the basic model with asymptomatic individuals. Any harvesting rate between h1=0.4262 and h2=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 4
Critical harvesting rates in the basic model with asymptomatic individuals. Too little harvesting results in disease persistence. Too much eradicates the population. Note that the optimum harvesting interval is tighter when asymptomatic individuals are included.Deer population in the WMU 150 area for the basic model with asymptomatic individuals when (a)h<h1, resulting in disease persistence; (b) h1<h<h2, resulting in disease control; (c) h>h2, resulting in eradication of the population. (solid curve; susceptible; diamonds: exposed; squares: asymptomatic; dashed curve: infected).
(a)
(b)
(c)Figure6 shows the level curves R0=1 and K=0 as functions of gender-specific harvesting parameters hf and hm for the full model. Any point (hf,hm) between R0=1 and K=0 implies that the disease is controlled and the population survives. In Figure 7, we take three different points ((hf,hm)=(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.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.Deer population in the WMU 150 area for the full model. ((a) and (b)): Ifh<h1, then the disease persists, although it is stronger in males. ((c) and (d)): If h1<h<h2, 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>h2, 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).
(a)
(b)
(c)
(d)
(e)
(f)


### 6.1. Sensitivity Analysis

Due to the degree of uncertainty in the parameters values, we considered a range of parameters to examine the dependence ofR0 on parameter variation for the basic model. We used Latin hypercube sampling and partial rank correlation coefficients (PRCCs) to identify which parameters R0 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.Figure8 plots PRCCs for each input parameter. This demonstrates that R0 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.Figure 8
Partial rank correlation coefficients indicate that harvesting has the greatest effect onR0. Parameters with PRCC >0 will increase R0 when they are increased, while parameters with PRCC <0 will decrease R0 when they are increased.Monte Carlo simulations for 1000 runs drawn from parameter ranges using Latin hypercube sampling. Ifβ is sufficiently small or h is sufficiently large, then the disease can be controlled. Note that h has been constrained to ensure population survival. However, if h is in the approximate range 0.6<h<0.8, then the disease is controlled and the population survives.
(a)
(b)
(c)
(d)
(e)
(f)
(g)It should be noted that variations inh will change R0 from values much greater than 1 to very low values, resulting in significant dependence of R0 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.


## 6.1. Sensitivity Analysis

Due to the degree of uncertainty in the parameters values, we considered a range of parameters to examine the dependence ofR0 on parameter variation for the basic model. We used Latin hypercube sampling and partial rank correlation coefficients (PRCCs) to identify which parameters R0 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.Figure8 plots PRCCs for each input parameter. This demonstrates that R0 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.Figure 8
Partial rank correlation coefficients indicate that harvesting has the greatest effect onR0. Parameters with PRCC >0 will increase R0 when they are increased, while parameters with PRCC <0 will decrease R0 when they are increased.Monte Carlo simulations for 1000 runs drawn from parameter ranges using Latin hypercube sampling. Ifβ is sufficiently small or h is sufficiently large, then the disease can be controlled. Note that h has been constrained to ensure population survival. However, if h is in the approximate range 0.6<h<0.8, then the disease is controlled and the population survives.
(a)
(b)
(c)
(d)
(e)
(f)
(g)It should be noted that variations inh will change R0 from values much greater than 1 to very low values, resulting in significant dependence of R0 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.


## 7. 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 disease-free 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 density-dependent 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
(65)K=b-h-μκ(μ+h),
and the basic reproduction number
(66)R0=βσ(b-μ-h)κ(μ+h)(μ+σ+h)(μ+μd).
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.Figure 10
Basic reproduction number and critical harvesting rates for a variant of the model that considers density-dependent births and density-independent deaths. The disease persists whenR0>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.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 frequency-dependent 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 ofR0 calculated for the models without gender structure is relatively high (in fact, most likely too high), the value for R0 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 gender-structured 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].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 reduceR0 and increase the critical transmission rate. This reduction in R0 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 thatR0 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 whenR0 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 R0=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.



---

*Source: 802450-2012-09-20.xml*