key: cord-0783443-vxcz2ogt
authors: Zhang, Juan; Jin, Zhen; Sun, Gui-Quan; Sun, Xiang-Dong; Ruan, Shigui
title: Modeling Seasonal Rabies Epidemics in China
date: 2012-03-01
journal: Bull Math Biol
DOI: 10.1007/s11538-012-9720-6
sha: 91096eede9d3809948df0fa851951341b12998c6
doc_id: 783443
cord_uid: vxcz2ogt

Human rabies, an infection of the nervous system, is a major public-health problem in China. In the last 60 years (1950–2010) there had been 124,255 reported human rabies cases, an average of 2,037 cases per year. However, the factors and mechanisms behind the persistence and prevalence of human rabies have not become well understood. The monthly data of human rabies cases reported by the Chinese Ministry of Health exhibits a periodic pattern on an annual base. The cases in the summer and autumn are significantly higher than in the spring and winter. Based on this observation, we propose a susceptible, exposed, infectious, and recovered (SEIRS) model with periodic transmission rates to investigate the seasonal rabies epidemics. We evaluate the basic reproduction number R (0), analyze the dynamical behavior of the model, and use the model to simulate the monthly data of human rabies cases reported by the Chinese Ministry of Health. We also carry out some sensitivity analysis of the basic reproduction number R (0) in terms of various model parameters. Moreover, we demonstrate that it is more reasonable to regard R (0) rather than the average basic reproduction number [Formula: see text] or the basic reproduction number [Formula: see text] of the corresponding autonomous system as a threshold for the disease. Finally, our studies show that human rabies in China can be controlled by reducing the birth rate of dogs, increasing the immunization rate of dogs, enhancing public education and awareness about rabies, and strengthening supervision of pupils and children in the summer and autumn.

Keywords Rabies · SEIRS model · Basic reproduction number · Periodic solution · Vaccination

Rabies, a fatal disease of humans and many other mammals, is caused by a virus which is associated with the bite and virus-containing saliva of an infected host (CDC 2010a) . Rabies may affect all mammals, including livestock and pets. In most African and Asian countries dogs continue to be the main hosts and are responsible for most of the human rabies deaths (WHO 2010a) .

Infective dogs can bite other dogs and humans and spread the rabies virus. After entering the body, the rabies virus travels quickly along the neural pathways to the noncentral nervous system, from there the virus further spreads to other organs and causes morbidity by intruding many tissues. Infected individuals can experience incubation before showing symptoms. The incubation period of the disease is usually a few months, but can be as long as years, depending on the distance the virus must travel to reach the central nervous system. If the bitten location is near the head, the incubation is relatively shorter. The first symptoms of rabies may be very similar to those of the flu expressing fever or headache, progressing to symptoms of cerebral dysfunction, anxiety, confusion, fearing light and water within days. With the disease exacerbating, the infected individual may experience delirium, abnormal behavior, hallucinations, and insomnia (CDC 2010b) . Once the rabies virus reaches the central nervous system and symptoms of the disease develop, the course of the disease is less than 10 days (CDC 2010b) and the mortality rate reaches up to 100%.

Rabies transmission in dogs can be prevented by vaccination. Treatment of humans after exposure, known as post-exposure prophylaxis (PEP), is highly successful in preventing the disease if administered promptly. Although it is a vaccinepreventable disease, rabies still remains a neglected, untreatable public-health problem in many countries in Asia and Africa where 95% of human deaths occur (WHO 2010b) .

Recently, mathematical models have been used to study the rabies epidemics in dogs and the transmission dynamics of rabies from dogs to humans. For example, Hampson et al. (2007) observed rabies epidemics cycles with a period of 3-6 years in dog populations in Africa and built a susceptible, exposed, infectious, and vaccinated model with an intervention response variable to show significant synchrony. Zinsstag et al. (2009) classified both dog and human populations in susceptible, exposed, infectious, and immunized classes and proposed a model for dog-human transmission dynamics and economics of rabies control in an African city.

Human rabies is a major public-health problem in China. Since 1950, human rabies has been classified as a class II infectious disease in the National Stationary Notifiable Communicable Diseases and the annual data of human rabies have been archived by the Chinese Center for Disease Control and Prevention. In the last 60 years there had been 124,255 human rabies cases reported by the Chinese Ministry of Health (MOHC 2009 (MOHC , 2011 , an average of 2,037 cases per year. However, the factors and mechanisms behind the persistence and prevalence of human rabies have not been well understood. Most recently, Zhang et al. (2011) developed a deterministic model to study the transmission dynamics of rabies in China. The model, consisted of susceptible, exposed, infectious, and recovered subpopulations of both dogs and humans and described by eight ordinary differential equations, was used to simulate the human rabies data from 1996 to 2010 reported by the Chinese Ministry of Health. It was shown that reducing dog birth rate and increasing dog immunization coverage rate are the most effective methods for controlling rabies in China and large scale culling of susceptible dogs can be replaced by immunization of them. Hou et al. (2012) proposed a similar dog-human interaction model considering both domestic and wild dogs and used the model to simulate the rabies data from Guangdong Province, China. It was shown that the quantity of stray dogs also plays an important role in the transmission of rabies.

After the outbreaks of Severe Acute Respiratory Syndromes (SARS) in 2003, the Chinese Ministry of Health started to publish reported cases about the National Stationary Notifiable Communicable Diseases every month. We observe that the monthly data of human rabies cases reported by the Chinese Ministry of Health since January 2004 exhibit a periodic pattern on an annual base. The cases in the summer and autumn are significantly higher than in the spring and winter (MOHC 2009 (MOHC , 2011 . Song et al. (2009) also reported that the main seasons for rabies epidemics in China are summer and fall. Moreover, the infected areas are mainly distributed in the south provinces such as Sichuan, Hunan, Guangxi, Guangdong, Anhui, Fujian (Song et al. 2009 ), which demonstrates that rabies transmission depends on the weather.

It is well-known that many diseases exhibit seasonal fluctuations, such as whooping cough, measles, influenza, polio, chickenpox, mumps, etc. (Bjornstad et al. 2002; Dowell 2001; London and Yorke 1973) . Seasonally effective contact rate (Dushoff et al. 2004; Schwartz 1992; Schwartz and Smith 1983; Smith 1983) , periodic changing in the birth rate (Ma and Ma 2006) and vaccination program (Earn et al. 2000) are often regarded as sources of periodicity. In this paper, we take periodic transmission rate into account based on the following facts. (1) In the summer and fall, people wear light clothing and are lack of protection for the bites or scratches of dogs. Also, in summer and fall people, in particular farmers, have more frequent outdoor activities which increase the chance of human-dog interaction. (2) In these seasons, dogs are more maniacal and apt to attack each other and humans. (3) From July to September schools are closed for summer vacations and children are out of supervision and enjoy tantalizing dogs. In fact, it was reported (Song et al. 2009 ) that 25.7% of human rabies cases in China are students and unattended children. (4) In addition, temperature may be related to the fluctuation of diseases. Under high temperature in summer, rabies virus can survive easily and its infectivity is stronger.

The purpose of this paper is to propose a susceptible, exposed, infectious, and recovered (SEIRS) model with periodic transmission rates to investigate the seasonal rabies epidemics. We will evaluate the basic reproduction number R 0 , analyze the dynamical behavior of the model, and use the model to simulate the monthly data of human rabies cases reported by the Chinese Ministry of Health from January 2004. We will also carry out some sensitivity analysis of the basic reproduction number R 0 in terms of various model parameters. Moreover, we will show that it is more reasonable to regard R 0 rather than the average basic reproduction numberR 0 or the basic reproduction numberR 0 of the corresponding autonomous system as a threshold for the disease. Finally, we will explore some effective control measures for the rabies epidemics in China.

The article is organized as follows. In Sect. 2, we introduce the model and present the expression of the basic reproduction number. Then we study the global asymptotic stability of the disease-free equilibrium and the existence of positive periodic solutions. Simulations of the model and sensitivity analysis of the basic reproduction number are performed in Sect. 3. In Sect. 4, we give a brief discussion. The detailed calculation of the basic reproduction number using the spectral radius of the operator is presented in the Appendix.

We denote the total numbers of dogs and humans by N(t) and N 1 (t), respectively, and classify each of them into four subclasses: susceptible, exposed, infectious and recovered, with the numbers of dogs denoted by S(t), E(t), I (t), and R(t), and human sizes denoted by S 1 (t), E 1 (t), I 1 (t), and R 1 (t), respectively. The transmission dynamics associated with these subpopulations are illustrated in Fig. 1 .

The transmission rate between S(t) and I (t) is β(t), the transmission rate between S 1 (t) and I (t) is β 1 (t), and humans do not spread rabies to each other. We can write transmission rate β(t) in the general form β(t) = λ 0 β (N )β /N , where N is the total number of dogs, β (N ) is the number of dogs that a susceptible dog comes across per unit time, β is the probability of getting bitten after interacting with the susceptible dog, and λ 0 is the probability of being infected after bitten for the susceptible dog. We can express β 1 (t) similarly. As discussed in the Introduction, in the summer and fall there are more frequent interactions among dogs and between dogs and humans and these coefficients are more likely to change as season changes. Thus we use the periodic functions β(t) = a[1 + b sin( π 6 t + 5.5)] and β 1 (t) = a 1 [1 + b 1 sin( π 6 t + 5.5)] proposed by Schenzle (1984) to describe the transmission rates among dogs and from dogs to humans, where a and a 1 are the baseline contact rates and b and b 1 are the magnitudes of forcing.

The birth numbers of dogs and humans per unit time are constant. Vaccination is often applied to seemingly healthy dogs (S(t) and E(t)) and people bitten by dogs (E 1 (t)). Particularly, we need to interpret that k 1 and k are the products of the vaccination coverage rate and the vaccination effective rate. However, there is a protection period for rabies vaccine. Thus, we import loss rates of immunity λ and λ 1 . Because not all the exposeds will develop clinical outbreak, clinical outcome rates γ and γ 1 are presented. Natural death rates are m and m 1 , and disease-related death rates are μ and μ 1 , respectively. The model is a system of ordinary differential equations:

where all parameters are positive, the interpretations and values of parameters are described in Table 1 , β(t) = a[1 + b sin( π 6 t + 5.5)] and β 1 (t) = a 1 [1 + b 1 sin( π 6 t + 5.5)].

Notice that from the equations in model (1), we have

(2) Let X = {(S, E, I, R, S 1 , E 1 , I 1 , R 1 )|S, E, I, R, S 1 , E 1 , I 1 , R 1 ≥ 0, 0 < S + E + I + R ≤ A m , 0 < S 1 + E 1 + I 1 + R 1 ≤ B m 1 }. 

The region X is positively invariant with respect to system (1).

It is easy to see that system (1) has one disease-free equilibrium

We can evaluate the basic reproduction number R 0 for system (1) following the definition of Bacaer and Guernaoui (2006) and the general calculation procedure in Wang and Zhao (2008) , which is defined as z 0 such that g(z 0 ) = 1, where

The detailed computations are given in the Appendix. Before analyzing the disease-free equilibrium, we make a claim.

For an arbitrary positive number θ , there is t 2 > 0 such that for all t > t 2 , S ≤Ŝ + θ .

Proof From the last equation of system (1), we have

Thus, for a positive number

The disease-free equilibrium P 0 is globally asymptotically stable when Wang and Zhao (2008) . We can choose θ > 0 small enough such that ρ(Φ F −V +M θ (ω)) < 1, where

Considering the region X and using Lemma 2.2, we know that S 1 (t) ≤Ŝ 1 = B m 1 and

Consider the following comparison system:

that is,

By Lemma 2.1 in Zhang and Zhao (2007) , it follows that there exists a positive

Therefore, we have h(t) → 0 as t → ∞, which implies that the zero solution of system (4) is globally asymptotically stable. Applying the comparison principle (Smith and Waltman 1995) , we know that for system (1), E(t) → 0, I (t) → 0, E 1 (t) → 0 and I 1 (t) → 0 as t → ∞. By the theory of asymptotic autonomous systems (Thieme 1992) , it is also known that S(t) →Ŝ, R(t) →R, S 1 (t) →Ŝ 1 and R 1 (t) → 0 as t → ∞. So P 0 is globally attractive when R 0 < 1. It follows that P 0 is globally asymptotically stable when R 0 < 1.

Define

Let P : X → X be the Poincaré map associated with system (1), i.e.,

where ω is the period. Applying the fundamental existence-uniqueness theorem (Perko 2000) , we know that u(t, x 0 ) is the unique solution of system (1) with u(0, x 0 ) = x 0 . From Theorem 2.1, we know that X is positively invariant and P is point dissipative.

When R 0 > 1, then there exists a δ > 0 such that when Wang and Zhao (2008) .

If not, then

Without loss of generality, we assume that d(P m (S 0 , E 0 , I 0 , R 0 , S 0 1 , E 0 1 , I 0 1 , R 0 1 ), P 0 ) < δ for all m ≥ 0. By the continuity of the solutions with respect to the initial values, we obtain

For any t ≥ 0, let t = mω + t 1 , where t 1 ∈ [0, ω] and m = [ t ω ], which is the greatest integer less than or equal to t ω . Then we have

Next we consider the linear system

Once again by Lemma 2.1 in Zhang and Zhao (2007) , it follows that there exists a positive ω-periodic functionĝ(t) such that g(t) = e ptĝ (t) is a solution of system (7)

Applying the comparison principle (Smith and Waltman 1995) , we know that when E(0) > 0, I (0) > 0, E 1 (0) > 0 and I 1 (0) > 0,

This is a contradiction. The proof of the lemma is complete.

Theorem 2.5 System (1) has at least one positive periodic solution.

Proof We first prove that {P m } m≥0 is uniformly persistent with respect to (X 0 , ∂X 0 ). First of all, we explain that X 0 and ∂X 0 are positively invariant. In fact, for any (S 0 , E 0 , I 0 , R 0 , S 0 1 , E 0 1 , I 0 1 , R 0 1 ) ∈ X 0 , solving the equations of system (1), we derive that

and

Similarly, S 1 (t) > 0, R 1 (t) > 0, E 1 (t) > 0 and I 1 (t) > 0. So, X 0 is positively invariant. Clearly, ∂X 0 is relatively closed in X. Set

It is easy to show that

Note that

we only need to prove that

That is, for any (S 0 , E 0 , I 0 , R 0 , S 0 1 , E 0 1 , I 0 1 , R 0 1 ) ∈ ∂X 0 , we have

If there exists an m 1 ≥ 0 such that

by replacing the initial time 0 with m 1 ω and following the processes as in (8)-(11), it can be seen that

So, the equality (12) holds, which implies that E 0 is the only fixed point of P and acyclic in ∂X 0 .

Moreover, Lemma 2.4 implies that P 0 = (Ŝ, 0, 0,R,Ŝ 1 , 0, 0, 0) is an isolated invariant set in X and W S (P 0 ) ∩ X 0 = ∅. By the acyclicity theorem on uniform persistence for maps (Theorem 1.3.1 and Remark 1.3.1 in Zhao 2003) , it follows that P is uniformly persistent with respect to (X 0 , ∂X 0 ). Now Theorem 1.3.6 in Zhao (2003) implies that P has a fixed point

From the first equation of system (1) we have

The periodicity of S * (t) implies S * (t) > 0 for all t > 0. Following the processes as in inequalities (8)

is a positive ω-periodic solution of system (1).

In this section, we first use model (1) to simulate the reported human rabies data of China from January 2004 to December 2010, predict the trend of the disease and seek for some control and prevention measures. The data, concerning human rabies from 2004 to 2010, are obtained mainly from epidemiologic bulletins published by the Chinese Ministry of Health (MOHC 2011). We need to estimate the parameters of model (1), most of which can be obtained from the literature or assumed on the basis of common sense. However, we have to estimate β(t), β 1 (t) and A by using the least-square fitting of I 1 (t i ) through discretizing the ordinary differential system (1) as follows:

The least-square fitting is to minimize the objective function

which is implemented by the instruction lsqnonlin, a part of the optimization toolbox in MATLAB. The values of parameters are listed in Table 1 . We obtain the annual number of human population using the annual birth and death data from the National Bureau of Statistics of China (NBSC 2009). Then we calculate the average and divide it by 12 to derive the monthly human birth population B = 1,340,000. We need the initial values to perform the numerical simulations of the model. The number of the initial susceptible human population at the end of 2003, S 1 (0), is obtained from the China Statistical Yearbook and the number of the initial infective humans I 1 (0) is from epidemiological bulletins published by the Chinese Ministry of Health. However, the numbers of the initial exposed humans E 1 (0) and the recovered humans R 1 (0) cannot be obtained. We derive E 1 (0) reversely by the parameter γ and R 1 (0) is estimated roughly. Regarding the initial values for dogs, we only know that there are about 75 millions dogs in 2009 from online news. So, S(0), E(0), I (0), and R(0) are calculated reversely by the corresponding parameters r 1 , k 1 and data fitting. The numerical simulation of the model on the number of human rabies cases is shown in Fig. 2 . We observe that the data of 2005, 2008 and 2009 are slightly different from the solution as observed in Zhang et al. (2011) . We think this is because of large scale culling of dogs in these years. However, culling of dogs is not considered in model (1).

Moreover, with these parameter values, we can roughly estimate that the basic reproduction number R 0 = 1.03 under the current circumstances in China. From Fig. 3 , we can see that when R 0 < 1, the number of infected humans I 1 (t) tends to 0. On the contrary, when R 0 > 1, I 1 (t) tends to a stable periodic solution.

We can also predict the general tendency of the epidemic in a long term according to the current situation, which is presented in Fig. 4 . From these figures we can see that the epidemic of rabies can be relieved in a short time, but cannot be eradicated with the current prevention and control measures. Fig. 3 The tendency of the human rabies infectious cases I 1 (t) in a long time with different values of R 0 . When A = 220,000 (lower curve) and 300,000 (upper curve), and the values of other parameters in Table 1 do not change, R 0 = 0.97 and 1.32, respectively Fig. 4 The tendency of the human rabies infectious cases I 1 (t) in short and long times

The initial conditions adopted in model fitting are mostly assumed and backextrapolated by parameters. So it is necessary to study the influence of initial conditions on the rabies epidemics which are showed in Figs. 5 and 6. From Figs. 5 and 6, we can see that the initial value S(0) has a stronger influence on I 1 (t) and other initial conditions have little or almost no effect on I 1 (t). It implies that the increasing number of dogs is really an important factor for the prevalence and persistence of rabies in China.

Finally, we perform some sensitivity analysis to determine the influence of parameters A, k, γ and a on R 0 . From Fig. 7(a) , it is obvious that when A is less than 226,920, R 0 can be less than 1. However, the annual birth population of dogs can achieve 400,000 or more in China. This indicates that human rabies in China cannot be eradicated if the birth number of dogs cannot be controlled under 2 million.

Post-exposure prophylaxis (PEP) is used for most situations for human rabies. In model (1), it is embodied in the terms k and k 1 and it can affect γ and γ 1 . We observe that R 0 is a concave function of k from Fig. 7(b) . So k has an obvious effect on R 0 . We also know that immunization is an effective measure to control rabies. Next, we consider the effect of γ on R 0 , which is depicted in Fig. 7(c) . We can observe that it is linear in γ . Most people, especially in the rural and remote areas, have little knowledge about rabies and even do not know what to do after being bitten by dogs. Song et al. (2009) reported that 66.3% of rabies victims did not seek medical services at all and 27.6% of the cases received inadequate PEP. Although the effect of γ on R 0 is less than k, we can enhance people the awareness and knowledge about rabies and the emergency measure and treatment after they are bitten and scratched by dogs to decrease the rate of clinical outbreak of rabies. Now we discuss how a effects R 0 in Fig. 7(d) . Although they are linear, a is very small and a slight change of a can lead to large variations of R 0 . Since β(t) = a[1 + b sin( π 6 t + 5.5)] = λ 0 β (N )β /N , we can manage a by controlling β (N ), i.e., the number of dogs a susceptible dog runs into per unit time, which is to strengthen the management of dogs, especially stray dogs, in case they run wild and bite each other and humans.

Finally, we consider the combined influence of A and k, and a and k on R 0 in Fig. 8 , respectively. From the contour surfaces, we can see that when vaccination, management of dogs and controlling the birth rate of dogs are combined, controlling rabies will be more effective. Moreover, the effect of a is greater than A by comparing the two figures. In conclusion, controlling the population of dogs, reducing the birth rate of dogs, increasing the immunization rate of dogs, improve the management of dogs, enhancing the awareness of people about rabies, and combining these measures are effective measures to control rabies in China. In addition, because the monthly data of human rabies cases exhibits a periodic pattern on an annual base and the human rabies cases in the summer and autumn are higher, it will be useful to take extra measures from May to July every year before the infection peaks, such as extra supervision of children and students out of school.

The transmission of rabies has been a growing concern in China. The data of human rabies cases reported by the Chinese Ministry of Health exhibit seasonal characteristics that the morbidity rates in the summer and autumn are much higher than in the winter and spring. In order to study the transmission dynamics of rabies in China, seasonality of the spreading of the rabies was incorporated into an SEIRS mathematic model with periodic transmission rates. Firstly, we calculated the basic reproduction number R 0 (Diekmann et al. 1990 ) and analyzed the dynamics of the model including the global stability of the disease-free equilibrium and the existence Table 1 do not change Fig. 8 The graph of R 0 in terms of (a) A and k and (b) a and k. Other parameter values in Table 1 do not change of periodic solutions. R 0 was calculated following the definition of Bacaer and Guernaoui (2006) , namely R 0 = ρ(L), where L is the next infection operator, which has been employed in some other studies (Bai and Zhou 2011; Liu 2010; Liu et al. 2010; Fig. 9 The fitting result of the tendency of the infected human rabies cases I 1 (t) using the corresponding autonomous system. A = 232,000, β = 9.56 × 10 −8 , β 1 = 3.68 × 10 −11 and other parameter values are the same as those in the nonautonomous system. The oscillatory curve is the data reported and the skewed curve is the fitting result by the autonomous system Nakata and Kuniya 2010; Wang and Zhao 2008; Zhang and Zhao 2007) . In particular, Wang and Zhao (2008) generalized the techniques and results of van den Driessche and Watmough (2002) to periodic ODE models and provided a recipe to calculate the basic reproduction number. Their results also indicate that R 0 is a threshold value for determining the local stability of the disease-free periodic solution. We then used our model to simulate the monthly data on the number of infected human cases from January 2004 to December 2010 in China reported by the Chinese Ministry of Health and predicted the general tendency of disease in China. Moreover, we carried out some sensitivity analysis of parameters on R 0 . The demographic data were estimated from National Bureau of Statistics of China (NBSC 2009) . The values of most parameters in our model were obtained from the literature or by assumptions. The values of β(t), β 1 (t) and A were estimated through least-square fitting of I 1 (t i ) by discretizing the ordinary differential system as in Chowell et al. (2006) and Stafford et al. (2000) .

Moreover, we discussed and compared the basic reproduction numbers under different conditions. In Sect. 3, we evaluated that R 0 = 1.03 for the nonautonomous model. We can also define and calculate the basic reproduction number for the corresponding autonomous system (see Diekmann et al. 2010; van den Driessche and Watmough 2002) . In this case, the fitting results and parameter values which also are obtained by the least-square method are shown in Fig. 9 . Comparing Fig. 9 with Fig. 2 , we can see that the new model gives a better fit to the monthly human rabies cases.

The basic reproduction number of the corresponding autonomous system iŝ

It can be seen thatR 0 is slightly less than R 0 . From the fitting results and the basic reproduction number, we think that the nonautonomous system is better and biologically more realistic than our previous autonomous model (Zhang et al. 2011) . When studying the transmission dynamics of periodic epidemic models, some researchers use the average basic reproduction numberR 0 , namely the basic reproduction number of the time-averaged autonomous system of the periodic epidemic Fig. 10 The tendency of the human rabies infectious cases I 1 (t) when A = 220,000 and the values of other parameters do not change. Here, R 0 = 0.97 andR 0 = 1.17 model over a time (Greenhalgh and Moneim 2003; Ma and Ma 2006; Moneim 2007; Wesley and Allen 2009; Williams 1997) . However, the average basic reproduction numberR 0 can overestimate or underestimate infection risks. Using the expression of R 0 , we can give the corresponding average basic reproduction number in the nonautonomous system,R

.

Using the parameter values in Table 1 , we haveR 0 = 1.24, which is larger than R 0 . By numerical simulations of our model, we found that the average basic reproduction numberR 0 overestimates infection risks. When A = 220,000 and other parameters values do not change, R 0 = 0.97 andR 0 = 1.17. Note thatR 0 = 1.17 predicts that the disease should be prevalent. However, the curve in Fig. 10 tends to zero. Furthermore, we would like to point out that Zhang et al. (2011) estimated that the basic reproduction number R 0 = 2, but in this paper, we calculated it to be 1.03 with similar model structure. This can be explained as follows. In Zhang et al. (2011) , the data adopted by model fitting were from 1996 to 2010. In this paper we only used the data from 2004 to 2010. From Fig. 11 , we can see that the numbers of infected human cases increased dramatically from 159 cases in 1996 and were fierce from 1996 to 2004. After 2004, the spread of rabies began to slow down. Moreover, the data in 2005 and after 2007 decreased. So the fact that R 0 in this paper is less than in previous paper is reasonable. We also would like to calculate the basic reproduction number for the nonautonomous system from 1996 to 2010. However, we cannot obtain the monthly data before 2004 since only after the SARS outbreaks in 2003 the Chinese Ministry of Health began to publish monthly data on various infectious diseases.

In summary, firstly we have proposed a nonautonomous system with periodic transmission rates to study the transmission dynamics of rabies in China, which is Fig. 11 The data of the number of human rabies cases from 1996 to 2010 reported by the Chinese Ministry of Health more realistic and described the monthly rabies data from China more accurately. Moreover, the nonautonomous differential system is different from and better than the autonomous model used in Zhang et al. (2011) . The dynamics, especially the basic reproduction number, are more complex. Thirdly, by comparing the basic reproduction numbers, we found that it is more realistic to regard R 0 rather thanR 0 or R 0 as a threshold to determine whether the disease can establish or vanish eventually. Finally, we concluded that human rabies in China can be controlled by reducing the birth rate of dogs, increasing the immunization rate of dogs, enhancing the awareness of people about rabies, improving supervision of children and students in the summer, and increasing the medical service and treatment such as PEP after bitten by dogs.

Finally, we would like to comment that most human rabies cases in China are distributed the south provinces such as Guangdong, Guangxi, Guizhou, Hunan, Sichuan. However, the monthly data about human rabies cases we used in our simulations, reported by the Chinese Ministry of Health, were national data and homogeneous mixing was assumed in the model. In future studies, it will be very interesting to obtain regional rabies data (such as provincial data) and more realistic to make inhomogeneous mixing assumption. Reaction-diffusion equation models Wu 2009 and Zhang et al. 2012) or multi-patch models may be more suitable to study the spatial transmission dynamics of rabies in China. We leave these for future consideration. and Zhao (2008) . It is easy to see that system (1) has one disease-free equilibrium P 0 = (Ŝ, 0, 0,R,Ŝ 1 , 0, 0, 0), whereŜ = (m+λ)A m(m+λ+k) ,R = kA m(m+λ+k) andŜ 1 = B m 1 . We rewrite the variables of system (1) as a vector x = (E, E 1 , I, I 1 , S, S 1 , R, R 1 ). Following Wang and Zhao (2008) , we have

So we derive

We introduce some notation. Let Y (t, s), t ≥ s, be the evolution operator of the system

That is, the 4 × 4 matrix Y (t, s) satisfies

for any t ≥ s, Y (s, s) = I , where I is the 4 × 4 identity matrix. Now we introduce the linear ω-periodic system

with parameter z ∈ R. Let W (t, s, z), t ≥ s, be the evolution operator of system (16) on R 4 . Clearly, Φ F −V (t) = W (t, 0, 1), ∀t ≥ 0. Following the method in Wang and Zhao (2008) , we let φ(s) be ω-periodic in s and the initial distribution of infectious individuals. So F (s)φ(s) is the rate of new infections produced by the infected individuals who were introduced at time s. When t ≥ s, Y (t, s)F (s)φ(s) gives the distribution of those infected individuals who were newly infected by φ(s) and remain in the infected compartments at time t. Naturally,

is the distribution of accumulative new infections at time t produced by all those infected individuals φ(s) introduced at time previous to t.

Let C ω be the ordered Banach space of all ω-periodic functions from R to R 4 , which is equipped with the maximum norm · and the positive cone C + ω := {φ ∈ C ω : φ(t) ≥ 0, ∀t ∈ R + }. Then we can define a linear operator L : C ω → C ω by (Lφ)(t) = ∞ 0 Y (t, t − a)F (t − a)φ(t − a) da, ∀t ∈ R + , φ ∈ C ω . L is called the next infection operator and the spectral radius of L is defined as the basic reproduction number

for the periodic epidemic model. To determine the threshold dynamics, we use Theorems 2.1 and 2.2 in Wang and Zhao (2008) . First of all, we verify the seven assumptions in the theorems.

(1)-(5) The first five conditions can be verified by observing F , V + and V − . (6) ρ(Φ M (ω)) < 1, where ρ(Φ M (ω)) is the spectral radius of Φ M (ω) and Φ M (t) is the monodromy matrix of the linear ω-periodic system dq dt = M(t)q with

When M is a constant matrix, M is stable if and only if ρ(Φ M (ω)) < 1. So, we just need to show that M is stable, that is, all eigenvalues are negative. It is obvious that −(m 1 + λ 1 ) and −m 1 are the eigenvalues of M and are negative. Then we consider the eigenvalues of −(m + k) λ k −(m + λ) .

By calculating, −(m + k) − (m + λ) < 0 and [−(m + k)] × [−(m + λ)] − kλ = m 2 + (k + λ)m > 0. So the remaining two eigenvalues are also negative. We can conclude that M is stable, namely ρ(Φ M (ω)) < 1. (7) ρ(Φ −V (ω)) < 1, where Φ −V (t) is the monodromy matrix of the linear ωperiodic system dy dt = −V (t)y with

Because −V is a constant matrix, we need to show that −V is stable. We can see that the eigenvalues of −V are the diagonal elements and negative. So this condition is satisfied.

Using (ii) in Theorem 2.1 in Wang and Zhao (2008) , we derive We see that all eigenvalues of the above matrix are real and the largest one is

So the basic reproduction number is z 0 such that g(z 0 ) = 1.

Rabies knowledge for 20 questions

The epidemic threshold of vector-borne diseases with seasonality

Threshold dynamics of a Bacillary Dysentery model with seasonal fluctuation

Dynamics of measles epidemics: Estimating scaling of transmission rates using a time series SIR model

Rabies-What are the signs and symptoms of rabies?

Rabies answer of knowledge and hot question

Transmission dynamics of the great influenza pandemic of

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

The construction of next-generation matrices for compartmental epidemic models

Seasonal variation in host susceptibility and cycles of certain infectious diseases

Dynamical resonance can account for seasonality of influenza epidemics

A simple model for complex dynamical transitions in epidemics

SIRS epidemic model and simulations using different types of seasonal contact rate

Synchronous cycles of domestic dog rabies in Sub-Saharan Africa and the impact of control effort

Dynamics of rabies epidemics and the impact of control efforts in Guangdong Province

Threshold dynamics for a HFMD epidemic model with periodic transmission rate

A Tuberculosis model with seasonality

Recurrent outbreaks of measles, chickenpox and mumps. i. Seasonal variation in contact rates

Epidemic threshold conditions for seasonally forced SEIR models

Ministry of health of the People's Republic of China, the status of prevention and control of rabies in China (Zhongguo Kuangquanbing Fangzhi Xiankuang

Ministry of health of the People's Republic of China

The effect of using different types of periodic contact rate on the behaviour of infectious diseases: A simulation study

Global dynamics of a class of SEIRS epidemic models in a periodic environment

National Bureau of Statistics of China

Differential equations and dynamical systems

Modeling spatial spread of communicable diseases involving animal hosts

An age-structured model of pre-and pose-vaccination measles transmission

Small amplitude, long periodic out breaks in seasonally driven epidemics

Infinite subharmonic bifurcation in an SIER epidemic model

Multiple stable subharmonics for a periodic epidemic model

The theory of the chemostat

Epidemiological investigations of human rabies in China

Modeling plasma virus concentration during primary HIV infection

Convergence results and a Poincaré-Bendixson trichotomy for asymptotically autonomous differential equations

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Threshold dynamics for compartmental epidemic models in periodic environments

The basic reproduction number in epidemic models with periodic demographics

Human rabies

Infectious disease persistence when transmission varies seasonally

A periodic epidemic model in a patchy environment

Analysis of rabies in China: Tranmission dynamics and control

Spatial spread of rabies in China

Dynamical systems in population biology

Transmission dynamic and economics of rabies control in dogs and humans in an African city

We evaluate the basic reproduction number R 0 for system (1) following the definition of Bacaer and Guernaoui (2006) and the general calculation procedure in Wang