The effect of predation on the dynamics of Chronic Wasting Disease in deer

TL;DR

This paper models the spread of Chronic Wasting Disease in deer and explores how predation by wolves and population control measures can effectively limit disease transmission, providing analytical and simulation insights.

Contribution

It introduces a bifurcation analysis of CWD dynamics and evaluates wolves as a biological control strategy through modified Gillespie simulations.

Findings

Wolves can effectively reduce CWD spread in cervid populations.

Population control measures like hunting aid in disease management.

Analytical probabilities of disease spread and elimination are provided.

Abstract

Chronic Wasting Disease (CWD) is a neurological disease impacting deer, elk, moose, and other cervid populations and is caused by a misfolded protein known as a prion. CWD is difficult to control due to the persistence of prions in the environment. Prions can remain infectious for more than a decade and have been found in soil as well as other environmental vectors, such as ticks and plants. Here, we provide a bifurcation analysis of a mathematical model of CWD spread in a cervid population, and use a modification of the Gillespie algorithm to explore if wolves can be used as an ecological control strategy to limit the spread of the disease in several relevant scenarios. We then analytically compute the probability that the disease spreads given one infected member enters a fully healthy population and the probability of elimination, given a fully susceptible population and remaining prions in the environment. From our analysis, we conclude that wolves can be used as an effective control strategy to limit the spread of CWD in cervid populations, and hunting or other means of lowering the susceptible population are beneficial to controlling the spread of CWD, although it is important to note that inferring biologically relevant parameters from the existing data is an ongoing challenge for this system.

Tables

Table 1. Table 1: a) birth rate at population =0 (=1.8 fawns per female deer per year) b) death rate = 0.1 (probability of survival per year = 0.9).

Parameter	Definition	Value	Reference
$b$	birth ratea	0.6 yr-1
	at population =0
$d$	death rateb	0.1 yr-1
$r=b-d$	effective birthrate	0.5 yr-1
${\gamma }_i$	direct contact infection rate	0.1 yr-1(infected animal density)-1	[30]
${\gamma }_e$	prion transmission rate	1yr-1(prion density)-1	wlog
${\mu }_i$	CWD death rate	0.6 yr-1	[30]
$ϵ$	Rate of excretion	0.1 yr-1(prion density)(deer density)-1	[30]
	of infectious material
${\mu }_e$	Rate of decay	0.1-0.2 yr-1	[30] [40]
	of infectious material
$K$	Carrying capacity	0-100(100 km${}_{}{}^2){}^{-1}$
$K^{\prime }$	=$K/\left(1-\frac{d}{b}\right)$	0-100(100 km${}_{}{}^2){}^{-1}$

Equations

$$\frac{dS}{dt}$$

$$\frac{dI}{dt}$$

$$\frac{dE}{dt}$$

$$\frac{dS}{dt}$$

$$\frac{dI}{dt}$$

$$\frac{dE}{dt}$$

$$\frac{dS}{dt}$$

$$\frac{dI}{dt}$$

$$\frac{dE}{dt}$$

$$S=K(1-\frac{{\rho }_s}{r}W),I=0,E=0.$$

$$S$$

$$I$$

$$E$$

$$H(\kappa )={\kappa }^2-\kappa ({\mu }_e^2+3{\mu }_e{\mu }_i+{\mu }_i^2)-{\mu }_e{\mu }_{i}r({\mu }_e+{\mu }_i)=0,$$

$$H_W$$

$$R_1:$$

$$R_2:$$

$$R_3:$$

$$R_4:$$

$$R_5:$$

$$R_6:$$

$$R_7:$$

$$R_1:$$

$$R_2:$$

$$R_3:$$

$$R_4:$$

$$R_5:$$

$$R_6:$$

$$R_7:$$

$$E(t)=\left\{\begin{array}[]{cc}\beta(1-\exp(-\mu_{e}t)),&0<t<t_{d}\\ \beta(1-\exp(-\mu_{e}t_{d}))\exp(-\mu_{e}(t-t_{d}),&t_{d}<t<\infty\end{array}\right.,$$

$$\frac{d{p}_s}{dt}=({\gamma }_{e}sE(t)+{\gamma }_i\frac{A_0}{A}s(i==1))(1-p_s),$$

$$p_s(t)=1-\exp (-{\gamma }_{e}s{\int }_0^{t}E(\sigma )d\sigma -{\gamma }_{i}s\frac{A_0}{A}(t_d-H(t_d-t)(t-t_d))),$$

$$p_s(\infty )=1-\exp (-{\gamma }_{e}s{\int }_0^{\infty }E(\sigma )d\sigma -{\gamma }_{i}s\frac{A_0}{A}{t}_d).$$

$$p_s(\infty )=1-\exp (-{\gamma }_{e}s\beta t_d-{\gamma }_{i}s\frac{A_0}{A}{t}_d).$$

$$p_s=\mu {\int }_0^{\infty }(1-\exp (-{\gamma }_{e}s\beta t_d-{\gamma }_i\frac{A_0}{A}{t}_d))\exp (-\mu t_d)d{t}_d=\frac{1}{\frac{\mu }{{\gamma }_{e}s\beta +{\gamma }_{i}s\frac{A_0}{A}}+1}.$$

$$\frac{\mu }{{\gamma }_{e}s\beta +{\gamma }_{i}s\frac{A_0}{A}}=\frac{A{\mu }_e({\mu }_i+{\rho }_{i}W)}{s{A}_0({\gamma }_eϵ+{\gamma }_i{\mu }_e)}.$$

$$\frac{d{p}_{s,0}}{dt}=-{\gamma }_{e}sE(t)p_{s,0}.$$

$$E(t)=E_0\exp (-{\mu }_{e}t),$$

$$\frac{d{p}_{s,0}}{dt}=-{\gamma }_{e}s{E}_0\exp (-{\mu }_{e}t)p_{s,0}.$$

$$\ln (p_{s,0}(t))=\frac{{\gamma }_e}{{\mu }_e}s{E}_0(\exp (-{\mu }_{e}t)-1),$$

Full text

The effect of predation on the dynamics of Chronic Wasting Disease in deer

Cody E. FitzGerald1 and James P. Keener2111deceased

1 Department of Engineering Sciences and Applied Mathematics, Northwestern University,

Evanston, IL 60208, USA

and

2 Department of Mathematics, University of Utah,

Salt Lake City, Utah, 84112 USA

Keywords— Chronic Wasting Disease, population ecology, control, management, bifurcation theory, wolves, dynamical systems, disease ecology

1 Abstract

Chronic Wasting Disease (CWD) is a neurological disease impacting deer, elk, moose, and other cervid populations and is caused by a misfolded protein known as a prion. CWD is difficult to control due to the persistence of prions in the environment. Prions can remain infectious for more than a decade and have been found in soil as well as other environmental vectors, such as ticks and plants. Here, we provide a bifurcation analysis of a simple mathematical model of CWD spread in a cervid population and use a modification of the Gillespie algorithm to explore if predators can be used as an ecological control strategy to limit the spread of the disease in several relevant scenarios. We then use several analytical probabilistic arguments to show that lowering the susceptible population is also beneficial to controlling CWD outbreaks. Finally, we consider a more complex model of CWD spread in a cervid population in which predators are assumed to be dynamic. Here, we find that, again, predators may be used to control CWD outbreaks, assuming they selectively prey upon infected cervids at a sufficiently high rate. From our analysis, we conclude that wolves may be used as an effective control strategy to limit the spread of CWD in cervid populations, and hunting or other means of lowering the susceptible population are beneficial to controlling CWD.

2 Introduction

Chronic Wasting Disease (CWD) is a prion-mediated neurodegenerative disease that has been observed in cervids, such as elk, deer, and moose, and was first described in 1980 [47]. Prions are misfolded and infectious proteins that can cause neurological disease, such as CWD. CWD is transmitted both directly between susceptible animals and indirectly through environmental exposure. Cervids can become infected through contact with saliva, blood, and other bodily fluids of an infected animal or through indirect environmental infections from water sources, plants, soil, or potentially insects. In recent years, wild populations of cervids have exhibited high levels of CWD prevalence. The 2024 Chronic Wasting Disease Surveillance Report from the Wyoming Game and Fish Department reported that the prevalence of CWD in some populations of mule deer is as high as 66.3% 222The report can be found here https://wgfd.wyo.gov/wyoming-wildlife/wildlife-disease-and-health/chronic-wasting-disease using the drop-down Report menu at the bottom of the website.. CWD has proven difficult to control or eradicate, due to the time scale associated with the persistence of CWD prions in the environment. Prions decay slowly and can remain latent and infectious for more than a decade [14, 45]. Many mathematical models have been developed to help understand the disease spread in deer populations, evaluate the efficacy of control strategies, and infer the modes of transmission. We briefly review several of the published models, but our summary should not be considered exhaustive.

One of the earlier discrete population models examined the spread of CWD and control strategies [31, 17], but the modeling assumptions, in particular frequency-dependent transmission, were reevaluated in the follow-up paper [39]. Miller and co-authors developed and attempted parameter estimation of competing SIR-type models that describe transmission of CWD in a population of cervids through different mechanisms—though at least one of the models considered appears to exhibit a lack of structural identifiability, which we briefly comment on in the Discussion [30]. They found support for indirect transmission of CWD using a model selection approach [30]. A different study also suggests that the disease will spread via indirect transmission due to the persistence of prions in the environment [4]. The model presented in [30] was further modified for further mathematical and control-based investigations.

Sharp and Pastor modified the model presented in [30] to include logistic growth in the susceptible population and performed a mathematical analysis of the system [40]. Wild and co-authors later extended the model presented in [30] to include the effects of predation from wolves [46] and analyzed the system in a non-oscillatory parameter regime, using parameter values estimated in [30] and other studies. An earlier report to the National Park Service by Hobbs also analyzed wolf-predation as a control strategy using a similar deterministic model in an oscillatory parameter regime [18]. Barbera and Pollino performed a mathematical analysis of a PDE [5] that was inspired by the model presented in [40]. Reyes and co-workers designed a spatiotemporal causal inference scheme to examine the effects of culling strategies using a PDE-based approach [25]. Potapov and co-authors used a complex age-structure differential equation model to examine the influence of harvesting strategies on population dynamics and disease suppression [36]. Vasylieva and co-authors performed a mathematical analysis of a differential equation model that describes CWD spread in a population of deer, including both indirect and direct transmission [43]. They also attempted a parameter estimation of the system, using upper and lower bounds for parameter estimation based on estimates of parameters found in [30]. A different study examined the influence of culling strategies to control CWD spread using a differential equation model that incorporated seasonality into a population model of cervids [35]. Others created an integrodifference equation model that incorporated direct and indirect transmission and long-distance dispersal from juvenile deer to understand CWD spread in a population of white-tailed deer [28]. Several groups have examined harvest strategies using mathematical models that delineate the deer population by sex [3, 2, 37]. Control investigations into CWD transmission in deer populations are summarized in [42].

Recently, agent-based models, statistical models, and machine learning approaches have been used to examine CWD transmission in deer populations. An agent-based model has been used to examine the spread of CWD in deer populations and possible management strategies to limit CWD [41]. A recent long-term statistical analysis used movement data collected from 596 deer to examine how seasonal behaviors, habitat selection, and home range size influence CWD transition in southwest Wisconsin [15]. Machine learning approaches have also been applied to predicting CWD incidence at the county level [1]. While machine learning approaches are often useful for accurate forecasting and prediction, they typically do not shed light on the underlying mechanisms that underpin the system and its dynamics, i.e., what drives CWD outbreaks and how such outbreaks may be controlled. Here, we attempt to study the control question using a simple, interpretable mechanistic model.

In this paper, we present and analyze a relatively simple model to understand the CWD spread in a population of deer and explore the role of predators, such as wolves, as an ecological control strategy. In our initial model, we treat the predator density as constant. Our model is inspired by previous mathematical models [40, 46, 18]. We first perform a bifurcation analysis to demonstrate the effect of predation on the dynamics of the disease and then simulate the model stochastically using a modified form of the Gillespie algorithm [16]. We show that the presence of wolves in all cases lessens the severity of a CWD outbreak and can lead to the elimination or at least significant reduction of the severity of the disease in the deer population. This is true even if the control intervention is applied after CWD has spread in the population. The mechanism of this reduction is easy to understand, as predation increases the death rate of infected animals and thereby decreases the amount of spread of prions from the infected animals. We then make analytic probabilistic arguments to show that reducing the susceptible population, not just the infected population, is also beneficial to limit future outbreaks. Finally, we study a more complex, realistic model, in which the predator is assumed to be dynamic. In this case, we find that predators can still be an effective control strategy for CWD outbreaks, but this intervention may not be expected to fully eradicate the disease. Our analysis shows that CWD outbreaks in deer may be controlled or limited by a combination of predator-based control strategies, which lower the infected population, and hunting or other harvesting strategies to limit the size of the susceptible population.

3 A Continuous Variable Population Model

The population dynamics model that we begin with was first presented and studied in [40], and is given by

[TABLE]

Here, $S$ is the density of susceptible and uninfected animals, $I$ is the density of infected animals, and $E$ is the density of prions in the environment (primarily in the soil). The meaning of the terms of the model is that $S$, in the absence of prions, grows according to logistic dynamics, but becomes infected at a law of mass action rate from prions. The infected animals have a death rate $\mu_{i}$. Prions are shed by infected animals at rate $\epsilon$, and degrade at rate $\mu_{e}$.

A slightly earlier model [46] included the effects of predation using the equations

[TABLE]

where $\rho_{s}(S,I)={S+I\over S+vI}$, $\rho_{i}(S,I)=v(1-c){S+I\over S+vI}$.

The analysis of the model in [40] showed that there are three regions of parameter space with distinct dynamic behaviors: a region with no infected animals or prions, a region in which the disease is endemic, but steady, and a region in which there are periodic (oscillatory) outbreaks of the disease, during which the animal population falls catastrophically low. They did not study the effect of predation on the dynamics of the disease. The model in [46] was used to study the effects of predation through the parameter, $\delta$, but failed to notice the existence of periodic solutions and the effect of predation on those oscillations. Separately, a rigorous mathematical analysis of a closely related model is given in [24].

The model considered here uses a slightly different density dependence on birth and death and a simplified predation effect, motivated by simple predator-prey interactions,

[TABLE]

where $W$ is the density of the predator (wolf) population, assumed here to be constant, and $\rho_{s}$ and $\rho_{i}$ are both constant, quantifying the effect of predation. Here $b$ is the zero population birthrate, $K^{\prime}$ is the population density at which the birthrate is zero, $d$ is the natural death rate of healthy animals, and $r=b-d$, $K=(1-{d\over b})K^{\prime}$, are both assumed to be positive. Further, the model incorporates the assumption that infected animals contribute to the carrying capacity, but not to the birth of susceptible animals. The disease is assumed to spread in one of two ways, through contact with prions in the environment (rate $\gamma_{e}$) and through direct contact between infected and susceptibles (rate $\gamma_{i}$). Prions are deposited into the environment because of shedding by infected animals (at rate $\epsilon_{i}I$) or death of infected animals (at rate $\epsilon_{d}\mu_{i}I$). Notice that the units for $\epsilon_{i}$ (yr*-1*(prion density)(deer density)-1) are different than those for $\epsilon_{d}$ ( (prion density)(deer density)-1). Here, our goal was to develop and analyze a model that is biologically motivated, contains both indirect and direct transmission pathways, contains monomial terms rather than rational terms for ease of use in analytic bifurcation theory, is simple to understand, and is interpretable. Parameter values used for this study are shown in Table 1, although we briefly mention possible concerns about structural identifiability of models used in [30] in Section 6. Without loss of generality, $\gamma_{e}$ can be taken as one by rescaling the units on $E$.

It is likely that $\rho_{s}<\rho_{i}$, since wolves prey primarily on aged and disabled deer, and the percentage of such among the susceptible population is certainly smaller than among the CWD-infected deer.

There are three possible steady solutions for this model: the extinct state with $S=0$, the disease-free state, and the disease endemic state. The disease-free state has

[TABLE]

This solution requires that $\rho_{s}W<r$, i.e., that predation of the wolves on healthy animals is not too large.

The disease endemic state is given by

[TABLE]

where $\gamma=\epsilon\gamma_{e}+\gamma_{i}\mu_{e}$. Clearly, this solution only exists provided $r(K\gamma-\mu_{e}\mu_{i})-(K\gamma\rho_{s}+\mu_{e}r\rho_{i})W>0$.

One can also show that the disease-free state is stable if $r(K\gamma-\mu_{e}\mu_{i})-(K\gamma\rho_{s}+\mu_{e}r\rho_{i})W<0$, and unstable if not. It is immediately apparent that the presence of predators decreases the likelihood of an endemic disease state existing and stabilizes the disease-free state. Indeed, if $W$ is large enough, there is only the disease-free state. Furthermore, in the disease-free state, the presence of predators decreases the $S$ population (see (10)), whereas in the endemic disease state, if it exists, the presence of predators increases the $S$ population, while it decreases the level of the infected population $I$ as well as the environmental load of prions, $E$ (see (11)).

Using XPP, Sharp and Pastor [40] showed that there is a curve in parameter space at which there is a Hopf bifurcation for the equations (1-3). It is a direct calculation (i.e., find the location in parameter space at which the Jacobian of the system (1-3) has purely imaginary eigenvalues [20]) to show that this curve is given by

[TABLE]

where $\kappa=K\epsilon\gamma$. Numerical simulations verify that there are stable periodic solutions for parameter values for which $H(\kappa)>0$. Since $H(\kappa)=0$ has a unique positive root, this implies that there are periodic solutions for $\kappa$ sufficiently large.

We can also find the curve of Hopf bifurcation points for equations (7-9). It is given by $H_{W}=0$, where

[TABLE]

where $\beta=M^{2}+(3\mu_{e}+r-\rho_{s}W)M+\mu_{e}^{2}$, $\gamma=\epsilon\gamma_{e}+\gamma_{i}\mu_{e}$, $M=\mu_{i}+\rho_{i}W$.

These features of predation are depicted in Fig. 1. It is shown that parameter space is divided into four regions: I, in which there is no disease; II, in which the disease is endemic; III, in which the disease is oscillatory; and IV, in which the deer population is extinct because of over-predation. The boundaries between these regions are shown plotted in $\rho_{i}W-K$ parameter space, showing the stabilizing effect of predators.

Not only does the presence of a predator decrease the likelihood of disease, but it also reduces the prion load of the environment and decreases the amplitude, hence the severity, of oscillatory outbreaks when they occur. This can be seen in numerical simulations, as demonstrated in Fig. 3. Specifically, in these plots are shown the numerical simulation of the ODE system for the fixed value of $K=30$ for different values of $\rho_{i}W=0,0.1,0.15,0.5$, placing these in regions II and III (see asterisks in Fig. 1). Noticeably, the effect of predation is to reduce the severity of the disease, at the cost of reducing the overall population size below its carrying capacity.

4 Stochastic simulation

Because deer populations are discrete and not too large, a continuous variable model may be inappropriate, and there might be insight to be gained from stochastic simulations of this process. To this end, we consider the following seven reactions:

[TABLE]

where $s$ and $i$ represent the number of susceptibles and infected, respectively, at any given time. Since we are taking the wolf population to be constant, it is not necessary to include wolf dynamics here. We need to get the reaction rates correct for the stochastic simulation. We start with some representative area $A$ that we want to consider in our simulation, and realize that the unit of densities is set by some area $A_{0}$, in this case, $A_{0}=100$km2. This means, for example, that if $K$ is the carrying capacity density, then $K$ is also the number size of the carrying capacity in an area of size $A_{0}$. Therefore, the size of the carrying capacity for an area of size $A$ is $K{A\over A_{0}}$. Now, let $s$ and $i$ be the integers that the densities $S$ and $I$ represent. Then, if $S$ is the density of susceptibles in the area $A_{0}$, the number of susceptibles in the area $A$ is $s=S{A\over A_{0}}$, and similarly for infected animals, $i=I{A\over A_{0}}$.

The reaction rates for the stochastic process are, therefore,

[TABLE]

Here, $E$ remains a continuous, not discrete, variable, and so is governed by the differential equation $\frac{dE}{dt}=\frac{\epsilon A_{0}}{A}i-\mu_{e}E$. We would like to use these reaction rates with the Gillespie algorithm [16] to do stochastic simulations. However, since $E$ is continuously changing in time, the Gillespie algorithm is not directly applicable. Instead, we use “Poisson thinning” to account for this fact, as follows[23, 34, 38]: Since between reactions, $E(t)$ is a monotone function, bounded between $E(t_{0})$ and $\epsilon\frac{iA_{0}}{\mu_{e}A}$, where $t_{0}$ is the last reaction time, we let $E^{*}={\rm max}[E(t_{0}),\epsilon\frac{iA_{0}}{\mu_{e}A}]$, and introduce a modified reaction $R_{3}^{*}$ (since $R_{3}$ depends on $E$), with reaction rate $r_{3}^{*}=\gamma sE^{*}$ (which is the fastest possible reaction rate), and then proceed with the usual Gillespie determination of next reaction time and next reaction. However, whenever the next reaction is to be $R_{3}^{*}$ at time $t^{*}$, the reaction is split into two, only one of which is implemented: a null reaction for which no change to $s$ or $i$ is made with probability $\frac{E^{*}-E(t^{*})}{E^{*}}$, and the reaction that reduces $s$ by one and increases $i$ by one with probability $\frac{E(t^{*})}{E^{*}}$.

Because $K$ is in the range of 0-100(100 km2)-1, so $A_{0}=100$km2, to have a reasonable number of deer, it must be that $A=1000$km2, at least. This also needs to take into account the territory size of a wolf pack, which can vary widely from less than 100 km2 to more than 4000 km2 [29]. In our stochastic simulations, we take $A=2000$km2 unless otherwise noted.

The first stochastic simulation shown (Fig. 4) is with no predators ($\rho_{i}W=0$). For these parameter values, the deterministic model has oscillatory solutions. For this model with the stochastic simulation, the population $s$ always goes extinct in finite time [10], and so to prevent that from happening, $s$ is artificially never allowed to become zero. Our goal is to maintain stochastically oscillating populations to assess the possibility of controlling CWD outbreaks using wolves. However, it should be noted that the system cannot realize the extinction state $s=i=E=0$ when we make this modeling choice. Fig. 4 (Right) shows the prion density trajectories for 50 independent trials, and it is seen here that most of the invasions with one infected cervid result in an oscillatory prion epidemic. The Figs. 5-7 show a sample trajectory (Left) and the prion density for 50 trajectories (Right) for the three values of $\rho_{i}W=0.1,0.5$, and 0.75.

There are several observations to make from these simulations. First, the probability of an epidemic getting started from one infected individual is a decreasing function of the number of wolves present at the time of the appearance of the infected individual (see Figure 8). Second, it could be that the diseased animals are eradicated, but the disease reemerges because of the remnant prions in the environment. However, third, there is a possibility that the disease is permanently eradicated even after it initially spreads, and this probability is an increasing function of the number of wolves present. For example, in Fig. 6 roughly 80% of the trajectories show disease elimination (with $\rho_{i}W=0.5$), while with $\rho_{i}W=0.75$ (Fig. 7) about 98% of the trajectories show elimination of the disease. We define disease eradication as $E<0.001$ at the end of the simulation. This, of course, is somewhat arbitrary. To our knowledge, the smallest density of prions required for infection is currently unknown and may depend on the transmission pathway. See also related analytic calculations in Sections 4.2-4.3. More synthetic data for the probability of disease eradication is shown in Fig. 8, where the probability of long-term disease survival as a function $\rho_{i}W$ is shown, determined using 200 trajectories for each value of $\rho_{i}W$.

4.1 Disease Control

For the previous simulations, it was assumed that there was a preexisting population of predators at the time the infected animal was introduced. Instead, it may be that wolves can be introduced only after the disease has been introduced and has begun to spread. In many regions333For example, see the 2024 Chronic Wasting Disease Surveillance Report from the Wyoming Game and Fish Department Report., the disease has already spread, and wildlife managers may wonder if introducing a population of predators could bring the disease under control. Here, we attempt to model this situation by introducing wolves once the infected animals comprise 5% of the population carrying capacity. The 5% threshold was chosen, as we felt that it was high enough that wildlife managers could feasibly detect in their annual surveys, whereas lower thresholds may be challenging to recover in a statistically accurate manner from hunter-harvested deer. Certainly, other thresholds could be used. The result without the introduction of predators is shown in Fig. 4. The result of this delayed introduction is summarized in Figs. 9-10. The phase portrait course of the disease for 50 trajectories with no intervention is shown in Fig. 9(Left), and with the introduction of wolves with $\rho_{i}W=0.8$ in Fig. 9(Right) and with $\rho_{i}W=1.5$ in Fig. 10(Left). In Fig. 10(Right) is shown the probability of disease survival as a function of time for these three scenarios. Apparent from these plots is that with the introduction of wolves, the severity of the disease is decreased in two ways, namely, the maximal level of infected animals is decreased, and the survival of uninfected animals is increased, giving improved opportunity for the herd to recover from the invasion.

4.2 Probability of Spread

We can get an analytical understanding of some of these results as follows: Suppose there is one infected cervid with a death rate $\mu(=\mu_{i}+\rho_{i}W)$ introduced into a healthy population of size $s$. What is the probability that prions shed by the one introduced infected individual and its interactions with healthy animals will lead to at least one other cervid becoming infected?

The density of prions resulting from shedding of one infected cervid satisfies the differential equation $\frac{dE}{dt}=\frac{\epsilon A_{0}}{A}-\mu_{e}E$ as long as the infected deer is alive $0<t<t_{d}$, and $\frac{dE}{dt}=-\mu_{e}E$ after its death.

[TABLE]

where $t_{d}$ is the time of death of the infected invader, and $\beta=\frac{\epsilon A_{0}}{\mu_{e}A}$. The probability $p_{s}$ of the disease spreading to at least one susceptible individual in a total population of $s$ by time $t$ is given by (ignoring all other births or deaths)

[TABLE]

from which it follows that the probability of spread by time $t$ is given by

[TABLE]

where $H$ is the Heaviside function, and the probability of eventual spread is

[TABLE]

A direct calculation using (30) yields

[TABLE]

Now, the probability that the infected animal dies at time $t_{d}$ is given by the probability density function $p_{d}(t_{d})=\mu\exp(-\mu t_{d})$, so the probability of spread, taking into account that the time of death $t_{d}$ is random, is

[TABLE]

In terms of original parameters,

[TABLE]

Notice that $p_{s}$ in (35) is a decreasing function of $\mu$ and an increasing function of $s$. Thus, predation of infected cervids is a beneficial thing, since anything that can increase $\mu$, the death rate of infected cervids, will decrease the probability of spread of the infection. In other words, the presence of wolves in a disease-free population helps to prevent the invasion of the disease by an immigrant infected animal. Further, observe that overpopulation of susceptible cervids increases the probability of the spread of the disease, while predation of the susceptible population, which decreases $s$ through $\rho_{s}W$, helps prevent the spread of the disease. Thus, the presence of predators is beneficial for at least these two reasons.

4.3 Probability of Elimination

Because prions survive in the environment, even if all infected animals are eliminated by some control mechanism, there is the possibility that the disease will reemerge some time later. To calculate this probability, let $p_{s,0}(t)$ be the probability that there are $s$ susceptible deer and no infecteds at time $t$, but there is some contamination of the environment with prions, $E$. It follows that

[TABLE]

(Here, again, we are ignoring births and deaths of healthy animals). $E(t)$, of course, is the density of prions in the environment, which, as long as there are no infecteds, satisfies $\frac{dE}{dt}=-\mu_{e}E(t)$. Consequently,

[TABLE]

where $E_{0}$ is the initial density of prions at time $t=0$, so that

[TABLE]

It follows that

[TABLE]

and

[TABLE]

is the probability that there will not be another outbreak. Obviously, the probability that there will not be another outbreak goes to zero as ${\gamma_{e}\over\mu_{e}}sE_{0}$ becomes large. This once again shows the risk associated with overpopulation of susceptible deer (large $s$), and the benefit of having wolves or other harvesting measures to help control the population size of healthy deer, even when there are no infected animals present. This is borne out by a recent statistical study that showed that deer populations that experienced high harvest pressure had lower prevalence of CWD [32].

5 Control in the full dynamical model

In our analysis so far, we took $W$ to be constant, to keep the analysis more analytically tractable (see Section 3). Here, we explore the control implications of CWD outbreaks when $W$ is dynamic. Consider the following augmented form of Eqs. (7)-(9)

[TABLE]

Here, we assume that $W$ grows at some rate $r_{w}$ and the population is governed by a carrying capacity, $K_{w}$. As in Eq. (7), we take $r_{w}=b_{w}-d_{w}$, where $b_{w}$ is the $d_{w}$ are zero population birthrate for wolves and the natural mortality rate for wolves, respectively. We take $b_{w}=0.5$ and $d_{w}=0.2$, as we thought it reasonable that wolves should have a lower birthrate than deer but a higher natural mortality rate than deer. As before in Eq. (7) , $K_{w}=K_{w}^{\prime}(1-\dfrac{d_{w}}{b_{w}})$. We selected the carrying capacity to be $K_{w}=3$ (100 km${}^{2})^{-1}$. Clearly, it must be true that $K_{w}<K$, and this value seemed reasonable with estimates of wolf densities [29], though the reported densities have a large spread, and other choices for this parameter could have been made. We also assume predation increases the population size, through some conversion efficiency, $\alpha$, but very weakly. We assume that $\alpha<<1$, as most predation events should not significantly contribute to the growth of the wolf population. In Section 3, we assumed that $\rho_{s}<\rho_{i}$, and set $\frac{\rho_{s}}{\rho_{i}}=0.5$, as predators may selectively target weaker prey. Here, we set $\rho_{i}=\omega\rho_{s}$, where $\omega\geq 1$ (Eqs. (42)-(45)). Biologically, $\omega$ represents the degree of the predation bias towards infected deer over susceptible deer. We will treat $\omega$ as our bifurcation parameter. At steady state, Eqs. (42)-(45) have five physical equilibria:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

Due to the nonlinear structure of Eqs. (42)-(45), computing analytic expressions that describe the location of the Hopf bifurcation in parameter space becomes intractable. The analytic expression becomes so large and unwieldy that computer algebra systems struggle to compute the curve. However, we can still investigate the system numerically to determine if predators can be used to control a CWD outbreak when the system is governed by Eqs. (42)-(45). We see that Eqs. (42)-(45) emit an oscillatory solution, which biological represents oscillatory CWD outbreaks, for the parameter values listed in Table 1, $\rho_{s}=0.01$, $r_{w}=0.3$, $K_{w}=3$, and sufficiently low $\alpha$ and $\omega$ ($\alpha=0.05$, and $\omega=1$, see Figure 11, top). Increasing $\omega$ dampens the oscillation (see Figure 11, middle), and further increasing $\omega$ drives the system into the endemic state (Eq. 50) where oscillatory outbreaks are eliminated in the deterministic model. Computationally, the transition to the endemic state occurs around $\omega\approx 5$ for the parameter values we used in our study. Moving the system from oscillatory CWD outbreaks to the endemic state comes at the cost of reducing the overall deer population well below the carrying capacity. It may be possible for wolves or other predators to control or at least temper CWD outbreaks, as we see here, but $\omega$ needs to be sufficiently high for this strategy to work. Biologically, this means that predators need to selectively prey upon infected deer at a sufficiently high rate as compared with the predation rate of susceptible deer. As before, we also study a stochastic version of the system, as the populations are expected to be relatively small.

We stochastically simulate the system for variable $\omega$, using the modified Gillespie algorithm (see Section 4). As wolves are no longer constant, we need to modify existing reactions and include several new reactions. In particular, we have

[TABLE]

where $s$ and $i$ represent the number of susceptibles and infected, respectively, at any given time, as before, and similarly, $w$ represents the number of wolves at any given time. These reactions occur at the rates:

[TABLE]

As before (see Sections 3-4), we see that it is possible to use wolves to force the system into the endemic disease state, but the effectiveness of the strategy crucially depends on the parameter $\omega$, which describes the selectively of predation on the infected deer population (see Figure 12). We also see that when the predator is assumed to be dynamic, it may be less likely that wolves eradicate the disease completely, and instead we may expect wolves to push the system into the endemic state (Eq. 50) rather than the disease-free state (Eq. 48) (see Figure 13). To ensure that the disease is eradicated when the predator is dynamic requires $\omega>>1$, which may not be biologically reasonable (see Figure 13). Further biological research is needed to determine estimates of $\omega$.

6 Discussion

Here, we performed a bifurcation analysis of a simple deterministic mathematical model (Eq. (7)-(9)) that describes the direct and indirect transmission of CWD in a population of deer subject to predation by wolves, which we originally assume are constant (see Figures 1-3, Section 3). We then stochastically simulated the system using the Gillespie Algorithm and applied Poisson-thinning to the prion density, as it varies continuously in time (see Figures 4-8). This model investigation demonstrates that introducing predators, such as wolves, may be an effective way to control diseases such as CWD, for which there are no other known effective control strategies. We also demonstrated that this control strategy works even if CWD is already prevalent in the population, as is the case for many cervid populations in the American West and Upper Midwest (see Figures 9-10, Section 4.1). We also compute the probability of the disease spreading to at least one susceptible individual in a total population of $s$ by time $t$ (Section 4.2) and the probability that no future outbreak will occur given a completely susceptible population and $E$ in the environment (Section 4.3). These calculations make clear the need to control not only the size of the infected population, but also the size of the susceptible population.

Finally, we study a more complex version of the system (Eqs. (42)-(45)) in which the predator is assumed to be dynamic, rather than constant, as in (Eqs. (7)-(9)). Our analysis of the fully dynamic model (Eqs. (42)-(45)) shows that wolves still may be effective for controlling CWD, but the strategy depends on strong selective predation of infected deer over susceptible deer, which we describe using the parameter $\omega$ (see Figures 11, 12, and Section 5). For the parameter regime we used, we found that full eradication of the disease required $\omega>>1$, which may not be biologically reasonable, and so one might expect that introducing wolves may ultimately push the system towards the endemic state (Eq. 50) rather than the disease-free state (Eq. 48). Further biological research is needed to determine estimates of $\omega$, and specifically if $\omega>1$, though some evidence supports the idea that some predators may selectively target infected deer at least compared with human hunters [21]. A recent study using a more complex, age-structured deterministic model also found that predators may help to limit CWD outbreaks, but the result also depended on predation selectivity [8]. Taken together, our analysis highlights wolves, or other predators, as a possible effective control strategy for CWD and makes clear the role for human-hunting or other harvesting measures to limit the susceptible population size.

Our analysis has several limitations relating to the parameter values, model complexity, and modeling assumptions. We found it somewhat difficult to find published parameter values to parameterize our mathematical models with. Ultimately, we used parameter estimates from [30] and previous modeling efforts [40], and otherwise used intuition to make reasonable guesses about the magnitude of the parameters. We noticed that some of the models used to generate parameter estimates may lack structural identifiability [30]. Consider the “Indirect + Direct” model examined in [30], which has a similar form to Eq. (7)-(9). The “Indirect + Direct” model has the form

[TABLE]

The authors attempted to fit this model using observations of $C$ from two data sets of CWD outbreaks. However, notice that the system has a scaling symmetry 444We used STRIKE-GOLDD to find the scaling symmetry [27]. The MATLAB implementation we used can be found here https://github.com/afvillaverde/strike-goldd. The “Indirect” model in [30] may also have the same scaling symmetry. Other models were not examined. , which causes a structural identifiability issue (see {\color[rgb]{0,0,1}\psi} below). The structure of the model prohibits the discovery of a unique parameter estimate because the dynamics of the system are the same for any choice of {\color[rgb]{0,0,1}\psi}.

[TABLE]

A structural identifiability analysis 555We used StructuralIdentifiability.jl to perform the identifiability analysis of the system [11]. The implementation we used can be found here https://docs.sciml.ai/StructuralIdentifiability/stable/. For a structural identifiability tutorial in an epidemiological context, see [26]. of the model structure reveals that the parameters $\gamma$ and $\epsilon$ are structurally non-identifiable when only $C$ is observed, though the other kinetic parameters are globally identifiable. A recent movement ecology analysis showed that susceptible and infected deer can be detected from GPS tracks [6], possibly providing a cheap method for temporal measurement of $S$ and $I$. Given this finding and temporal estimates of $S$ and $I$ from state CWD surveillance programs, it is tempting to think that if one could measure $S$, $I$, and $C$, instead of just $C$, the structural identifiability issue may be resolved. This is not true. The parameters $\gamma$ and $\epsilon$ are still globally non-identifiable even if $S$, $I$, and $C$ are observed, but $E$ is unobserved. However, if one is able to observe $C$ and $E$, but not $S$ and $I$, all kinetic parameters are globally identifiable, but measuring $E$ is a key challenge for the system. The prions responsible for CWD are known to exist in soil, and can be measured from soil samples [22]. In practice, estimating $E$ through soil samples would likely be expensive, infeasible on a large scale, and may miss some of the sources of the prions in the environment. The prions responsible for CWD can be found in multiple environmental reservoirs beyond soil. It was recently shown that ticks can harbor levels of prions that could cause transmission of CWD [19]. Prions were also recently shown to exist in plants at relevant levels for transmission of CWD [9]. For these reasons, getting estimates of the prion density will be difficult at best or otherwise infeasible.

An alternative approach to ensure that unique parameters can be found is to reparameterize the model to remove the structural identifiability issue ahead of parameter estimation [27, 11] (for an example in a different modeling context, see [13]). In general, this process can impact the interpretability of state variables of the model or rate parameters. One way to resolve the issue in the case above is to take $\psi=\gamma$, which rescales the parameters $\epsilon$ and the state variable $E$, is enough to break the scaling symmetry and ensure that the model is structurally identifiable if only $C$ can be observed666This can be verified independently using StructuralIdentifiability.jl [11].The reparameterized “Indirect + Direct” model would then have the form

[TABLE]

where $\tilde{E}=\gamma E$, and $\tilde{\epsilon}=\gamma\epsilon$.

Of course, whether unique parameter estimates are needed depends on how the model will be used. For example, in the analysis presented in this paper, unique parameter estimates are not needed, as we are concerned with qualitative changes in the system’s behavior. However, in different modeling situations (i.e., generating precise forecasts using a model, or evaluating the impact of managerial controls quantitatively), unique parameter estimates may be important, or at least modelers should be conscientious of possible identifiability issues, especially in the context of wildlife management decisions.

We sought to limit the complexity of the model considered here, largely to keep the analytic calculations for the bifurcation analysis tractable (Section 3) and allow the model to be easily interpretable. In actuality, some of the model parameters, for example, the carrying capacities, $K$ and $K_{w}$, may vary over time, as carrying capacities are a dynamic interaction between populations and the environment. We assumed constant parameters to retain simplicity in our model and allow for analysis of the system. Additionally, given the newer findings about the wide range of environmental reservoirs in which prions are found [9, 19], our model may be overly simple. Both ticks and plants have complex, multi-year seasonal life cycles that exhibit phases of dormancy that are cued by environmental signals, potentially inducing a periodic and stochastic component to indirect environmental infection. These complex dynamics are beyond the scope of our model, but may be important to understanding the spread of CWD in deer at the level of an ecosystem as more is understood about the relative importance of different modes of transmission.

In our analysis, we made several assumptions related to the implementation of the modified Gillespie algorithm and our modeling framework. We took the rate at which susceptible deer, $s$, and wolves, $w$, are born as in Eqs. (62) and (71), respectively. Note that this is a relatively common approximation, but not exact, as it allows for the same “molecule” to react with itself. In Eq. (71) we take the rate as $R_{10}=b_{w}w(1-{A_{0}\over K_{w}^{\prime}A}w)$ instead of $R_{10}=b_{w}w(1-{A_{0}\over K_{w}^{\prime}A}(w-1))$ and similarly for Eq. (62). Near the carrying capacity, this is a fine approximation, but less accurate for smaller populations. Relatedly, the population $s$ will always go extinct in finite time for this model and stochastic simulation [10], and so to prevent that from happening, $s$ is artificially never allowed to become zero. Our goal was to maintain stochastically oscillating populations to assess the possibility of controlling CWD outbreaks using wolves. However, it should be noted that the system cannot realize the extinction state when we make this modeling choice, though we would not expect that wolves would hunt a large deer population to extinction.

In terms of our modeling framework, we chose to work with an ODE model and did not consider the role of possible spatial inhomogeneity. However, the role of space may be important in understanding the spread of CWD in deer. Currently, little is known about the spatial distribution of CWD prions. One would imagine the spatial distribution of the prion is initially somewhat “patchy” in space because an infected animal carcass contains a high concentration, which would then diffuse in the local environment through a variety of mechanisms. For instance, it is known that prions can remain infectious after passing through a crow’s digestive system [12, 44] and a coyote’s digestive system [33], which may lead to the spread of $E$ in the environment. The same may be expected of other predators, such as a wolf, and so it may be expected that wolves would spread the prions throughout the environment through feces. However, a different study showed that mountain lions may neutralize CWD prions during digestion [7]. As more is understood about the spread or elimination of the prion by predators, one might consider modifying Eq. (44) in our simple model to include a term that describes the spread ($\xi>0$) or elimination ($\xi<0$) of $E$ by the predator, such as

[TABLE]

but more biological research is needed to determine how dominant a factor the spread or elimination of the CWD prion by predators may be. Currently, less is known about the transmission of the disease through insects and plants, but these are other mechanisms by which prions can spread throughout the environment in a spatial sense. More research is needed to understand the spatial distribution of prions in the environment and how this contributes to the indirect transmission of CWD. As more is understood about the dynamics of environmental reservoirs of prions, such as ticks and plants, the spatial distribution of prions in the environment, and, more broadly, an ecosystem-scale understanding of CWD transmission, our model can be modified to understand these dynamical processes, or it may be more appropriate to use a spatially explicit model, such as a PDE.

This modeling analysis shows that predators may be used as an effective control strategy for managing the spread of CWD in deer populations and highlights a role for hunting or other harvesting of the susceptible population. We are careful to note that for wolves, or other predators, to act as an effective control strategy for CWD outbreaks, strong selective predation of infected deer over susceptible deer is required in our analysis. However, this modeling work should be understood as a proof-of-concept and not used directly for management purposes, given the uncertainty that exists in the parameters of the system. Given the available biological data for this system and associated identifiability challenges, care is needed when using mathematical models for wildlife management purposes.

7 Funding

C. E. F. is supported in part by The James S. McDonnell Foundation Postdoctoral Fellowship Award in Complex Systems (https://doi.org/10.37717/2020-1591) and by the NSF-Simons Center for Quantitative Biology at Northwestern University (NSF: 1764421 and Simons Foundation/SFARI 597491-RWC).

8 Author Contribution

C.E.F. and J.P.K. designed the research. J.P.K. performed the model analysis and simulations. C.E.F. reviewed the analysis and performed the identifiability analysis. J.P.K. and C.E.F. wrote and edited the paper.

9 Code

The code used for this analysis can be found here https://github.com/cefitzg/cwd_code.git. ChatGPT-5.5 was used to check the code for bugs, and Grammarly was used to catch typos in the manuscript.