SPATIAL SPREAD OF RABIES IN CHINA ∗

The reported data on human rabies in Mainland China in the past 10 years demonstrate that the number of counties with human rabies cases has increased dramatically and the disease has spread from southern provinces to northern provinces. The aim of this paper is to study how the movement of dogs impacts the spatial spread of rabies and to find effective strategies for controlling rabies. We firstly construct a reaction-diffusion model for dogs and humans and obtain the minimal wave speed. Then we consider the existence of traveling waves and the influences of parameters on human rabies cases and the minimal wave speed by numerical simulations. Lastly, we discuss the local stability of the endemic equilibrium. Our analysis indicates that dog movement leads to the traveling wave of dog and human rabies and has a large influence on the minimal wave speed. Thus restricting the movement of dogs, in particular the exposed and infected dogs, is an effective way to prevent the spatial spread of rabies.


Introduction
Rabies roots in its enzootic environment (animal host) and is often transmitted through the bite of a rabid animal.It causes severe threats to public health worldwide [19].Though rabies has a nearly 100% mortality, it is a preventable viral disease of mammals [4].Human rabies can be effectively treated if the exposure is recognized before the symptoms develop.More detailed information about rabies can be founded in [19,3,17,18,4].Numerous wildlife species are reservoirs of rabies virus and sources of human rabies.In Mainland China, however, domestic dogs, especially in the countryside, play a main role in human rabies.In fact, 85%-95% of human rabies cases are due to dog bites and 50%-70% of human cases are reported in the countryside [10].In rural areas, nearly every family owns several dogs and most of them are free-roaming.In recent years, rural communities and The figures are derived from [16].
areas invaded by rabies are gradually enlarged.Fig. 1(a) shows that there were only 120 counties reporting human rabies cases in 1999.However, the number of infected counties in 2008 was nearly seven times as that in 1999.Moreover, we can see the tendency of rabies diffusion from Fig. 1(b).In 1996, the 160 infected cases were sporadically distributed in the south provinces.As time goes on, the situation of rabies in south provinces has become more serious and it has spread to central and north provinces.
Though the situation of rabies in Mainland China is so severe, there are very few studies on modeling the transmission dynamics of rabies, in particular the spatial spread of the disease, in Mainland China [9,20,21].Meng et al. [9] adopted a statistical method to study the rabies virus.Zhang et al. [20,21] constructed SEIRS models to study the transmission of rabies in Mainland China from 1996 to 2010.Hou et al. [7] proposed a deterministic model for the dog-human transmission od rabies, taking into account both domestic and stray dogs, and applied the model to simulate the reported human cases in Guangdong Province, China.Now we want to explore the factors behind the spatial spread of rabies in Mainland China.There are some studies on investigating animals rabies and using partial differential equations by adding diffusion terms to ordinary differential equations [8,11,6,2,15,14].Kallen [8] studied rabies transmission in fox population by differential equations with diffusion.He used the deterministic model to simulate rabies epizootic in foxes crossing continental Europe and proved the existence of traveling waves.Murray et al. [11] 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 mean to control the spatial spread of the disease.Artois et al. [2] employed a one-dimensional discrete deterministic model to analyze the evolution of a given population which is infected by rabies.Ou and Wu [14] considered the spatio-temporal patterns of rabies spread involving structured fox populations and delayed nonlocal interaction on the disease spread.
In this paper, we add diffusion to the dog population in a SEIRS model considered by Zhang et al. [20] to obtain a reaction-diffusion equation model.For reaction-diffusion equations, the asymptotic speed and traveling waves are attract-ing more attention and usually need to be investigated and discussed in order to understand the dynamics.The concept of asymptotic speeds of spread was first introduced in [1] for reaction-diffusion equations.There is an intuitive interpretation for the spreading speed c * in a spatial epidemic model: if one runs at a speed c > c * , then one will leave the epidemic behind; whereas if one runs at a speed c < c * , then one will eventually be surrounded by the epidemic [5].Following the method adopted by [12,13], we determine the minimum wave speed connecting the disease-free equilibrium to the endemic equilibrium and illustrate the existence of traveling waves by numerical simulations.It was shown that the endemic equilibrium of the ODE model is stable [20].Here we want to study the diffusive model by linear analysis to see whether the endemic equilibrium is still stable and whether spatial patterns can appear.Moreover, we carry out some sensitivity analysis on parameters to determine factors that affect the spatial spread of rabies in Mainland China.
The article is organized as follows.In Sections 2, we present the model without and with diffusion respectively.In Section 3, for the model with diffusion we calculate the minimum wave speed.The existence of traveling waves, sensitivity analysis of parameters about the human rabies cases, and the minimal wave speed c * are studied numerically in section 4. In section 5, we discuss the stability of the endemic equilibrium by considering the real parts of the eigenvalues in term of different parameters.A brief discussion is given in section 6.

Model for rabies
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 Mainland China.So the model will include the populations of humans and dogs.We denote the total number of dogs and humans by N d (t), N h (t), and classify each of them into four subclasses: susceptible, exposed, infectious and vaccinated, with dog population denoted by S d (t), E d (t), I d (t), and R d (t), and human population denoted by S h (t), E h (t), I h (t), and R h (t) respectively.The model we study is presented in Eqs.(1) [20].The parameters are described in Table 1.
We focus on a SEIRS model governing the temporal evolution of rabies as follows.
The dynamics of model (2.1) have been studied in [20].The main results can be summarized as follows.
The disease-free equilibrium of system (2.1) is The basic reproduction number is [20] This is the threshold value to determine whether rabies can invade a susceptible population and is important for designing control strategies.
If R 0 > 1, we can derive the unique endemic equilibrium: where Here, In [20], we show that if R 0 < 1, then the disease-free equilibrium E 0 is globally asymptotically stable.If R 0 > 1, E 0 becomes unstable, and the endemic equilibrium E * of system (2.1) is locally asymptotically stable.
In this paper, we consider spatial spread of rabies in China.Let S d (x, t), E d (x, t), I d (x, t), and R d (x, t) to denote the density of susceptible, exposed, infectious, and recovered dogs at location x and time t, respectively.The corresponding human sub-populations are denoted by S h (x, t), E h (x, t), I h (x, t), and R h (x, t), respectively.Suppose that the spatial diffusion of rabies is a result of the movement of dogs.So we add the diffusion of dogs to system (2.1) and obtain the following reaction-diffusion equation model: for t > 0, x ∈ Ω with null Neumann conditions: where, ∆ = ∂ 2 ∂ 2 x and d 1 , d 2 , d 3 , d 4 are the non-negative diffusion rates whose values are given in Table 1.
Since the last four equations are independent of the first four equations and only the first four equations have diffusion terms, we only need to study the subsystem: (2.4) In order to investigate traveling waves, system (2.4) is rewritten in term of a coordinate frame to the right with speed c.Let z = x − ct, we can rewrite system (2.4) as the following form: (2.5) We have the corresponding first-order ordinary differential equations with respect to the variable z of above system: , , , which is the system we will study in next section.

The minimum wave speed
The purpose of this section is to determine the minimum wave speed c * connecting the disease-free equilibrium E 0 to the endemic equilibrium E * .The solutions corresponding to the minimum wave speed describe the observed epidemic waves and are called traveling waves.The traveling waves, if they exist, satisfy the following conditions: = (S 0 d , 0, 0, 0, 0, 0, R 0 d , 0) and ).As for system (2.6), we need to consider the eigenvalues of its Jacobian matrix at the equilibrium E 0 to find the minimum wave speed for traveling waves.The disease-free equilibrium E 0 has two zeros corresponding to dogs, which must not be negative.Considering this constraint, we impose that solutions of the linear system must not be oscillate, i.e., the eigenvalues corresponding to E 0 must be real values.The Jacobian matrix at the equilibrium E 0 is as follows: .
As for this matrix, the eigenvalues are the roots of the polynomial equation: where Firstly, we investigate the roots of f (x).Let By analysing the coefficients of f 1 (x), we know that f 1 (x) has four intersection points with the x axis.It is easily obtained that f 2 (x) is over the x axis from f 2 (x) = λk > 0.Moreover, we obtain f 1 (0) = (m + λ)(m + k) > λk = f 2 (x).So we show their relations in Fig. 2 and see that f 1 (x) and f 2 (x) have four intersection points, meaning that f (x) has four real roots.
Next, we study g(x) which is more complex.Let Because of its complexity, we consider the special situation that d 2 = d 3 .When d 2 = d 3 , g 1 (x) has three extreme points So, we know that x 1 is the maximum point and calculate the maximum value If g 1 (x 1 ) > g 2 (x) = β dd S 0 d σr, there are four intersection points for g 1 (x) and g 2 (x).If g 1 (x 1 ) < g 2 (x) = β dd S 0 d σr, there only appear two intersection points for g 1 (x) x 1 Figure 3.A general view of g(x).and g 2 (x).So, g 1 (x 1 ) = g 2 (x) = β dd S 0 d σr is the condition we need to obtain the minimum wave speed c * .We describe the relation between g 1 (x) and g 2 (x) as in Fig. 3.
Therefore, when Putting the values of parameters in Table 1 into the above formula, we obtain c * = 0.1792.

Traveling wave and sensitivity analysis
In this section, we give some numerical results about the existence of traveling waves, sensitivity analysis of parameters about the human rabies cases and the minimum wave speed.The initial data we take in numerical simulations of system (2.2) are as follows: From the two-dimension figures on the population number in every subclasses in the one-dimension space (Figs. 4 and 5), we can see that with the movement of dogs, not only there really exist traveling waves in every subclasses of dogs but also there appears epidemic waves in the infected human population.with parameters given in Table 1.The solutions are plotted when t=2, 8, 16, 24, 32, 40.Moreover, from Fig. 5 we can see that the diffusion of dogs not only spreads the disease to other adjacent areas but also makes the rabies situation more serious, which agrees with the fact that not only rabies has spread from the south to the central and north areas but also the situation of rabies in the south is getting more serious.
We have known that the traveling waves exist in the dog population.We now explore the influencing factors for the traveling waves.Taking this into account, Fig. 6 intuitively illustrates how the parameters impact the transmission and diffusion of rabies.We can visually see that the influences of d 1 , d 2 , d 3 , d 4 are smaller than that of A and k.Moreover, d 1 , d 2 , d 3 , d 4 only affect the speed of traveling waves, while other parameters not only affect the speed of waves but also change the magnitude of waves.This illustrates that strengthening quarantine during the dog movement only can control the diffusion spread of rabies to other areas but cannot eliminate rabies.In addition, in Fig. 6(a) we can observe that the bottle-green curve is convergent, which is because R 0 < 1 when A = 2 × 10 6 . (a)

Steady states and stability analysis
In [20], we presented the steady states and discussed their stability without diffusion.Now we want to know whether the endemic equilibrium E * is still stable when the system is diffusive, that is, if there appears any spatial pattern.Since the last four equations are independent of the first four equations, we only need to study the point . Firstly, we linearize the subsystem about this point.Set s(x, t) = s 0 e αt+inx , e(x, t) = e 0 e αt+inx , i(x, t) = i 0 e αt+inx , (5.2) r(x, t) = r 0 e αt+inx , where s 0 , e 0 , i 0 and r 0 are constants, e inx is periodic and bounded (|e inx | = 1) and n is the wavenumber that indicates the wavelength of the emergent spatial pattern.The sign of the real part of α is crucially important to determine the stability of the steady state.
Substituting (5.2) into (2.4) and ignoring nonlinear terms, we obtain the following equations where I is the 4 × 4 identity matrix and −A =     (a) where Using the Routh-Hurwitz criterion, we know that all roots of (5.3) have negative real parts if and only if  1 and n = 5, all roots of (5.3) have negative real parts, that is, the endemic equilibrium is stable.
We also have the largest real parts of roots in terms of ) is stable, that is, E * is stable and spatial patterns cannot appear.Moreover, in Fig. 10, we can intuitively see that the endemic equilibrium with small perturbations is still stable.

Discussion
The reported data of human rabies in Mainland China demonstrate that the number of counties in 2008 with reported cases was nearly seven times as that in 1999.It not only becomes more serious in south provinces but also has spread to the north and northwest provinces.Many wildlife species are reservoirs of rabies virus and sources of human rabies.In Mainland China, however, dogs are the main source of human rabies.
To investigate how the movement of dogs leads to and impacts the spatial spread speed of rabies, we constructed a reaction-diffusion equation model to describe the interactions between dogs and humans, where dogs can disperse.We calculated the minimum wave speed and considered the existence of traveling wave by numerical simulations.Moreover, we discussed the stability of the endemic equilibrium with diffusion to illustrate that the spatial pattern cannot appear.The key is to check the magnitude of influence of parameters on human rabies cases and the minimal wave speed c * .From the formula of c * , we can see that only d 2 , d 3 (= d 2 ) can influence c * , but d 1 , d 4 do not.Therefore, only the movement of exposed and infected dogs affects the spatial spread of rabies.From Fig. 7(b), we know that the influence of d 2 , d 3 on c * is larger which illustrates that quarantine of dogs during movement is important to control the spatial spread of rabies.
The existence of traveling waves in the model indicates that the spatial spread of rabies in both dogs and humans is caused by the dispersal of dogs.Thus, management of dogs and dog movement is important to reduce and prevent such spread of rabies.To study the long distance inter-province spread of rabies, multi-patch models might be helpful and we leave this for future consideration.

Figure 1 .
Figure 1.(a) The data on the numbers of counties reported rabies cases in Mainland China from 1999 to 2008.The figure is derived from [10].(b) Geographic distributions of human rabies in Mainland China in 1996, 2000 and 2006.The figures are derived from [16].

Figure 5 .
Figure 5.The density of infected human I h in the PDE model (2.2) with the diffusion condition on the dog population.The solutions are plotted when t=8, 16, 24, 32.

5 Figure 6 .
Figure 6.The waves of human rabies cases I h in term of x under different parameter values.The solutions are plotted when t = 8, 16, 24, 32.(a) with different A; (b) with different k; (c) with different d 1 , d 2 , d 3 , d 4 .

Figure 7 .
Figure 7.The minimum wave speed c * in term of different parameters.

Figure 8 .
Figure 8.The largest real parts of roots in term of parameters d 1 , d 2 , d 3 , d 4 and n.
h 1 year −1 Human 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 where Figure10.The temporal and spatial solution of I d (x, t), R d (x, t), I h (x, t) and R h (x, t) for the parameters given in Table1.The initial values are the endemic equilibrium values with small perturbations.So, under the parameter values in Table a 1 = 9.4585 > 0, a 2 = 14.882 > 0, a 3 = 6.6762 > 0, a 4 = 1.4153 > 0, a 1 a 2 − a 3 = 134.0885>0,(a 1 a 2 − a 3 )a 3 − a21 a 4 = 768.5886>0.h