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

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.

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 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.

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 𝑁: 𝑁^′=𝐵(𝑁)𝑁−𝐷(𝑁)𝑁−ℎ𝑁, (1) where 𝐵 and 𝐷 are the per capita birth and death rates, respectively, and ℎ 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 𝐵(𝑁)−𝐷(𝑁) is a monotone decreasing function, (1) will exhibit logistic-like behaviour.

Next, we split the population into compartments according to disease state as susceptible (𝑆), exposed (𝐸), and infectious (𝐼). We assume that all individuals are born susceptible; that is, there is no vertical transmission. Infectious individuals experience an additional disease-induced mortality 𝜇_𝑑. 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) where 𝑁=𝑆+𝐸+𝐼. If infectious individuals are unable to reproduce, then the birth term in the 𝑆-equation reads 𝐵(𝑁)(𝑆+𝐸). 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 𝜇_𝑑 as the combined effect of disease-induced death and disease-specific culling. When the disease is not easily detected, then we can replace 𝜇_𝑑 with 𝜇_𝑑+ℎ 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 𝜇_𝑑+ℎ≈𝜇_𝑑. 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 (𝐴). We allow for differential transmission rates from the asymptomatic and infectious class (𝛽_𝐴,𝛽_𝐼) and denote the sojourn times in the exposed and asymptomatic class as 1/𝜎_𝐸 and 1/𝜎_𝐴, 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 𝑆^′=𝐵(𝑁)𝑆+𝐸+𝜂𝐴−𝐷(𝑁)𝑆−ℎ𝑆−𝛽_𝐴𝐴+𝛽_𝐼𝐼𝐸𝑆,_′=𝐵(𝑁)1−𝜂𝐴−𝐷(𝑁)𝐸+𝛽_𝐴𝐴+𝛽_𝐼𝐼𝑆−𝜎_𝐸𝐴+ℎ𝐸,_′=𝜎_𝐸𝐸−𝐷(𝑁)𝐴−𝜎_𝐴𝐼+ℎ𝐴,_′=−𝐷(𝑁)𝐼+𝜎_𝐴𝐴−𝜇_𝑑𝐼. (3)

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 classes 𝑁=𝑁_𝑓+𝑁_𝑚, where 𝑁_𝑓 is the female population and 𝑁_𝑚 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 𝑁_𝑓 and 𝑁_𝑚 as 𝑁^′_𝑓(𝑡)=𝐵𝑁_𝑓,𝑁_𝑚𝑁_𝑓−𝐷(𝑁)𝑁_𝑓−ℎ_𝑓𝑁_𝑓,𝑁^′_𝑚(𝑡)=𝐵𝑁_𝑓,𝑁_𝑚𝑁_𝑚−𝐷(𝑁)𝑁_𝑚−ℎ_𝑚𝑁_𝑚, (4) 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)

In this case, 𝑏 is the fecundity, which is half the maximum per-female number of offspring at low densities, and 𝑘⁻¹ 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 𝛾_𝑞(𝑞∈{𝑚,𝑓,𝑓𝑚,𝑚𝑓}) and transmissibility 𝛽^𝑖_𝑗(𝑗∈{𝐴,𝐼}) can depend on gender. For simplicity, we assume that infectious individuals cannot reproduce and denote the effective population size for reproduction as 𝑁_𝑓=𝑆_𝑓+𝐸_𝑓+𝐴_𝑓 and 𝑁_𝑚, accordingly. This assumption is justified for CWD. We have the following system of equations: 𝑆^′_𝑓(𝑡)=𝐵𝑁_𝑓,𝑁_𝑚𝑆_𝑓+𝐸_𝑓+𝜂𝐴_𝑓−𝐷(𝑁)𝑆_𝑓−ℎ_𝑓𝑆_𝑓−𝐹_𝑓𝑆_𝑓,𝑆^′_𝑚(𝑡)=𝐵𝑁_𝑓,𝑁_𝑚𝑆_𝑓+𝐸_𝑓+𝜂𝐴_𝑓−𝐷(𝑁)𝑆_𝑚−ℎ_𝑚𝑆_𝑚−𝐹_𝑚𝑆_𝑚,𝐸^′_𝑓(𝑡)=𝐵𝑁_𝑓,𝑁_𝑚1−𝜂𝐴_𝑓−𝜎_𝐸𝐸_𝑓−𝐷(𝑁)𝐸_𝑓−ℎ_𝑓𝐸_𝑓+𝐹_𝑓𝑆_𝑓,𝐸^′_𝑚(𝑡)=𝐵𝑁_𝑓,𝑁_𝑚1−𝜂𝐴_𝑓−𝜎_𝐸𝐸_𝑚−𝐷(𝑁)𝐸_𝑚−ℎ_𝑚𝐸_𝑚+𝐹_𝑚𝑆_𝑚,𝐴^′_𝑓(𝑡)=𝜎_𝐸𝐸_𝑓−𝜎_𝐴𝐴_𝑓−𝐷(𝑁)𝐴_𝑓−ℎ_𝑓𝐴_𝑓,𝐴^′_𝑚(𝑡)=𝜎_𝐸𝐸_𝑚−𝜎_𝐴𝐴_𝑚−𝐷(𝑁)𝐴_𝑚−ℎ_𝑚𝐴_𝑚,𝐼^′_𝑓(𝑡)=𝜎_𝐴𝐴_𝑓−𝐷(𝑁)𝐼_𝑓−𝜇_𝑑𝐼_𝑓,𝐼^′_𝑚(𝑡)=𝜎_𝐴𝐴_𝑚−𝐷(𝑁)𝐼_𝑚−𝜇_𝑑𝐼_𝑚. (6) Here, the force of infection is given by 𝐹_𝑓=𝛽^𝑓_𝐴𝛾_𝑓𝐴_𝑓+𝛾_𝑚𝑓𝐴_𝑚+𝛽^𝑓_𝐼𝛾_𝑓𝐼_𝑓+𝛾_𝑚𝑓𝐼_𝑚,𝐹_𝑚=𝛽^𝑚_𝐴𝛾_𝑓𝑚𝐴_𝑓+𝛾_𝑚𝐴_𝑚+𝛽^𝑚_𝐼𝛾_𝑓𝑚𝐼_𝑓+𝛾_𝑚𝐼_𝑚. (7) The model is illustrated in Figure 1, with variables and parameters listed in Table 1.

Table 1 

Parameter definitions.

Figure 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 𝑅₀, also known as the basic reproduction number [21, 22]. The model reads 𝑆^′(𝑡)()()𝐸=𝑏𝑆+𝐸−𝜇1+𝛼𝑁𝑆−ℎ𝑆−𝛽𝐼𝑆,_′(𝑡)()()𝐼=−𝜇1+𝛼𝑁𝐸+𝛽𝐼𝑆−𝜎+ℎ𝐸,_′(𝑡)=−𝜇(1+𝛼𝑁)𝐼+𝜎𝐸−𝜇_𝑑𝐼, (8)(9)(10) 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 𝑁. Adding the steady state equations for 𝑆 and 𝐸, we obtain (𝑏−ℎ)𝑆_𝑒+𝐸_𝑒−𝜇1+𝛼𝑁_𝑒𝑆_𝑒+𝐸_𝑒−𝜎𝐸_𝑒=0. (11) Here, we use 𝑆_𝑒, 𝐸_𝑒, 𝐼_𝑒, and 𝑁_𝑒 to represent the equilibrium values.

Substituting 𝑆+𝐸=𝑁−𝐼 and solving for 𝐸 in terms of 𝐼 and 𝑁 from (10) with 𝐼^′=0, we obtain (𝑏−ℎ)𝑁_𝑒−𝜇1+𝛼𝑁_𝑒𝑁_𝑒−𝑏−ℎ+𝜇_𝑑𝐼_𝑒=0. (12) We can express the steady-state values of 𝑆_𝑒, 𝐸_𝑒,  and  𝐼_𝑒 in terms of 𝑁_𝑒 as 𝐼_𝑒=𝑏−ℎ−𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝑏−ℎ+𝜇_𝑑,𝐸_𝑒=𝑏−ℎ−𝜇1+𝛼𝑁_𝑒𝜇_𝑑+𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝜎𝑏−ℎ+𝜇_𝑑,𝑆_𝑒=𝑁_𝑒−𝐼_𝑒−𝐸_𝑒=𝑁_𝑒𝜇_𝑑+𝜇1+𝛼𝑁_𝑒𝜎−𝑏+ℎ+𝜇1+𝛼𝑁_𝑒𝜎𝑏−ℎ+𝜇_𝑑, (13) under the condition 𝑏−ℎ+𝜇_𝑑≠0. Using (9) with 𝐸^′=0, we obtain 𝑁_𝑒𝑁_𝑒−𝐾𝑛₂𝑁²_𝑒+𝑛₁𝑁_𝑒−𝑛₀=0, (14) where 𝐾=(1/𝜇𝛼)(𝑏−ℎ−𝜇) and 𝑛₂𝑛=𝛽𝜇𝛼,₁=𝜎−𝑏+ℎ+𝜇𝛽−𝜇𝛼𝑏−ℎ+𝜇_𝑑,𝑛₀=𝑏−ℎ+𝜇_𝑑𝜎+ℎ+𝜇>0. (15) There are solutions 𝑁=0, corresponding to the absence of population, and 𝑆=𝑁=𝐾, 𝐸=𝐼=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 (𝐾≠0), we require 𝑏−ℎ−𝜇>0. (16) Finally, we obtain another solution 0<𝑁_𝑒<𝐾, 𝑆_𝑒,𝐸_𝑒,𝐼_𝑒>0, where 𝑁_𝑒=12𝑛₂−𝑛₁+𝑛²₁+4𝑛₂𝑛₀. (17) This solution is positive under the condition 𝑏−ℎ+𝜇_𝑑>0.

3.2. Stability

The zero equilibrium of system (8)–(10) is unstable if condition (16) holds; that is, 𝑏>ℎ+𝜇. We start by analyzing the stability of the DFE (𝐾,0,0). The corresponding Jacobian matrix 𝐽()=𝛽𝐾,0,0−𝑏+𝜇+ℎℎ+𝜇−𝑏+𝜇+ℎ−𝛽𝜇𝛼𝑏−𝜇−ℎ0−𝑏−𝜎𝜇𝛼𝑏−𝜇−ℎ0𝜎−𝑏+ℎ−𝜇_𝑑 (18) has one negative eigenvalue 𝜆=−𝑏+𝜇+ℎ (see (16)) and the two eigenvalues of the matrix 𝐽₁=𝛽−𝑏−𝜎𝜇𝛼𝑏−𝜇−ℎ𝜎−𝑏+ℎ−𝜇_𝑑. (19) Since all parameters are positive, the trace of 𝐽₁ is negative, due to (16). If the determinant of 𝐽₁ is positive, then the eigenvalues of 𝐽₁ are negative. Equivalently, we require 𝑅₀≡𝛽𝜎𝑏−𝜇−ℎ𝜇𝛼(𝜎+𝑏)𝑏−ℎ+𝜇_𝑑<1. (20)

Therefore, we have proved the following proposition.

Proposition 1. The equilibrium (𝐾,0,0) is locally asymptotically stable for 𝑅₀<1 and unstable for 𝑅₀>1.

Note that (20) is equivalent to ()𝐾<𝜎+𝑏𝑏−ℎ+𝜇_𝑑𝛽𝜎≡𝐾_𝐶, (21) where 𝐾_𝐶 is the critical carrying capacity. As a result, we have the following corollary.

Corollary 2. The equilibrium (𝐾,0,0) is stable for 𝐾<𝐾_𝐶 and unstable for 𝐾>𝐾_𝐶.

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 𝑁_𝑒=𝐾, we obtain the equation 1=1𝜇𝛼𝑏−ℎ−𝜇2𝛽_𝐶𝜇𝛼𝜇𝛼𝑏−ℎ+𝜇_𝑑−𝛽_𝐶+𝜎−𝑏+ℎ+𝜇(𝜇𝛼(𝑏−ℎ+𝜇_𝑑)−𝛽_𝐶(𝜎−𝑏+ℎ+𝜇)²+4𝛽_𝐶𝜇𝛼𝑏−ℎ+𝜇_𝑑𝜎+ℎ+𝜇. (22)This implies 𝜎𝑏−ℎ−𝜇=𝜇𝛼𝛽_𝐶𝑏−ℎ+𝜇_𝑑(𝑏+𝜎), (23) from which we conclude 𝛽=𝜇𝛼𝑏−ℎ+𝜇_𝑑(𝑏+𝜎)𝜎𝑏−ℎ−𝜇≡𝛽_𝐶. (24) Thus 𝛽=𝛽_𝐶 if and only if 𝑁_𝑒=𝐾. Also, it is clear that 𝛽=𝛽_𝐶 if and only if 𝑅₀=1. Note also that, at this critical rate of transmission, we have 𝐸_𝑒=𝐼_𝑒=0 and 𝑆_𝑒=𝐾; see (13).

As a result, we have the following two corollaries.

Corollary 3. Suppose 0<𝑁_𝑒≤𝐾. Then the following are equivalent: (i) 𝛽=𝛽_𝐶, (ii) 𝑁_𝑒=𝐾, (iii) 𝐾=𝐾_𝐶, (iv) 𝑅₀=1 and (v) (𝑆_𝑒,𝐸_𝑒,𝐼_𝑒)=(𝐾,0,0).

Corollary 4. Suppose 0<𝑁_𝑒≤𝐾 holds. Then(1)for 𝛽<𝛽_𝐶, we have 𝑅₀<1 so the disease-free equilibrium (𝐾,0,0) is locally stable;(2)for 𝛽>𝛽_𝐶, we have 𝑅₀>1 so the disease-free equilibrium (𝐾,0,0) is unstable.

Remarks. Suppose (16) holds.(1)The EE exists if and only if 𝑅₀>1.(2)Note that 𝛽>𝛽_𝐶 implies 𝑁_𝑒<𝐾 and (𝐾,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 𝑅₀>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 𝜆=−𝑏+𝜇+ℎ,−𝑏−𝜎 and 0, when 𝛽=𝛽_𝐶. Also, the zero eigenvalue has right eigenvector (𝑢₁,𝑢₂,𝑢₃) and left eigenvector (𝑣₁,𝑣₂,𝑣₃), where 𝑢₁=1−𝑏−ℎ−𝜇_𝑑𝜎−𝑏−ℎ−𝜇_𝑑𝑢𝑏−𝜇−ℎ³,𝑢₂=𝑏−ℎ+𝜇_𝑑𝜎𝑢₃ (25) and 𝑢₃>0 is free. Also, 𝑣₁𝑣=0,₂=𝜎𝑣𝑏+𝜎³ (26) and 𝑣₃>0 is free. We have ^𝑛_𝑘,𝑖=1𝑣_𝑘𝑢_𝑖𝜕²𝑓_𝑘𝜕_𝑥𝑖𝜕_𝛽𝐸₀,𝛽_𝐶=𝐾𝑣₂𝑢₃>0, (27)where 𝐸₀ is the DFE. This implies that (i) or (iv) in Theorem 4 [24] is applicable. However, for 𝛽<𝛽_𝐶, 𝐸₀ 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 𝛽<𝛽_𝐶 to 𝛽>𝛽_𝐶, 𝐸₀ 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 𝑆^′(𝑡)=𝑏𝑆+𝐸+𝜂𝐴−𝜇(1+𝛼𝑁)𝑆−ℎ𝑆−𝛽_𝐴𝐴𝑆−𝛽_𝐼𝐸𝐼𝑆,_′(𝑡)=𝑏1−𝜂𝐴−𝜇(1+𝛼𝑁)𝐸+𝛽_𝐴𝐴𝑆+𝛽_𝐼𝐼𝑆−𝜎_𝐸𝐴+ℎ𝐸,_′(𝑡)=𝜎_𝐸𝐸−𝜇(1+𝛼𝑁)𝐴−𝜎_𝐴𝐼+ℎ𝐴,_′(𝑡)=−𝜇(1+𝛼𝑁)𝐼+𝜎_𝐴𝐴−𝜇_𝑑𝐼. (28)(29)(30)(31)

4.1. The Equilibria

The equilibria of (28)–(31) satisfy 𝐼_𝑒=𝑏−ℎ−𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝑏−ℎ+𝜇_𝑑𝐴_𝑒=𝜇_𝑑+𝜇1+𝛼𝑁_𝑒𝜎_𝐴𝑏−ℎ+𝜇_𝑑×𝑏−ℎ−𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝜎_𝐴𝑏−ℎ+𝜇_𝑑𝐸_𝑒=𝜎_𝐴+ℎ+𝜇1+𝛼𝑁_𝑒𝜎_𝐸𝜎_𝐴𝑏−ℎ+𝜇_𝑑×𝜇_𝑑+𝜇1+𝛼𝑁_𝑒𝑏−ℎ−𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝜎_𝐸𝜎_𝐴𝑏−ℎ+𝜇_𝑑×𝑏−ℎ−𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝜎_𝐸𝜎_𝐴𝑏−ℎ+𝜇_𝑑𝑆_𝑒=𝜇_𝑑+𝜇1+𝛼𝑁_𝑒𝑁_𝑒𝜎_𝐸𝜎_𝐴𝑏−ℎ+𝜇_𝑑×𝜎_𝐸𝜎_𝐴+𝑏−ℎ−𝜇1+𝛼𝑁_𝑒×−𝜎_𝐸−𝜎_𝐴−ℎ−𝜇1+𝛼𝑁_𝑒. (32) Equation (29) with 𝐸^′=0 implies 𝑁_𝑒𝑁_𝑒−𝐾𝑛₄𝑁⁴_𝑒+𝑛₃𝑁³_𝑒+𝑛₂𝑁²_𝑒+𝑛₁𝑁_𝑒+𝑛₀=0, (33) with 𝑛₄=𝛽_𝐴𝜇𝛼³,𝑛₃=𝜇𝛼²𝛽_𝐴𝜇+𝜇_𝑑+𝛽_𝐼𝜎_𝐴+𝛽_𝐴𝜇𝛼²−𝑏+2ℎ+2𝜇+𝜎_𝐸+𝜎_𝐴,𝑛₂=𝛽_𝐴𝜇𝛼𝜎_𝐸𝜎_𝐴+𝑏−ℎ−𝜇−𝜎_𝐸−𝜎_𝐴−ℎ−𝜇+𝜇𝛼−𝑏+2ℎ+2𝜇+𝜎_𝐸+𝜎_𝐴𝛽_𝐴𝜇+𝜇_𝑑+𝛽_𝐼−𝜎_𝐴𝑏−ℎ+𝜇_𝑑𝜇𝛼²,𝑛₁=𝜎_𝐸𝜎_𝐴+𝑏−ℎ−𝜇−𝜎_𝐸−𝜎_𝐴×𝛽−ℎ−𝜇_𝐴𝜇+𝜇_𝑑+𝛽_𝐼−𝜎_𝐴𝜇𝛼𝑏−ℎ+𝜇_𝑑2𝜇+2ℎ+𝜎_𝐸+𝜎_𝐴,𝑛₀=𝜎_𝐴𝑏−ℎ+𝜇_𝑑×𝑏1−𝜂𝜎_𝐸−𝜇+ℎ+𝜎_𝐸𝜇+ℎ+𝜎_𝐴. (34)(35)(36)(37)(38) Therefore, the equilibrium solutions are 𝑁=0, corresponding to the absence of the population, and 𝑁=𝐾=𝑆, 𝐸=𝐴=𝐼=0, corresponding to the DFE when (16) holds. Finally, when 𝑛₄𝑁⁴_𝑒+𝑛₃𝑁³_𝑒+𝑛₂𝑁²_𝑒+𝑛₁𝑁_𝑒+𝑛₀=0, (39) we have 0<𝑁_𝑒<𝐾, 𝑆_𝑒,𝐸_𝑒,𝐴_𝑒,𝐼_𝑒>0.

We define the polynomial 𝑃𝑁_𝑒≡𝑛₄𝑁⁴_𝑒+𝑛₃𝑁³_𝑒+𝑛₂𝑁²_𝑒+𝑛₁𝑁_𝑒+𝑛₀, (40) which satisfies(1)𝑃(0)<0 for 𝜂 close to 1, and(2)𝑃(𝑁_𝑒)→∞ as 𝑁_𝑒→∞.

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, 𝑏>ℎ+𝜇. For the DFE (𝐾,0,0,0), the Jacobian matrix𝐽()=𝛽𝐾,0,0,0−𝑏+𝜇+ℎℎ+𝜇ℎ+𝜇−^𝐴𝛽𝜇𝛼𝑏−𝜇−ℎ−𝑏1−𝜂+ℎ+𝜇−^𝐼𝜇𝛼𝑏−ℎ−𝜇0−𝑏−𝜎_𝐸𝑏+𝛽1−𝜂^𝐴𝛽𝜇𝛼𝑏−𝜇−ℎ^𝐼𝜇𝛼𝑏−ℎ−𝜇0𝜎_𝐸−𝑏−𝜎_𝐴000𝜎_𝐴−𝑏+ℎ−𝜇_𝑑 (41)has the negative eigenvalue 𝜆=−𝑏+𝜇+ℎ (see (16)) and three eigenvalues of the matrix 𝐽₁=−𝑏−𝜎_𝐸𝑏+𝛽1−𝜂^𝐴𝛽𝜇𝛼𝑏−𝜇−ℎ^𝐼𝜎𝜇𝛼𝑏−ℎ−𝜇_𝐸−𝑏−𝜎_𝐴00𝜎_𝐴−𝑏+ℎ−𝜇_𝑑, (42) which are the zeros of the characteristic polynomial 𝜆³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀=0, with coefficients 𝑎₂=𝜎_𝐸+3𝑏+𝜎_𝐴−ℎ+𝜇_𝑑𝑎>0,₁𝜎_𝐸+𝑏2𝑏+𝜎_𝐴−ℎ+𝜇_𝑑−𝜎_𝐸𝛽_𝐸,𝑎𝜇𝛼𝑏−ℎ−𝜇₀=𝜎_𝐸+𝑏𝜎_𝐴+𝑏𝑏−ℎ+𝜇_𝑑−𝜎_𝐸𝑏−ℎ+𝜇_𝑑×𝑏+𝛽1−𝜂^𝐴−𝜎𝜇𝛼𝑏−ℎ−𝜇_𝐸𝜎_𝐴𝛽_𝐼.𝜇𝛼𝑏−ℎ−𝜇 (43)(44) Note that if 𝑎₀=0, then 𝑎₁>0, which means that there are two negative eigenvalues and one zero eigenvalue. Also, 𝑎₀=0 implies 𝑅₀≡𝑏−𝜇−ℎ𝜎_𝐸𝛽_𝐴𝑏−ℎ+𝜇_𝑑+𝛽_𝐼𝜎_𝐴𝜇𝛼𝑏−ℎ+𝜇_𝑑𝜎_𝐸+𝑏𝜎_𝐴++𝑏𝜇𝛼𝜎_𝐸1−𝜂𝑏𝑏−ℎ+𝜇_𝑑𝜇𝛼𝑏−ℎ+𝜇_𝑑𝜎_𝐸+𝑏𝜎_𝐴+𝑏=1. (45) In order to have three negative eigenvalues, we require 𝑎₀>0, which is equivalent to 𝑅₀<1.

Therefore, we have proved the following proposition.

Proposition 6. The equilibrium (𝐾,0,0,0) is stable for 𝑅₀<1 and unstable for 𝑅₀>1.

Condition (45) is equivalent to 𝐾<𝑏−ℎ+𝜇_𝑑𝜎_𝐸+𝑏𝜎_𝐴+𝑏−𝜎_𝐸1−𝜂𝑏𝜎_𝐸𝛽_𝐴𝑏−ℎ+𝜇_𝑑+𝛽_𝐼𝜎_𝐴≡𝐾_𝐶, (46) the critical carrying capacity. As a result, we have the following corollary.

Corollary 7. The equilibrium (𝐾,0,0,0) is stable for 𝐾<𝐾_𝐶 and unstable for 𝐾>𝐾_𝐶.

Remark 8. Note that if 𝛽_𝐴→0 and 𝜎_𝐴→∞ (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 𝑁_𝑒=𝐾 in (29) with 𝐸^′=0 and solving for 𝛽_𝐴(𝑏−ℎ+𝜇_𝑑)+𝛽_𝐼𝜎_𝐴, we obtain 𝛽_𝐴𝑏−ℎ+𝜇_𝑑+𝛽_𝐼𝜎_𝐴=𝜇𝛼𝑏−ℎ+𝜇_𝑑𝜎_𝐸+𝑏𝜎_𝐴+𝑏−𝜎_𝐸1−𝜂𝑏𝜎_𝐸𝑏−ℎ−𝜇≡𝛽_𝐶. (47) We have thus proved that 𝛽_𝐴(𝑏−ℎ+𝜇_𝑑)+𝛽_𝐼𝜎_𝐴=𝛽_𝐶 if and only if 𝑁_𝑒=𝐾. Also, 𝛽_𝐴(𝑏−ℎ+𝜇_𝑑)+𝛽_𝐼𝜎_𝐴=𝛽_𝐶 if and only if 𝑅₀=1. Moreover, at this critical rate of transmission, we have 𝐸_𝑒=𝐴_𝑒=𝐼_𝑒=0 and 𝑆_𝑒=𝐾.

As a result, we have the following two corollaries.

Corollary 9. Suppose 0<𝑁_𝑒≤𝐾. Then the following are equivalent: (i) 𝛽_𝐴(𝑏−ℎ+𝜇_𝑑)+𝛽_𝐼𝜎_𝐴=𝛽_𝐶, (ii) 𝑁_𝑒=𝐾, (iii) 𝐾=𝐾_𝐶, (iv) 𝑅₀=1, and (v) (𝑆_𝑒,𝐸_𝑒,𝐴_𝑒,𝐼_𝑒)=(𝐾,0,0,0).

Corollary 10. Suppose 0<𝑁_𝑒≤𝐾. Then(1)for 𝛽_𝐴(𝑏−ℎ+𝜇_𝑑)+𝛽_𝐼𝜎_𝐴<𝛽_𝐶, we have 𝑅₀<1; so the disease-free equilibrium (𝐾,0,0,0) is locally asymptotically stable;(2)for 𝛽_𝐴(𝑏−ℎ+𝜇_𝑑)+𝛽_𝐼𝜎_𝐴>𝛽_𝐶, we have 𝑅₀>1; so the disease-free equilibrium (𝐾,0,0,0) is unstable.

Remarks. Suppose (16) holds.(1)For 𝑅₀>1 and 𝜂→1, the EE exists. Since 𝑃(𝐾)>0 and 𝑃(0)<0, the result follows from applying the intermediate value theorem on the polynomial 𝑃 in [0,𝐾]. However, for 𝑅₀<1, the existence of the EE is not clear.(2)The condition 𝛽>𝛽_𝐶 implies 𝑁_𝑒<𝐾, 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 𝑅₀ is a decreasing function of 𝜂 while 𝛽_𝐶 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.

Let 𝐸_𝑓=𝐸_𝑚=𝐼_𝑓=𝐼_𝑚=𝐴_𝑓=𝐴_𝑚=0. Then (6)-(7) imply 𝑆^′_𝑓=𝑏𝑆_𝑓𝑆_𝑚𝑘𝑆_𝑓+𝑆_𝑚−𝜇_𝑓1+𝛼_𝑓𝑁𝑆_𝑓−ℎ_𝑓𝑆_𝑓,𝑆^′_𝑚=𝑏𝑆_𝑓𝑆_𝑚𝑘𝑆_𝑓+𝑆_𝑚−𝜇_𝑚1+𝛼_𝑚𝑁𝑆_𝑚−ℎ_𝑚𝑆_𝑚. (48) The equilibria thus satisfy 𝑏𝑆_𝑓𝑆_𝑚𝑘𝑆_𝑓+𝑆_𝑚−𝜇_𝑓(1+𝛼𝑁)𝑆_𝑓−ℎ_𝑓𝑆_𝑓=0,𝑏𝑆_𝑓𝑆_𝑚𝑘𝑆_𝑓+𝑆_𝑚−𝜇_𝑚(1+𝛼𝑁)𝑆_𝑚−ℎ_𝑚𝑆_𝑚=0. (49)(50) We have two cases.

Case 1. If 𝑆_𝑓=𝑆_𝑚=0, we have the zero equilibrium.

Case 2. The nonzero disease-free equilibrium is calculated as follows. First, we multiply (49) by 𝜇_𝑚/𝑆_𝑓 and (50) by 𝜇_𝑓/𝑆_𝑚. Then we solve the resulting system for 𝑧≡𝑆_𝑚/𝑆_𝑓 and 𝐾=𝑆_𝑓+𝑆_𝑚. We obtain 𝑧=𝑏𝜇_𝑓−𝑘𝜇_𝑓ℎ_𝑚−𝜇_𝑚ℎ_𝑓𝑏𝜇_𝑚+𝜇_𝑓ℎ_𝑚−𝜇_𝑚ℎ_𝑓,1𝐾=𝛼𝜇_𝑓𝑏𝑧𝑘+𝑧−ℎ_𝑓−𝜇_𝑓=1𝛼𝜇_𝑚𝑏𝑘+𝑧−ℎ_𝑚−𝜇_𝑚. (51) Equality of the two expressions for 𝐾 results from the solution for 𝑧. For the equilibrium densities for male and females, we obtain 𝑆_𝑓=𝐾1+𝑧,𝑆_𝑚=𝑧𝐾.1+𝑧 (52) A positive equilibrium requires 𝐾>0, which implies 𝑏𝑘+𝑧−ℎ_𝑚−𝜇_𝑚>0,or,equivalently,𝑏𝑧𝑘+𝑧−ℎ_𝑓−𝜇_𝑓>0. (53)
To investigate the stability of the disease-free equilibrium, we use the method of the next-generation operator [25, 26]. Set 𝑟₁=𝜎_𝐸+𝑏𝑧𝑘+𝑧,𝑠₁=𝜎_𝐸+𝑏,𝑟𝑘+𝑧₂=𝜎_𝐴+𝑏𝑧𝑘+𝑧,𝑠₂=𝜎_𝐴+𝑏,𝑟𝑘+𝑧₃=𝑏𝑧𝑘+𝑧−ℎ_𝑓+𝜇_𝑑,𝑠₃=𝑏𝑘+𝑧−ℎ_𝑚+𝜇_𝑑. (54)
The Jacobian matrix for the equilibrium (𝐸_𝑓,𝐸_𝑚,𝐼_𝑓,𝐼_𝑚,𝐴_𝑓,𝐴_𝑚,𝑆_𝑓,𝑆_𝑚) = (0,0,0,0,0,0,𝐾/(1+𝑧),𝑧𝐾/(1+𝑧)) has the form𝐽𝑀=𝐹−𝑉𝑂₁𝐽₂, (55) where𝑧𝑏𝐹=001−𝜂𝑘+𝑧+𝛽^𝑓_𝐴𝛾_𝑓𝐾𝛽1+𝑧^𝑓𝐴𝛾_𝑚𝑓𝐾𝛽1+𝑧^𝑓𝐼𝛾_𝑓𝐾𝛽1+𝑧^𝑓𝐼𝛾_𝑚𝑓𝐾𝑧𝑏1+𝑧001−𝜂𝑘+𝑧+𝛽^𝑚_𝐴𝛾_𝑓𝑚𝑧𝐾𝛽1+𝑧^𝑚𝐴𝛾_𝑚𝑧𝐾𝛽1+𝑧^𝑚𝐼𝛾_𝑓𝑚𝑧𝐾𝛽1+𝑧^𝑚𝐼𝛾_𝑚𝑧𝐾,𝑟1+𝑧000000000000000000000000𝑉=¹000000𝑠₁0000−𝜎_𝐸0𝑟₂0000−𝜎_𝐸0𝑠₂0000−𝜎_𝐴0𝑟₃0000−𝜎_𝐴0𝑠₃,𝐽,𝑂=000000000000₂=𝑏𝑧²(𝑘+𝑧)²−𝑏𝑧−1𝑘+𝑧𝑟1+𝑧³𝑏𝑘(𝑘+𝑧)²−1𝑟1+𝑧³𝑏𝑧²(𝑘+𝑧)²−𝑧𝑠1+𝑧³𝑏𝑘(𝑘+𝑧)²−𝑏−𝑧𝑘+𝑧𝑠1+𝑧³. (56)The 2×6 matrix 𝐽₁ turns out to be irrelevant to our calculations. We have 𝑉⁻¹=1𝑟₁0100000𝑠₁𝜎0000_𝐸𝑟₁𝑟₂01𝑟₂0𝜎000_𝐸𝑠₁𝑠₂01𝑠₂𝜎00_𝐸𝜎_𝐴𝑟₁𝑟₂𝑟₃0𝜎_𝐴𝑟₂𝑟₃01𝑟₃00𝜎_𝐸𝜎_𝐴𝑠₁𝑠₂𝑠₃0𝜎_𝐴𝑠₂𝑠₃01𝑠₃. (57) The threshold parameter 𝑅₀ is the spectral radius of the matrix 𝐹𝑉⁻¹, provided that all eigenvalues of 𝐽₂ have negative real part [26]. Thus 𝑅₀=𝜌𝐹𝑉⁻¹1=max2𝑏₁+12𝑏²₁−4𝑏₂,12𝑏₁−12𝑏²₁−4𝑏₂, (58) where 𝑏₁=𝑧1−𝜂𝑏𝜎_𝐸/(𝑘+𝑧)+𝜎_𝐸𝛽^𝑓_𝐴𝛾_𝑓(𝐾/(1+𝑧))𝑟₁𝑟₂+𝜎_𝐸𝜎_𝐴𝛽^𝑓_𝐼𝛾_𝑓(𝐾/(1+𝑧))𝑟₁𝑟₂𝑟₃+𝜎_𝐸𝛽^𝑚_𝐴𝛾_𝑚(𝐾𝑧/(1+𝑧))𝑠₁𝑠₂+𝜎_𝐸𝜎_𝐴𝛽^𝑚_𝐼𝛾_𝑚(𝐾𝑧/(1+𝑧))𝑠₁𝑠₂𝑠₃,𝑏₂=𝑧𝜎_𝐸1−𝜂𝑏/(𝑘+𝑧)+𝜎_𝐸𝛽^𝑓_𝐴𝛾_𝑓(𝐾/(1+𝑧))𝑟₁𝑟₂+𝜎_𝐸𝜎_𝐴𝛽^𝑓_𝐼𝛾_𝑓(𝐾/(1+𝑧))𝑟₁𝑟₂𝑟₃×𝜎_𝐸𝛽^𝑚_𝐴𝛾_𝑚(𝐾𝑧/(1+𝑧))𝑠₁𝑠₂+𝜎_𝐸𝜎_𝐴𝛽^𝑚_𝑖𝛾_𝑚(𝐾𝑧/(1+𝑧))𝑠₁𝑠₂𝑠₃−𝜎_𝐸𝑧1−𝜂𝑏/(𝑘+𝑧)+𝛽^𝑚_𝐴𝛾_𝑓𝑚(𝑧𝐾/(1+𝑧))𝑟₁𝑟₂+𝜎_𝐸𝜎_𝐴𝛽^𝑚_𝐼𝛾_𝑓𝑚(𝑧𝐾/(1+𝑧))𝑟₁𝑟₂𝑟₃×𝜎_𝐸𝛽^𝑓_𝐴𝛾_𝑚𝑓(𝐾/(1+𝑧))𝑠₁𝑠₂+𝜎_𝐸𝜎_𝐴𝛽^𝑓_𝐼𝛾_𝑚𝑓(𝐾/(1+𝑧))𝑠₁𝑠₂𝑠₃. (59)(60) Note that if 𝑏²₁−4𝑏₂≥0, then 𝑅₀=(1/2)𝑏₁+(1/2)𝑏²₁−4𝑏₂. In order to prove that 𝐽₂ has only eigenvalues with negative real part, we prove that tr(𝐽₂)<0 and det(𝐽₂)>0. First, we note the inequalities 𝑏𝑧²(𝑘+𝑧)²−𝑏𝑧𝑘+𝑧=−𝑏𝑘𝑧(𝑘+𝑧)²≤0,𝑏𝑘(𝑘+𝑧)²−𝑏𝑘+𝑧=−𝑏𝑧(𝑘+𝑧)²≤0. (61) From (53), we obtain tr(𝐽₂)<0.
Now, using (53), (61), as well as the fact that 𝑧 and all parameters are positive, we have det𝐽₂=−𝑏𝑘𝑧(𝑘+𝑧)²+11+𝑧𝑏𝑧𝑘+𝑧−ℎ_𝑓−𝜇_𝑓×−𝑏𝑧+𝑧𝑘+𝑧𝑏1+𝑧𝑘+𝑧−ℎ_𝑚−𝜇_𝑚−𝑏𝑘(𝑘+𝑧)²−11+𝑧𝑏𝑧𝑘+𝑧−ℎ_𝑓−𝜇_𝑓×𝑏𝑧²(𝑘+𝑧)²−𝑧𝑏1+𝑧𝑘+𝑧−ℎ_𝑚−𝜇_𝑚=𝑏𝑘𝑧(𝑘+𝑧)²𝑧𝑏1+𝑧𝑘+𝑧−ℎ_𝑚−𝜇_𝑚+11+𝑧𝑏𝑧𝑘+𝑧−ℎ_𝑓−𝜇_𝑓𝑏𝑧+𝑘+𝑧𝑏𝑘(𝑘+𝑧)²𝑧𝑏1+𝑧𝑘+𝑧−ℎ_𝑚−𝜇_𝑚+11+𝑧𝑏𝑧𝑘+𝑧−ℎ_𝑓−𝜇_𝑓𝑏𝑧²(𝑘+𝑧)²>0. (62) As a result, we have proved the following theorem.

Theorem 11. The DFE is stable for 𝑅₀<1 and unstable for 𝑅₀>1.

Note that the value of 𝑅₀ 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 𝛽^𝑠_𝐴→0, 𝛽=𝛽^𝑠_𝐼 for 𝑠∈{𝑓,𝑚},𝜂→1, 𝜎=𝜎_𝐸 and 𝜎_𝐴→∞. Then we have 𝑏₁→𝑎₁ and 𝑏₂→𝑎₂, where 𝑎₁=𝜎𝐾𝛽1+𝑧^𝑓𝛾_𝑓(𝜎+(𝑏𝑧/(𝑘+𝑧)))𝑟₃+𝑧𝛽_𝑚𝛾_𝑚(𝜎+(𝑏/(𝑘+𝑧)))𝑠₃,𝑎₂=𝜎𝐾⋅1+𝑧𝑧𝛽_𝑓𝛽_𝑚𝛾_𝑓𝛾_𝑚−𝛾_𝑓𝑚𝛾_𝑚𝑓(𝜎+(𝑏𝑧/(𝑘+𝑧)))(𝜎+(𝑏/(𝑘+𝑧)))𝑟₃𝑠₃. (63) Moreover, if 𝛾_𝑓=𝛾_𝑓𝑚 and 𝛾_𝑚=𝛾_𝑚𝑓, then 𝑅₀=1/2𝑎₁. If 𝑘=0 and all parameters are identical for males and females, then 𝑧=1 and 𝑅₀=𝛽𝛾𝜎𝑏−ℎ−𝜇𝜇𝛼(𝜎+𝑏)𝑏−ℎ+𝜇_𝑑, (64) which is the value of 𝑅₀ for the basic model (20). In the same way, if 𝑘=0 and all parameters for females and males are identical, then we will have 𝑧=1 and 𝑅₀ for the full model (58) matches 𝑅₀ 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 km² [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.

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 ℎ_𝑓=0.0922, ℎ_𝑚=0.282, and ℎ=0.1512 (estimated using 𝑃(𝑡)=𝑃₀𝑒^−ℎ𝑡, where 𝑃(𝑡) is the population after 𝑡 years and 𝑃₀ 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 𝜇_𝑓=0.1, 𝜇_𝑚=0.1111, and 𝜇=0.1053 (𝜇=1/𝐿, where 𝐿 is the lifespan).(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 𝐸₀=24. We assume that there are no infected animals with clinical signs (𝐼⁰_𝑓=𝐼⁰_𝑚=0). All other animals are considered susceptible (𝑆⁰_𝑓+𝑆⁰_𝑚=626). Also, we assume that the disease prevalence in males is initially twice that in females [32], so we have 𝐸⁰_𝑓=8 and 𝐸⁰_𝑚=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 𝛽^𝑓_𝐼=𝛽^𝑓_𝐴=𝛽_𝑓, 𝛽^𝑚_𝐼=𝛽^𝑚_𝐴=𝛽_𝑚, and 𝛽_𝑓=𝛽−0.01, 𝛽_𝑚=𝛽+0.01. See Table 2.

For the basic model, the endemic equilibrium (𝑆_𝑒,𝐸_𝑒,𝐼_𝑒)=(100,56,17) is stable and the basic reproduction number is 𝑅₀=5.407828. For the basic model with asymptomatic individuals, the endemic equilibrium (𝑆_𝑒,𝐸_𝑒,𝐴_𝑒,𝐼_𝑒)=(47,27,12,11) is stable and the basic reproduction number is 𝑅₀=7.026539. Finally, for the two-gender nonasymptomatic model, with 𝜂=1, 𝜎_𝐴→∞ and 𝛽_𝐴→0, the endemic equilibrium (𝑆_𝑓,𝑆_𝑚,𝐸_𝑓,𝐸_𝑚,𝐼_𝑓,𝐼_𝑚)=(103,73,15,8,4,2) is stable and the basic reproduction number is 𝑅₀=1.56363. These numbers indicate that the disease persists in the absence of control measures. Note how prevalence and 𝑅₀ 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 𝑅₀ and 𝐾 as functions of harvesting for the basic model. Any harvesting rate between ℎ₁=0.4130 and ℎ₂=0.4647 implies that the disease goes to extinction, but the population will survive. In Figure 3, we take three values for harvesting (ℎ=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.

(a)

(b)

(c)

Figure 3 

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

Figure 4 shows the graph of 𝑅₀ and 𝐾 as functions of harvesting for the basic model with asymptomatic individuals. Any harvesting rate between ℎ₁=0.4262 and ℎ₂=0.4647 implies that the disease goes extinct, but the population survives. In Figure 5, we take three values for harvesting (ℎ=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.

(a)

(b)

(c)

Figure 5 

Deer population in the WMU 150 area for the basic model with asymptomatic individuals when (a) ℎ<ℎ1, resulting in disease persistence; (b) ℎ1<ℎ<ℎ2, resulting in disease control; (c) ℎ>ℎ2, resulting in eradication of the population. (solid curve; susceptible; diamonds: exposed; squares: asymptomatic; dashed curve: infected).

Figure 6 shows the level curves 𝑅₀=1 and 𝐾=0 as functions of gender-specific harvesting parameters ℎ_𝑓 and ℎ_𝑚 for the full model. Any point (ℎ_𝑓,ℎ_𝑚) between 𝑅₀=1 and 𝐾=0 implies that the disease is controlled and the population survives. In Figure 7, we take three different points ((ℎ_𝑓,ℎ_𝑚)=(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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Deer population in the WMU 150 area for the full model. ((a) and (b)): If ℎ<ℎ1, then the disease persists, although it is stronger in males. ((c) and (d)): If ℎ1<ℎ<ℎ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 ℎ>ℎ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).

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 𝑅₀ on parameter variation for the basic model. We used Latin hypercube sampling and partial rank correlation coefficients (PRCCs) to identify which parameters 𝑅₀ 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 𝑅₀ 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 on 𝑅0. Parameters with PRCC >0 will increase 𝑅0 when they are increased, while parameters with PRCC <0 will decrease 𝑅0 when they are increased.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 9 

Monte Carlo simulations for 1000 runs drawn from parameter ranges using Latin hypercube sampling. If 𝛽 is sufficiently small or ℎ is sufficiently large, then the disease can be controlled. Note that ℎ has been constrained to ensure population survival. However, if ℎ is in the approximate range 0.6<ℎ<0.8, then the disease is controlled and the population survives.

It should be noted that variations in ℎ will change 𝑅₀ from values much greater than 1 to very low values, resulting in significant dependence of 𝑅₀ on ℎ in Figure 8. However, sufficiently high values of ℎ result in population extinction, so we constrained ℎ to avoid such harvesting. Hence the effect of realistic values of ℎ 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 𝛼_𝑖), 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 𝛼_𝑖=0 and multiplied the birth terms by a factor (1+𝜅𝑁)⁻¹. In the simplest case of a basic model, this approach gives a disease-free equilibrium with 𝐾=𝑏−ℎ−𝜇𝜅,𝜇+ℎ (65) and the basic reproduction number 𝑅₀=𝛽𝜎𝑏−𝜇−ℎ𝜅𝜇+ℎ𝜇+𝜎+ℎ𝜇+𝜇_𝑑. (66) 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 when 𝑅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.

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 of 𝑅₀ calculated for the models without gender structure is relatively high (in fact, most likely too high), the value for 𝑅₀ 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 reduce 𝑅₀ and increase the critical transmission rate. This reduction in 𝑅₀ 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 𝑅₀ 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 𝑅₀ 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 𝑅₀=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.

Acknowledgments

The authors thank Daniel Krewski and Tamer Oraby for technical discussions. M. Al-Arydah was partially supported by funding from PrioNet. R. J. Smith? is supported by an NSERC Discovery Grant, an Early Researcher Award, and funding from MITACS. F. Lutscher is supported by an NSERC Discovery Grant and startup funds from the University of Ottawa. For citation purposes, please note that the question mark in “Smith?” is part of his name.