Spatiotemporal epidemic models for rabies among animals.

Rabies is a serious concern to public health and wildlife management worldwide. Over the last three decades, various mathematical models have been proposed to study the transmission dynamics of rabies. In this paper we provide a mini-review on some reaction-diffusion models describing the spatial spread of rabies among animals. More specifically, we introduce the susceptible-exposed-infectious models for the spatial transmission of rabies among foxes (Murray et al., 1986), the spatiotemporal epidemic model for rabies among raccoons (Neilan and Lenhart, 2011), the diffusive rabies model for skunk and bat interactions (Bonchering et al., 2012), and the reaction-diffusion model for rabies among dogs (Zhang et al., 2012). Numerical simulations on the spatiotemporal dynamics of these models from these papers are presented.

Introduction

Rabies is an acute, viral, and fatal zoonotic disease to mammals. It remains an important threat to public health and a concern of wildlife management worldwide. Human rabies still causes thousands of deaths annually in Asia and Africa (Fooks et al., 2014, Wunner and Briggs, 2010), and dogs are responsible for most of these deaths (CDC, 2011, WHO, 2010). Rabies virus is present among various mammal species, including red fox and raccoon dog in Europe; raccoon, red fox, skunk, and insectivorous bats in North America; domestic dogs, insectivorous and vampire bats in South America; and domestic dogs, bat, Chinese ferret badger, raccoon dog, rat, fox, and wolf in Asia (Sterner and Smith, 2006, Wang et al., 2014).

Rabies emerged in Eastern Europe after World War II and spread westward through the 1980s. The spread of rabies among foxes has inspired extensively studies on rabies, including mathematical modeling on analyzing the epidemiological characteristics and transmission dynamics of rabies and designing useful control measures. In a pioneer paper, Anderson, Jackson, May, and Smith (1981) developed a deterministic model consisting of three subclasses of fox, susceptible, infectious and recovered, to explain epidemiological features of rabies in fox populations in Europe. A susceptible, exposed, infectious, and recovered (SEIR) model was proposed by Coyne, Smith, and McAllister (1989), and lately was also used by Childs et al. (2000), to predict the local dynamics of rabies among raccoons in the United States. Clayton, Duke-Sylvester, Gross, Lenhart, and Real (2010) and Ding, Gross, Langston, Lenhart, and Real (2007) considered the optimal control of SEIRS models which describe the population dynamics of rabies in raccoons. Dimitrov, Hallam, Rupprecht, Turmelle, and McCracken (2007) presented a model for the immune responses to a rabies virus in bats and George et al. (2011) presented a mathematical model parametrized with data on rabies in big brown bats in Colorado. Besides these deterministic models, discrete models (Allen et al., 2002, Artois et al., 1997), individual-based models (Rushton, Shirley, MacDonald, & Reynolds, 2006), and stochastic models (Russell et al., 2006, Smith et al., 2002, Smith and Wilkinson, 2003) have also been employed to study the transmission dynamics of rabies. We refer to reviews by Sterner and Smith (2006) and Panjeti and Real (2011) for more references on different rabies models.

There have been some studies on modeling canine and human rabies, see, for example, Carroll et al., 2010, Hampson et al., 2007, Zinsstag et al., 2009, etc. Ruan (2017) reviewed some recent studies on modeling the transmission dynamics of human rabies in China by considering different characters and aspects.

To study the westward spread of rabies among foxes in Europe, mathematical models described by partial differential equations have been proposed. Källén (1984) and Källén, Arcuri, and Murray (1985) studied rabies transmission in fox population by differential equations with diffusion and used the deterministic model to simulate rabies epizootic in foxes crossing continental Europe and proved the existence of traveling waves. Murray, Stanley, and Brown (1986) and Murray and Seward (1992) also considered foxes rabies, calculated the speed of propagation of the epizootic front and the threshold for the existence of an epidemic, and quantified a means to control the spatial spread of the disease. Since then, spatiotemporal models have been developed to study the spatial spread of rabies among other animals. In this paper we give a mini-review on some reaction-diffusion models describing the spatial spread of rabies among animals. More specifically, we introduce the susceptible-exposed-infectious models for the spatial transmission of rabies among foxes (Murray et al., 1986), the spatiotemporal epidemic model for rabies among raccoons (Neilan & Lenhart, 2011), the diffusive rabies model for skunk and bat interactions (Borchering et al., 2012), and the reaction-diffusion model for rabies among dogs (Zhang et al., 2012). Numerical simulations on the spatiotemporal dynamics of these models from these papers are presented.

Spatial spread of rabies in fox (Murray et al., 1986)

Murray et al. (1986) studied the spatial spread of rabies among foxes and examined the rabies epidemic, started in 1939 in Poland and moved steadily westward at a rate of 30–60 km per year. The basic spatial model of Murray et al. (1986) is an extension of the ODE model developed in Anderson et al. (1981) by including the spatial spread of the disease, which is caused by the random dispersal of rabid foxes.

Let and denote densities of susceptible, infected but non-infectious, and infectious foxes, respectively, in the space-time coordinate The basic model assumptions made by Murray et al. (1986) are as follows: (i) The dynamics of the fox population in the absence of rabies is approximated by the logistic growth law with the birth rate the intrinsic death rate and the environmental carrying capacity . The seasonality of births and food supply are neglected. (ii) Rabies is transmitted from rabid to susceptible fox: interspecies transmission is neglected. Susceptible foxes become infected at an average rate per head which is proportional to the number of rabid foxes present. (iii) Infected foxes become infectious at an average rate per head where is the average incubation time. (iv) Infectious foxes die at an average per capita rate α ( is the average duration of clinical disease). (v) Infected and infectious foxes continue to pressure on the environment and to die of cause other than rabies, but they have a negligible number of healthy offspring. (vi) Foxes are territorial and divide their territories up into non-overlapping ranges. (vii) Rabies is transmitted by direct contact (usually by biting) between foxes. (viii) Rabies acts on the central nervous system inducing behavioral changes in foxes. About half of infected foxes have furious rabies and exhibit the ferocious symptoms typically associated with the disease, while with the rest the virus affects the spinal cord, causing gradual paralysis. Foxes with furious rabies may become aggressive and confused, losing their sense of direction and territorial behavior, and wandering randomly. So a diffusion term is added to the equation for the infectious foxes.

The base spatial model takes the following formwhere is the total fox population and D is the diffusion coefficient. The term in each equation represents depletion of the food supply by all foxes. The dimensional parameters in (1) are given in Table 1 and are taken from Murray et al. (1986).

Table 1

Description of parameters in model (1).

Parameters	Value	Description
a	1 fox per year	Fox average birth rate
b	0.5 per year	Fox average intrinsic death rate
1/α	5 days	Duration of clinical disease
1/σ	28 days	Incubation time
K	0.25−4.6 fox/km2	Fox carrying capacity
β	80 km2 per year	Rabies transmission coefficient
D	200 km2 per year	Fox diffusion coefficient

When model (1) becomes the spatial homogeneous model proposed by Anderson et al. (1981) who found that when rabies is introduced into a stable population of healthy foxes these equations predict three possible behaviors. By considering the basic reproduction number, they obtained a critical value of the carrying capacity of the system given by (see also Wang & Zhao, 2012)

If then rabies eventually disappears and the population returns to its initial size If then rabies become endemic and the population oscillates about a positive steady state These oscillations are damped if K is not too much larger than in which case the system approaches the steady state , whereas if K is sufficiently large the system approaches a limit cycle oscillating periodically about the steady state This critical value was estimated to be between 0.2 and 1.0 foxes km−2 (Anderson et al., 1981).

The case where indicates the persistence of the disease in a spatially homogeneous setting. The spatial diffusion then propagates the disease so that a small localized introduction of rabies evolves into a traveling wave with a certain wave speed, that is, a solution with with the wave variable In (Murray et al., 1986), graphs of the propagation of the initial rabies outbreaks are given for the model with foxes km−2 and foxes km−2, approximately the carrying capacities for foxes in continental Europe and England, respectively. The traveling waves in Fig. 1 (a) consists of the rabies front, in which the largest number of foxes die from the disease, followed by an oscillatory tail, in which each successive outburst of rabies is smaller than the preceding one. The oscillations gradually approach constant, non-zero values with the rabid and infected fox population zero. Fig. 1(b) illustrates the fluctuations in fox density for a traveling wave with parameters appropriate for England. The diffusion coefficient D is estimated to be 60 km2 per year, using the average territory of a fox and the mean time such a fox stays in its territory. This yields the minimal wave speed near 50 km per year, in good agreement with the empirical data from Europe (Murray et al., 1986).

Fig. 1

Typical fluctuations in the fox populations due to the passage of a rabies epidemic wave as calculated from model (1. (a) The fox density in the uninfected region is taken to be at a carrying capacity of 2.0 foxes km−2. (b) The fox density in front of the epidemic is at a carrying capacity of 4.6 foxes km−2. Figures were adopted from Murray et al. (1986).

Murray and Seward (1992) generalized model (1) to include a population of immune foxes and find that this aspect has little effect on the propagation speed of the initial wave of the rabies epidemic but it affects the behavior of the periodic outbreaks associated with the oscillating tail of the wave. They also used the modified model to estimate the width of a rabies break which would be required to contain the epidemic.

Optimal vaccine distribution for rabies among raccoons (Neilan & Lenhart, 2011)

One of the strategies in mitigating the spread of rabies among raccoons in the eastern US and Canada is to distribute oral rabies vaccine baits by hand and by aircraft. After eating a bait, a healthy raccoon will develop antibodies in weeks that will provide protection if the raccoon is exposed to an infectious raccoon. Ding et al. (2007) used optimal control theory to characterize optimal vaccination levels for a raccoon model with discrete time and space and Asano, Gross, Lenhart, and Real (2008) performed similar studies for continuous time and discrete space. Neilan and Lenhart (2011) used a reaction-diffusion model with continuous time and space domain to explore the influence of realistic heterogeneous spatial domains on raccoon vaccine distribution. Landscape features, such as rivers and heavy forest cover, and long-distance translocation (LDT) of raccoons are known to perpetuate irregular dynamics in the rabies wave front. They compared optimal strategies of vaccine bait placement on a homogeneous spatial domain with those on a heterogeneous domain incorporating a river, forest cover, and LDT.

Let represent a rectangular grid of size 30 km 20 km. Consider only diffusive movement in the x and y directions. Given a control representing the density of vaccine baits at location on week the corresponding susceptible ( if not previously exposed to rabies), infectious ( if able to transmit rabies), and immune ( if vaccinated) raccoon population densities satisfy the following equationsfor all with initial conditionsand no-flux boundary conditions

The parameter values and initial values are given in Table 2 which are taken from Neilan and Lenhart (2011).

Table 2

Description of parameters in model (2).

Parameters	Homogeneous	Heterogeneous	Description
Ω	30km×20km	same as homogeneous	Spatial domain
Ω_l	3km×2km	same as homogeneous	Spatial subdomain
S_0(x,y)	30/km2 for (x,y)∈Ω	0 for (x,y) on river	Initial susceptible raccoons
		10 for (x,y) in forest
		30 for (x,y) elsewhere
I_0(x,y)	30/km2 for (x,y)∈Ω_l	same as homogeneous	Initial infectious raccoons
	0/km2 for (x,y) elsewhere
R_0(x,y)	0/km2 for (x,y)∈Ω	same as homogeneous	Initial recovered raccoons
a_11(x,y)	0.5 km2/week for (x,y)∈Ω	0.01 for (x,y) on river	Diffusion coefficient
		0.20 for (x,y) in forest
		0.50 for (x,y) elsewhere
a_22(x,y)	0.5 km2/week for (x,y)∈Ω	0.20 for (x,y) in forest	Diffusion coefficient
		0.50 for (x,y) elsewhere
b(t)	0.096/week for 13≤t<28	same as homogeneous	Birth rate
	0/week otherwise
μ_1	0.026/week	same as homogeneous	Natural death rate
μ_2	0.490/week	same as homogeneous	Infectious death rate
β	0.03/(raccoons/km2·week)	same as homogeneous	Transmission rate
a	0.01/(raccoons/km2·week)	same as homogeneous	Vaccination rate
c	0.10 (raccoons/km)2/(vaccine)2	same as homogeneous	Balancing coefficient
T	20 weeks	same as homogeneous	Vaccine duration

Raccoons give birth during the spring of each year, March 20-June 21, a period of approximately 14 weeks (Clayton et al., 2010). Assuming a 50/50 sex rate within the population and that half the population are mature females, a reproductive rate of 1.34 year−1 is estimated (Coyne et al., 1989, Clayton et al., 2010). Dividing this yearly rate by the 14 weeks, we obtain week−1 for t within the birthing period. Let when t is not within the birthing period. For the simulations, assume the birthing period to be weeks 13 through 27. The constant year-long natural death rate, week−1, is calculated so that in absence of any disease or spatial spread, the susceptible population at and weeks are approximately equal (Clayton et al., 2010). Rabies-related death rate is estimated to be week−1 (Coyne et al., 1989, Clayton et al., 2010). The infection rate β is taken to be 0.03 (Russell et al., 2006, Smith et al., 2002, Clayton et al., 2010). The vaccine uptake rate a with units (vaccine · week)−1 is an indication of how successful the grounded baits are in vaccinating a raccoon. That is, to successfully vaccinate a raccoon, a bait must be first found and then eaten by a susceptible raccoon. This process can be inhibited by deterioration of the bait, human removal of the bait, or consumption of the bait by an animal other than a susceptible raccoon. Here (vaccine · week)−1 and how the parameter value influences the optimal control will be discussed later.

Neilan and Lenhart (2011) numerically approximated optimal vaccination strategies for (i) a homogeneous spatial domain with constant diffusion coefficients and a uniform initial susceptible population and (ii) a heterogeneous spatial domain with spatially dependent diffusion coefficients and a heterogeneous initial susceptible population. In heterogeneous case, movement within the forested area and across the river is inhibited and the initial susceptible population is assumed to be the largest in non-forested (urban) areas and absent on the river. See Table 2 for a list of all parameter values used in the homogeneous and heterogeneous examples. The boundary conditions imply that raccoons neither enter nor exit the domain. The set V of admissible controls consists of all measurable functions satisfying a.e. where is a large positive constant representing an upper bound on the density of baits placed at each location. Optimal control problem can be stated as follows: Find which minimizes the objective functionalsubject to system (2) and the initial and boundary conditions, where T is the number of weeks over which the control and observe population dynamics are applied.

Assume that the initial infection spreads for twenty-one weeks without intervention. Given the progression of the infection on week (Fig. 2), the optimal 20-week vaccination starting on week was computed. Fig. 3 display the results for the homogeneous and heterogeneous domains, respectively. Both schemes show control being applied at heaviest amounts initially and continuing with tapering amounts. Although, control is allowed to be applied for twenty weeks (), the simulations suggests very little or no vaccine is needed after ten weeks. Both vaccine strategies successfully eliminate rabies in the domain by week

Fig. 2

(a) In the absence of vaccination, infection starting in spreads as an expanding wave throughout the homogeneous spatial domain. (b) In the absence of vaccination, infection starting in spreads irregularly throughout the heterogeneous spatial domain. Figures were adopted from Neilan and Lenhart (2011).

Fig. 3

(a) For the homogeneous spatial domain, the optimal vaccination starting at week 21 is shown at weeks 21, 24 and 29. (b) For the heterogeneous spatial domain, the optimal vaccination starting at week 21 is shown at weeks 21, 25 and 29. Figures were adopted from Neilan and Lenhart (2011).

For the homogeneous case, the strategy is to immediately place an arc of vaccine bait in front of the infectious wave. The arc extends from the top to the bottom of the domain and its exact placement is determined by the vaccine uptake parameter In Fig. 3(a), and the arc is placed well in advance of the infectious wave. This placement allots ample time for susceptible raccoons to eat the baits and become immune. This immunity is sufficient to prevent rabies from spreading past the vaccine barrier. In additional simulations, it was found that by increasing the value of the optimal vaccine distribution remains in the shape of an arc but is placed closer to the infectious wave. In the corresponding heterogeneous case, the optimal vaccination is considerably curtailed, only forming a partial arc in front of the infectious wave and allowing the natural barriers to act in place of vaccination. Initially, the edges of forested areas and areas containing the river are fortified with a relatively small quantity of vaccine. With the partial vaccine arc and natural land features, the resulting population dynamics indicate that rabies does not cross to the right side of the river.

A spatiotemporal rabies model for skunk and bat interaction (Borchering et al., 2012)

The spread of rabies in areas such as Arizona and Texas with overlapping reservoir species (bats and skunks) is a unique problem. Bats are a major source of indigenously acquired human rabies infection in the United States and more than 2000 rabies-positive bats are collected annually. Cross-species transmission cases of rabies from bats to humans and other animals have been documented. Rabid skunks were diagnosed as infected with rabies virus of bat origin in Arizona.

Focusing on a geographic area of 300 km2 located in northeastern Texas, Borchering et al. (2012) propose a coupled system of nonlinear ordinary and partial differential equations to model the spatiotemporal dynamics of stripped skunks and eastern red bats interactions. Let and denote the numbers of susceptible, exposed, and infectious skunks, respectively, with the total number of skunks The bat population is divided into four groups, susceptible bats exposed bats infectious bats and recovered bats and the total number of bats is Logistic growth is assumed for both populations with appropriate birth rates ( and ) and carrying capacities ( and ). Skunks are susceptible to infection from skunks and bats. The term represents infected skunks produced per year resulting from contact between infected and susceptible skunks at a transmission rate Susceptible skunks progress into the exposed compartment after being inoculated with rabies virus due to the contact with infected skunks. The transmission function represents skunk infection resulting from contact with infected bats. The term represents the infection of susceptible bats by infected bats at a bat transmission rate After an average incubation period of exposed individuals move into the infected compartment. The incubation period for skunks is In the exposed compartments, individuals die from background mortality (terms and ). In the infected compartments, individuals die at a much higher rate that accounts for disease related mortality (terms and ). Recovered bat mortality is expressed by Diffusion terms ( and ) have been added to the infected compartments. The model of coupled ODEs/PDEs takes the form:

All parameter values are given in Table 3 and are taken from Borchering et al. (2012).

Table 3

Description of parameters in model (3).

Parameters	Value	Information	Description
K_s	20	Density 0.7−18.5 skunks/km^2	Skunk carrying capacity
β_s	2.5	Unknown	Skunk transmission rate
r_s	4	Litter size 3−9	Skunk birth rate
1/m_s	2.5	Lifespan 2−3 years	Skunk lifespan
1/m_rs	0.0274	Years (10 days)	Skunk disease progress rate
1/σ	0.164	Years (60 days)	Skunk incubation period
d_ss,d_es,d_is	10	km2 per year	Skunk diffusion coefficients
K_b	250	Density bats/km^2	Bat carrying capacity
β_b	0.12	Estimated	Bat transmission rate
r_b	0.4	Litter size 1−4	Bat birth rate
1/m_b	10	Lifespan 12 years	Bat lifespan
1/m_rb	0.0384	Years (14 days)	Bat disease progress rate
m_wb	1/10	Years (same as μ_b)	Bat mortality rate
1/σ_b	0.0384	Years (14 days)	Bat incubation period
1/ρ_b	0.5	Years	Bat recovery rate
γ	0.05	Unknown	Contact rate between skunks and bats
d_b	300	km2 per year	Bat diffusion coefficient

Bats and skunks account for the majority of rabies cases in northeastern Texas and the influence of bats on the spatial distribution and rabies dynamics of skunks are apparent when the map data of confirmed rabies cases are processed (Rabies Maps, 2011). The maps from 2003 to 2010 are altered to have a uniform size and orientation and the annual skunk and bat cases are plotted in a uniform format (Fig. 4). Numerical simulations using both the old model (skunk only) and the new model (skunks and bats) are given in Fig. 4(a). Gaussian distributions were used to instantiate the infected compartments with an approximation of the 2007 confirmed case data (Fig. 4(b)). The simulations indicate that the model with overlapping reservoir species more accurately reproduces the progression of rabies spread in northeastern Texas.

Fig. 4

(a) Simulations using both the old model (skunks only) and the new model (skunks and bats). (b) Gaussian distributions were used to instantiate the infected compartments with an approximation of the 2007 confirmed case data. Figures were adopted from Borchering et al. (2012).

Spatial models for rabies among dogs (Zhang et al., 2012)

All species of mammals are susceptible to rabies virus infection, but dogs remain the main carrier of rabies and are responsible for most of the human rabies deaths in China. To model the spatial spread of rabies among dogs, denote the total population density of dogs by and classify them into four subclasses: susceptible, exposed, infectious, and vaccinated classes, and their densities at time t and location are denoted by and respectively. The model is a dog-only subsystem of the one studied in Zhang et al. (2012), which is a reaction-diffusion SEIRS model of the following form:for where are the non-negative diffusion rates. All parameters are described in Table 4 are taken from Zhang et al. (2012).

Table 4

Description of parameters in model (4) (Zhang et al., 2012).

Parameters	Value	Unit	Description
A	3×10^6	year^−1	Dog birth population
λ	1	year^−1	Dog loss rate of immunity
i	16	year	Dog incubation period
σ	6	year^−1	1/i
γ	0.4	year^−1	Clinical outcome rate of exposed dogs
m	0.08	year^−1	Dog natural mortality rate
β_dd	1.58×10^−7	none	Dog-to-dog transmission rate
k	0.09	year^−1	Dog vaccination rate
μ	1	year^−1	Dog disease-related death rate
d_1	0.005	km·year^−1	Diffusion rate for the susceptible dogs
d_2	0.01	km·year^−1	Diffusion rate for the exposed dogs
d_3	0.01	km·year^−1	Diffusion rate for the infected dogs
d_4	0.005	km·year^−1	Diffusion rate for the vaccinated dogs

The dynamics of ODE version of model (4) have been studied in Zhang, Jin, Sun, Zhou, and Ruan (2011). It is known that there exists a disease-free equilibrium

If the basic reproduction number (Zhang et al., 2011)then there is a unique endemic equilibrium

The traveling wave solutions of system (4) is rewritten in term of a coordinate frame to the right with speed c; i.e., with . The traveling waves satisfy the boundary conditions:and

We now give some numerical results about the existence of traveling waves. The initial data are:

Through drawing two-dimension figures of the population number in every subclasses of population in one-dimension space (Fig. 5), it can be seen that with the movement of dogs there exist traveling waves in every subclasses of dogs. Thus, the dispersal of dogs induces the epidemic waves of rabies among the dog population.

Fig. 5

Traveling wave solutions of model (4) with parameters given in Table 4. The solutions are plotted when t = 2, 8, 16, 24, 32, 40. Figures were adopted from Zhang et al. (2012).

Discussion

We briefly reviewed some reaction-diffusion models describing the spatial spread of rabies among animals. More specifically, we introduced the susceptible-exposed-infectious models for the spatial transmission of rabies among foxes (Murray et al., 1986), the spatiotemporal epidemic model for rabies among raccoons (Neilan & Lenhart, 2011), the diffusive rabies model for skunk and bat interactions (Borchering et al., 2012), and the reaction-diffusion model for rabies among dogs (Zhang et al., 2012). Estimated parameter values and numerical simulations on the spatiotemporal dynamics of these models from these papers were presented.

The numerical simulations of traveling waves in these models indicate that the spatial spread of rabies is caused by the dispersal of the host animals. For the rabies model (1) among foxes, Yachi, Kawasaki, Shigesada, and Teramoto (1989) proved the existence of traveling wave solutions in the reaction-diffusion equations. It will be useful and interesting to study the dynamics of these rabies models (Wang, 2014), such as the threshold dynamics, the stability of steady states, and the existence of traveling solutions.

The diffusion model provides a useful framework to evaluate some spatially related control measures such as the possibility of stopping the spread of the disease by creating a rabies ‘break’ ahead of the front through vaccination to reduce the population to a level below the threshold for an epidemic to exist. It should be mentioned that there are other approaches to model the spatial spread of rabies among animals, such as multi-patch models (Russell et al., 2006, Dimitrov et al., 2008, Chen et al., 2015) and stochastic spatial models (Smith et al., 2002), in particular when the spread of rabies on heterogeneous landscapes is concerned.