Modeling transmission dynamics of rabies in Nepal.

Even though vaccines against rabies are available, rabies still remains a burden killing a significant number of humans as well as domestic and wild animals in many parts of the world, including Nepal. In this study, we develop a mathematical model to describe transmission dynamics of rabies in Nepal. In particular, an indirect interspecies transmission from jackals to humans through dogs, which is relevant to the context of Nepal, is one of the novel features of our model. Our model utilizes annual dog-bite data collected from Nepal for a decade long period, allowing us to reasonably estimate parameters related to rabies transmission in Nepal. Using our model, we calculated the basic reproduction number ( R 0 = 1.16 ) as well as intraspecies basic reproduction numbers of dogs ( R 0 D = 1.14 ) and jackals ( R 0 J = 0.07 ) for Nepal, and identified that the dog-related parameters are primary contributors to R 0 . Our results show that, along with dogs, jackals may also play an important role, albeit lesser extent, in the persistence of rabies in Nepal. Our model also suggests that control strategies may help reduce the prevalence significantly but the jackal vaccination may not be as effective as dog-related preventive strategies. To get deeper insight into the role of intraspecies and interspecies transmission between dog and jackal populations in the persistence of rabies, we also extended our model analysis into a wider parameter range. Interestingly, for some feasible parameters, even though rabies is theoretically controlled in each dog and jackal populations ( R 0 D < 1 , R 0 J < 1 ) if isolated, the rabies epidemic may still occur ( R 0 > 1 ) due to interspecies transmission. These results may be useful to design effective prevention and control strategies for mitigating rabies burden in Nepal and other parts of the world.

Introduction

Rabies, a viral zoonotic disease existing since at least 2300 BC, is a disease that can affect the nervous system (Yadav, 2012). Despite continuous effort for its control and advancement of medical technologies, rabies still remains a public-health concern in many countries around the world, particularly the countries of Africa and Asia, including Nepal (Devleesschauwer et al., 2016). Rabies has spread across more than 150 countries in all continents except Antarctica. All warm-blooded mammals can contract rabies when the virus-contaminated saliva enters the body through a cut in skin and/or via animal bites. The disease is caused by Lyssavirus, one of the three genera of viruses in the family Rhabdovirus (Center for Disease Control and Prevention Rabies, 2018; Devleesschauwer et al., 2016). If an infected host remains untreated, the virus can transfer to the brain, causing brain inflammation (Abera, Assefa, Belete, & Mekonen, 2015; Hemachudha, Laothamatas, & Rupprecht, 2002) and progression to coma, and eventually leading to death due to respiratory paralysis within ten days after the onset of symptoms (Asamoah, Oduro, Bonyah, & Seldu, 2017; Center for Disease Control and Prevention Rabies, 2018; Ruan, 2017; Yadav, 2012). If the symptoms appear in the infected hosts, the likelihood of death is almost 100% (Center for Disease Control and Prevention Rabies, 2018; Ruan, 2017). The common symptoms of rabies in human include fever, headache, general weakness or discomfort during early stage of the disease (Demirci, 2014), insomnia, anxiety, confusion, partial paralysis, excitation, hallucinations, agitation, hypersalivation (increase in saliva), difficulty swallowing, and hydrophobia (fear of water) during the late stage. Given the persistent rabies burdens in many communities and the lack of diagnostic tools for detecting rabies before the onset of the symptoms, further studies of rabies transmission dynamics are important for the control of disease epidemics.

Rabies usually circulates in two epidemiological cycles: an urban cycle involving dog populations and a sylvatic cycle involving wildlife such as jackals, wolves, foxes, cats, lions, monkeys, mongooses, raccoons, skunks, etc. (Center for Disease Control and Prevention Rabies, 2018; Devleesschauwer et al., 2016; Liu, 2013). In majority of the developing countries, dogs are the main reservoir of rabies virus for human infections (Ruan, 2017). Therefore, most of the previous studies have mainly focused on a single reservoir, i.e. dogs, regarding the urban cycle in isolation (Demirci, 2014; Ega, Luboobi, & Kuznetsov, 2015; Leung & Davis, 2017; Ruan, 2017). However, several wild animals are able to maintain the rabies in the wild independently from dogs (Bellan et al., 2012; Bingham, Foggin, Wandeler, & Hills, 1999; Bingham & Foggin, 1993; Devleesschauwer et al., 2016; Foggin, 1985, pp. 399–405; Singh et al., 2017). Importantly, there are many evidences of the transfer of rabies virus from dogs to wildlife and vice versa, connecting the two cycles of rabies epidemiology (Bellan et al., 2012; Bingham et al., 1999; Devleesschauwer et al., 2016; Gongal & Wright, 2011; Vodopija, Racz, & Pahor, 2016). In many parts of the world, including Nepal, India, Zimbabwe, and South Africa, the jackal in particular is the second most important vector after domestic dogs (Bingham et al., 1999; Vodopija et al., 2016; Yadav, 2012). Thus, studies that involve only a single vector may not be suitable to describe the transmission dynamics of rabies in countries like Nepal, where the wild animal vectors, such as jackals, can be significantly important for disease transmission.

As far as rabies in Nepal is concerned, the transmission of rabies from jackal to human through stray dogs as mediator, i.e. jackals → stray dogs → humans, makes the disease transmission dynamics particularly unique. The notation jackals → stray dogs → humans is used to mean the transmission from jackals to humans through dogs, but not to mean a linear progression. Rabies is endemic in Nepal (World Health Organization, 2014), where the most common biotype of the rabies virus is canine biotype. In Nepal, about 35,000 annual animal bites receiving post-exposure prophylaxis are reported in recent years (Devleesschauwer et al., 2016). The yearly human death due to rabies is documented to be about 100 but this number has been thought to be highly underestimated; the actual annual death may be more than 1,500 every year (Devleesschauwer et al., 2016; World Health Organization, 2014). In Nepal, the sylvatic cycle of rabies is maintained by jackals and foxes (Yadav, 2012). Unfenced forests, open border with India, and frequent contacts of wild animals with domesticated animals primarily increase the transmission risk of rabies from wild animals to stray dogs (Yadav, 2012), which may subsequently transmit the virus to humans. In addition, increase in population of stray dogs is a huge problem in urban areas, exacerbating the threat to rabies outbreak in these regions. Although the number of reported rabies death is relatively low, the rabies control program in Nepal is not very successful (Devleesschauwer et al., 2016). Unawareness about this disease among Nepalese people, frequent ignorance of bites from dogs and animals, passive rabies surveillance from the government, lack of technical facilities in the regional veterinary laboratories, and limited vaccination programs are some of the reasons for poor outcomes of control programs (Devleesschauwer et al., 2016). Also, since dogs are objects of religious worship and other signs of household association (Bögel & Joshi, 1990), controlling the dog population through culling may not be a preferred solution in Nepal. There is, therefore, a need to devise proper strategies for rabies control in Nepal, and mathematical modeling may help to predict outcomes of potential control strategies.

There are many existing mathematical models that describe transmission dynamics of rabies in various contexts (Abo & Twizell, 2006; Anderson, Jackson, May, & Smith, 1981; Asamoah et al., 2017; Clayton, Duke-Sylvester, Gross, Lenhart, & Real, 2010; Demirci, 2014; Ding, Gross, Langston, Lenhart, & Real, 2007; Ega et al., 2015; Hou, Jin, & Ruan, 2012; Leung & Davis, 2017; Liu, 2013; Neilan & Lenhart, 2011; Panjeti & Real, 2011; Rhodes, Atkinson, Anderson, & Macdonald, 1998; Ruan, 2017; Tohma et al., 2016; Zhang, Jin, Sun, Zhou, & Ruan, 2011). For example, Anderson et al. (Anderson et al., 1981) developed a mathematical model to describe rabies dynamics in the fox population and used their model to study rabies control by culling or vaccinating foxes. These simple models have been extended to include two populations, dogs and humans (Asamoah et al., 2017; Zhang et al., 2011), and further extended to include domestic dogs and stray dogs (Hou et al., 2012). Some models were developed to describe spatio-temporal dynamics of rabies spread (Abo & Twizell, 2006; Ding et al., 2007; Neilan & Lenhart, 2011). Some of these models have also been used to formulate optimal control problem to identify optimal strategies to control rabies in raccoon (Clayton et al., 2010). It is worth noting that Rhodes et al. (Rhodes et al., 1998) studied observations of the side-striped jackal (Canis adustus) in Zimbabwe and developed a stochastic model of the transmission dynamics of the virus and host demography. They found that in Zimbabwe, the side-striped jackal population alone cannot explain the observed rabies endemics because of the low jackal population density that is insufficient to maintain the chain of infection. However, they suggest that the disease is regularly introduced to jackals by rabid dogs from populations associated with human settlements, highlighting the fact that humans, dogs, and jackals need to be considered together for accurate modeling of rabies epidemics.

While previous models have provided important insights into many aspects of rabies epidemiology, almost all of the models considered interactions between dog and human populations only, ignoring secondary reservoirs, such as jackals. As discussed earlier, indirect transmission of rabies from jackals to humans through dogs as mediator is a unique feature that needs consideration in the context of rabies epidemics in Nepal, and such feature has not been incorporated in existing models. In this study, we develop a mathematical model that incorporates jackals → stray dogs → humans transmission to describe the rabies epidemics in Nepal. Our model, in the form of a system of nonlinear ordinary differential equations, considers human and dog as well as jackal populations. We obtain some of the model parameters from published literature (Animal Diversity web. Can, 2018; Center for Disease Control and Prevention Rabies, 2018; Kakati, 2012, pp. 1–21; Leung & Davis, 2017; World Health Organization, 2018) and estimate the remaining parameters from a survey data containing dog-bites and rabies outbreaks in Nepal (Yadav, 2012). We use our model to compute the basic reproduction number of rabies transmission in Nepal, and to predict the dynamics of disease prevalence in Nepal. Furthermore, we use our model to evaluate various control strategies, including the anti-rabies vaccine commonly practiced in Nepal. We also performed a detailed analysis to examine the potential contribution of intraspecies and interspecies transmissions to the persistence of rabies.

Methods

The mathematical model

We denote the total population of humans, stray dogs, and jackals as N, N, and N, respectively. Most of the dogs in Nepal roam freely in the street and are not kept inside the house, except a negligible number of dogs owned by some urban people. Therefore, stray dogs in our model include the majority of dogs in Nepal. Since pre-exposure rabies vaccines are not common in humans (Center for Disease Control and Prevention Rabies, 2018), we assume that only exposed humans may get vaccinated. We divide the human population into four mutually exclusive subpopulations: susceptible (S), exposed (E), vaccinated (V), and infected (I) humans. Since dogs are generally administered with both pre- and post-exposure vaccines, we assume that both uninfected and exposed dogs can get anti rabies vaccines. The dog populations are divided into susceptible (S), exposed (E), vaccinated (V), and infected (I) classes. As there is no vaccination for jackals in Nepal, the jackal population is divided into only three mutually exclusive subpopulations: susceptible (S), exposed (E), and infected (I) jackals. Using these notations we have N = S + E + V + I, N = S + E + V + I, and N = S + E + I.

The schematic diagram of the rabies transmission dynamics in Nepal is shown in Fig. 1. We describe the dynamics and interactions among the subpopulations of jackal, dog, and human, using the following system of ordinary differential equations.

Fig. 1

Schematic diagram of the rabies transmission dynamics in Nepal. The bold pointed solid lines represent transfer of individuals to another compartment due to infection or vaccination. The pointed dotted red lines represent the interaction with different compartments that lead to new infections. The shorter pointed solid lines represent recruitment, deaths (natural and disease related), and loss of immunity of vaccinated individuals.

In the model (1), the first three equations represent the rate of change of the susceptible jackal, exposed jackal, and infected jackal populations, respectively. Λ and μ are the rates of recruitment and natural death of jackals, respectively. The susceptible jackals get infected with rabies at intraspecies rate β when they interact with infected jackals and at interspecies rate β when they interact with infected dogs. The exposed jackals progress to the infectious compartment with a rate γ, and the infected jackals die due to disease with a rate δ.

The next four equations represent the rate of change of susceptible, exposed, vaccinated, and infected dog populations, respectively. The susceptible dogs are recruited with a rate Λ and die with a natural death rate μ. They get infected with rabies with rates β (interspecies) and β (intraspecies) after they are bitten by infected jackals and infected dogs, respectively. We denote the rate of dog vaccination (either susceptible or exposed) as u. Note that both susceptible and exposed dogs are vaccinated as dogs are not tested for rabies before the vaccines are administered. The exposed dogs progress to the infected class with rate γ, and the infected dogs also die at the disease induced death rate δ.

The last four equations describe the rate of change of susceptible, exposed, vaccinated, and infected human populations. The human susceptible population is recruited with rate Λ. Susceptible humans are infected through the bites of infected dogs at rate β. The exposed human progress to the infected class with rate γ. The model parameters u, δ, and μ represent rates of post-exposure human vaccination, death due to rabies, and the natural death, respectively.

For the purpose of comparison, we also consider commonly used rabies model, which takes only dog and human populations into account (Zinsstag et al., 2009). Note that with β = β = 0, jackal population related equations are decoupled in our model and the last eight equations become the previously used model with only dog to human interspecies transmission.

Data and parameter estimation

We estimate model parameters and initial populations using various published sources (Animal Diversity web. Can, 2018; Center for Disease Control and Prevention Rabies, 2018; Leung & Davis, 2017; MacLachlan & Dubovi, 2016; Plotkin, 2000) as well as dog-bites and rabies outbreak data (animal and human) from Nepal (Devleesschauwer et al., 2016; Leung & Davis, 2017; World Health Organization, 2018).

Data

In this study we used rabies related data published in Devleesschauwer et al. (Devleesschauwer et al., 2016). The data contains the annual dog bites in Nepal during a time period from 2004 to 2013 (Fig. 2). We acknowledge that, some dog bites, especially bites to children, may go unreported and not all reported dog bites may cause rabies. However, as identified by Devleesschauwer et al. (Devleesschauwer et al., 2016), development of fatal rabies was linked to dog bite incidents in Nepal and the most comprehensive source of epidemiological data on human rabies in Nepal are the annual reports of the Department of Health Services. We obtained these comprehensive available data compiled by Devleesschauwer et al. (Devleesschauwer et al., 2016), since these data are so far the most appropriate proxy of the annual rabies dog bites in Nepal. In addition, the data also includes that about 35,000 people in Nepal get post exposure vaccination annually (Devleesschauwer et al., 2016).

Fig. 2

The data (Devleesschauwer et al., 2016) (blue bars) and the model predictions (red bars) of annual dog bites in Nepal from 2004 to 2013.

Initial populations

To be consistent with the data, we set the beginning of the year 2003 as the initial time of our model dynamics. In Nepal, there are approximately 2 million dogs (Devleesschauwer et al., 2016). As we do not have dog vaccination data available for the year of 2003, we assume that about 20% (4 × 105) of the dog population were vaccinated in 2003. In addition, for our base case, we assume that about 0.5% dogs were exposed to rabies and 0.01% dogs were infected by rabies. The total area of Nepal is about 147, 000 km2, which is divided into mountains (65%), hills (18%), and Tarai (17%) regions (United States Agency for International Development). While very few and scattered jackal-related rabies cases are reported in the mountain and hilly regions, majority of such cases occur in Tarai region (The Himalayan Times, 2016, 2018; Yadav, 2012), i.e., the southern part of the country along the Indian border. The Tarai region contains dense forests, including the national forests, and is habitat of many wild animals. Debnath and Chaudhary (Debnath & Choudhury, 2013) reported that there was on average one jackal per square kilometer in southern Assam, India, which has similar geography and climate conditions as the Tarai region of Nepal. Therefore, taking one jackal per square kilometer in the Tarai region and a lesser density in other regions, we assume on average one jackal per two square kilometers in Nepal, which gives the total 0.5 × 147, 000 = 73, 500 jackals in Nepal. Again, as in the case of dog population, we assume 0.5% jackal are exposed and 0.01% were infected in 2003. The human population of Nepal in 2003 was about 25.31 million (Country Economy, 2019). Based on the data, we assume that about 15,000 human population were exposed, 1000 humans were infected and about 14,000 were vaccinated.

Demographic parameters

According to the WHO, the average human life expectancy in Nepal is 70.2 years (World Health Organization, 2018). Thus, we take the natural death rate as per year. Assuming the steady state level before the epidemic begins, the model implies Λ − μN(0) = 0, which allows us to estimate the recruitment rate for human as Λ = 0.01423 × N(0) individuals per year. The average lifespan of normal dogs is about 12 years (Michell, 1999; Proschowsky, Rugbjerg, & Ersbøll, 2003), but most of the stray dogs in Nepal live for less than 3 years (Leung & Davis, 2017) due to poor environmental conditions. Thus we take the average life expectancy for all dogs in our model as 5 years, implying the natural death rate of a dog as per year. Similarly, the life expectancy of jackals is about 8 years (Animal Diversity web. Can, 2018). Therefore, we take the natural death rate of a jackal is per year. As in the case of a human, we assume the steady state level before the epidemic begins for both stray dogs and jackals. As a result, we get Λ − μN(0) = 0, and Λ − μN(0) = 0. Thus, the recruitment rate for dog and jackal population is Λ = 0.033 × N(0) dogs per year and Λ = 0.125 × N(0) jackals per year, respectively.

Rabies related parameters

After the onset of rabies symptoms, an infected human dies in about 10 days (Center for Disease Control and Prevention Rabies, 2018), implying the rate of rabies-induced death to be per day which gives δ = 36.5 per year. A rabid stray dog can live at most 10 days after the onset of symptoms (Leung & Davis, 2017; Tepsumethanon et al., 2004). Using this information for rabid dogs, the rabies induced death rate for dogs is δ = 36.5 per year. We also assume the same rate for rabies related mortality rate for jackals. The incubation period is defined as the time from virus transmission until the symptoms appear. The incubation period of rabies in dog can vary from few days to more than a year (VCA Animal Hospital, 2020). Thus, we take the average incubation period in dogs to be 6 months, implying that the rate at which the rabies exposed dogs begin to show the symptoms is γ = 2 per year. Similarly, since the rabies incubation period in humans might be from few weeks to more than a year (World Health Organization, 2020a), we take the average of 6 months with γ = 2 per year. The incubation period of rabies in jackals has been found to be about 55 days (MacLachlan & Dubovi, 2016), giving γ = 6.64 per year.

Model fitting to the data

With the fixed parameters and initial conditions discussed above, we estimate the remaining parameters, transmission rates and basic vaccination rates in human and dog population, by fitting the model to the data. We solve the model system in MATLAB using subroutine ode45, and use the solution to minimize the following relative error:where n is the number of data points, is the number of dog-bites at time t documented in the experimental data, and y(t) is the model prediction of dog-bites at time t. For minimization, we use a combination of fmincon and multistart algorithm in Global Optimization Toolbox of MATLAB (Pantha, Cross, Lenhart, & Day, 2018; The MathWorks Inc Global, 2015a). The multistart algorithm uses several uniformly distributed starting points for performing fmincon optimization solver, allowing a thorough search for a global minimum.

Results

Model selection and parameter estimation through data-fitting

While our model includes both dog and jackal populations motivated from the evidences of rabies in both species, we also considered the commonly practiced single animal species (dog only) model to fit the available data from Nepal. For comparing these models, we implemented the Akaike’s Information Criterion (AIC) method, which has been widely used in previous studies, including the ones with differential equation models (Barker & Vaidya, 2020; Burnham, Anderson, & Huyvaert, 2011; Akaike, 1974; Portet, 2020). Specifically, we fitted both models to the data by using the techniques described in Section 2.2.5, and computed the corresponding values of AIC, which is given by the formula (Akaike, 1974; Barker & Vaidya, 2020; Burnham et al., 2011; Portet, 2020),

In this formula for AIC, n and n represent the number of data points and the number of parameters estimated through data-fitting, respectively, and SSR is the sum of the squared residuals, i.e., , where y(t) and are the model predicted value and the experimental value, respectively, of the number of dog-bites at time t. Our data fitting revealed that the dog-human model (single animal species model) produces significantly high AIC values (AIC = 325) compared to the jackal-dog-human model (two animal species model) (AIC = 312). Thus, the model containing both jackal and dog populations that we considered in this study can describe this data set from Nepal better and is more suitable to represent the rabies transmission dynamics in Nepal.

The fixed parameters and initial populations as well as estimated parameters through our model fitting (the jackal-dog-human model) to the data are given in Table 1. With these parameters, the model predictions agree well with the data collected from Nepal for a decade long period from 2004 to 2013 (Fig. 2). From the data fitting, we estimated the transmission rates as β = 3.79 × 10−5, β = 1.90 × 10−5, β = 1.52 × × 10−5, β = 2.74 × 10−5, and β = 1.71 × 10−6 per individual per year. Similarly, for dogs and humans, we estimated the vaccination rates to be u = 0.034 per year and u = 2.05 per year, respectively.

Table 1

The model variables and parameters.

1a. Model Variables and Initial Values
Variable	Description	Initial values	Source
SJ	Susceptible Jackals	73125	Assumed
EJ	Exposed Jackals	368	Assumed
IJ	Infected Jackals	73	Assumed
SD	Susceptible Dogs	15.898 × 105	Devleesschauwer et al. (2016)
ED	Exposed Dogs	104	Assumed
VD	Vaccinated Dogs	4 × 105	Assumed
ID	Infected Dogs	200	Assumed
SH	Susceptible	25.265 × 106	Country Economy (2019)
EH	Exposed	15534	Assumed
IH	Infected	1000	Assumed
VH	Vaccinated	14000	Assumed
1b. Demographic Parameters
Params.	Description	Value	Source
μJ	Jackal mortality rate	0.125	Animal Diversity web. Can (2018)
μD	Dog mortality rate	0.2	Leung and Davis (2017)
μH	Human mortality rate	0.0142	World Health Organization (2018)
ΛJ	Jackal recruitment rate	μJ × NJ(0)	Animal Corner (2018)
ΛD	Dog recruitment rate	μD × ND(0)	Kakati (2012)
ΛH	Human recruitment rate	μH × NH(0)	World Health Organization (2018)
1c. Rabies Related Parameters
Params.	Description	Value	Source
δJ	Jackal rabies related mortality rate	36.5	Assumed
δD	Dog rabies related mortality rate	36.5	Leung and Davis (2017)
δH	Human rabies related mortality rate	36.5	Center for Disease Control and Prevention Rabies, (2018)
γJ	Jackal rate of moving from exposed to infected	6.64	MacLachlan & Dubovi (2016)
γD	Dog rate of moving from exposed to infected	2	VCA Animal Hospital (2020)
γH	Human rate of moving from exposed to infected	2	World Health Organization (2020a)
1d. Transmission and Control Parameters
Params.	Description	Value	Source
βJJ	Transmission rate from Jackal to Jackal	3.79 × 10−5	Estimated
βDJ	Transmission rate from Dog to Jackal	1.90 × 10−5	Estimated
βJD	Transmission rate from Jackal to Dog	1.52 × × 10−5	Estimated
βDD	Transmission rate from Dog to Dog	2.74 × 10−5	Estimated
βDH	Transmission rate from Dog to Human	1.71 × 10−6	Estimated
uD	Dog vaccination rate(susceptible and exposed)	0.03	Estimated
uH	Human vaccination rate (PEP)	2.05	Estimated

Basic reproduction number

The basic reproduction number, , is defined as the expected number of secondary cases produced by a single (typical) infection in an entirely susceptible population (van den Driessche & Watmough, 2002). For infectious disease dynamics, is an important threshold which determines whether the outbreak occurs () or infection dies out () (van den Driessche & Watmough, 2002). We use the next generation matrix method (van den Driessche & Watmough, 2002) to derive expression for the basic reproduction number, , from our model. According to this method, we first compute the disease free equilibrium (DFE) of our system as

We then consider the infected equations of the system, i.e. a subsystem containing all of the equations except the equations for S, S, V, S, and V. This subsystem is then linearized about the DFE. From the resulting equations, we form two matrices, a matrix containing infection terms and a matrix containing transfer terms, as given belowandwhere

The basic reproduction number is then given by the spectral radius of the matrix . Therefore,where

Intraspecies basic reproduction numbers. We also consider the cases in which species (jackal, dog, humans) are isolated and interspecies transmissions are ignored, i.e. β = β = β = 0. In this case, the model system results in three decoupled subsystems for the dog, jackal, and human, and we refer the basic reproduction numbers corresponding to these subsystems as the intraspecies basic reproduction numbers. For the decoupled systems corresponding to the dog and jackal, we again use the next generation matrix method (van den Driessche & Watmough, 2002) to obtain the basic reproduction number for each species. We obtain and (defined above in equation (2)) as the intraspecies basic reproduction numbers for dogs and jackals, respectively.

Using parameters from Table 1 in the expression of , we estimate the basic reproduction number for rabies in Nepal as , while intraspecies basic reproduction numbers are , and . indicates the rabies spread in Nepal, consistent with the reported on-going epidemics, although some uncertainty remains in the data collected. The value of for our base case is maintained at a low level due to the current level of the dog vaccination. However, in the absence of such vaccination (i.e., u = 0) with all other conditions remaining the same, the basic reproduction number for Nepal could have reached as high as . This indicates that the currently implemented dog vaccination program is helpful in curbing the rabies incidence but is not sufficient enough to avoid the epidemic and/or eradicate the disease in Nepal. Note that the estimated value of is similar to the basic reproduction number provided in previous studies of rabies spread in other places (Ruan, 2017; Zhang et al., 2011).

To identify the parameters that have more impact on the basic reproduction number , we perform sensitivity analysis by calculating the sensitivity index S given bywhere α is the parameter for which the sensitivity is to be studied (Mutua, Barker, & Vaidya, 2017). The larger the magnitude of the value of S, the bigger the impact of the parameter α has on . The sign of S determines whether the parameter has a positive or negative impact on . The sensitivity index of can also be helpful to determine effective strategies for the control of rabies epidemics. The value of the sensitivity index for each of the model parameters is presented in Fig. 3.

Fig. 3

Sensitivity index, S, related to for each parameter. The parameter with a larger value of S has a bigger impact on . The sign of S determines whether an increase in the parameter value increase (positive) or decreases (negative) the values of .

We observe that the basic reproduction number is highly sensitive to dog recruitment rate Λ, dog transmission rate, β, rabies induced death rate, δ, dog natural mortality rate, μ, and dog vaccination rate, u (Fig. 3). Apart from these parameters with the largest sensitivity index, the second set of parameters to which is sensitive are: the rate at which exposed dogs shows clinical symptoms, γ and the dog vaccination rates, u. The higher sensitivity of to the parameters related to dog populations, as predicted by our model, indicates that the dog is the major vector in the human rabies outbreak. This also suggests that control measures targeted to alter the dog related parameters might be most effective in order to bring the basic reproduction number below one resulting in the rabies eradication in Nepal. However, some jackal related parameters may also have role in reducing , but in lesser magnitude.

Disease outcome: base case scenario

In this section, we present simulation results of our model with the baseline parameter values given in Table 1. The prevalence of rabies among dogs and jackals over time predicted by our model are shown in Fig. 4. Here, the prevalence is calculated as the percentage of the total infected individuals among corresponding species at a given time. For example, the rabies prevalence among the dog population at a time t is calculated as . We observe that the rabies prevalence among the dog population increases steadily and reaches the peak (0.09%) in about 10 years and then decreases until 17 years time point (Fig. 4). After this period, the prevalence eventually levels out at the steady state of about 0.08%. In the jackal population, the disease prevalence peaks at the level of 0.09% in around 7 years and then decreases until it reaches the steady state at 0.07% in about 17 year (Fig. 4). These results support the finding that rabies is currently endemic in Nepal (World Health Organization, 2014). Our model predicts that in the steady state, the total number of new human infections reaches about 55,000.

Fig. 4

The prevalence of rabies in Nepal. The rabies prevalence for 35 years among dog (a) and jackal (b) population predicted by our model with baseline parameter values from Table 1.

At the steady state endemic levels of rabies, we computed the proportion of new infections generated due to intraspecies and interspecies transmissions in animals (Fig. 5). We found that about 98% of the new rabies infections in dogs are due to infected dogs (intraspecies) while only about 2% of the new infections are due to infected jackals (interspecies). Higher intraspecies than interspecies transmission among dogs is expected due to the high density of stray dog population and their unlimited open interaction with other stray dogs. In the jackal population, about 94% of the new infections are due to interaction with the infected dogs (interspecies) while remaining 6% of the new infections are due to intraspecies interaction among jackals. Unlike in the dog population, higher interspecies infection in jackal population could be due to the fact that jackals live in their own territory or home range, causing low interactions with jackals in other territories and/or due to a higher number of infected dogs than jackals.

Fig. 5

New infection generated from intraspecies and interspecies transmission. The percentage contribution from interspecies (blue) and intraspecies (yellow) rabies transmission to annual new infections of dog and jackal population.

Analysis of parameter uncertainty

While our base parameters are estimated from thorough literature review and data-fitting, there remains a quite uncertainty on some of the parameters, especially because of lack of detailed study and data on jackal infection. To analyze how the parameter uncertainty affects model predictions, we vary these parameters over a wide range, and compute the rabies prevalence as well as the basic reproduction numbers. In particular, we generated 500 sets of parameters by choosing each parameter randomly from uniformly distributed values within ±10% of the base values. For each set of generated parameters, we simulated our model until it approached a steady state, and recorded the steady state dog and jackal prevalence (%). Also, we computed the values of the dog and jackal intraspecies reproduction numbers () as well as the overall reproduction number (). The resulting boxplots for the prevalence and the reproduction numbers are shown in Fig. 6. We obtain that the dog has the median prevalence of 0.08% with the first quartile 0.05% and the third quartile 0.10%. Similarly, jackal has the median prevalence of 0.08% with the first quartile 0.05% and the third quartile 0.09%. For the intraspecies reproduction number of dog, the median value remains 1.15 with first and third quartiles 1.08 and 1.2, respectively. Similarly, for the intraspecies reproduction number of jackal, the median value is 0.07 with first and third quartiles 0.07 and 0.08, respectively. In our computations, the overall basic reproduction number has median 1.17 with first and third quartiles 1.11 and 1.22, respectively. This implies that the model prediction does not change widely with variation of parameter values within ±10% range, establishing the robustness of our model.

Fig. 6

Analysis of parameter uncertainty. The boxplot of rabies prevalence among dogs (first), rabies prevalence among jackals (second), intraspecies reproduction number for dog (third), intraspecies reproduction number for jackal (fourth), and overall reproduction number (fifth), obtained from the model simulations with 500 sets of randomly generated parameters.

We note that we also performed sensitivity on the values of initial subpopulations. While there is a minor effect of varying initial values over a reasonable range on the beginning trend of the epidemic, the overall epidemic dynamics and the main conclusion of the paper remain unaffected by the different values of the initial subpopulations.

Prevention and control of rabies in Nepal

In this section, we use our model to evaluate potential prevention and control strategies for mitigating the rabies burden in Nepal. We primarily consider dog vaccination, dog sterilization, dog culling, and jackal vaccination through bait. By implementing these strategies into our model, we compute the expected disease outcome in terms of long-term rabies prevalence among dogs and jackals as well as annual human rabies cases in the steady state. We evaluated each control strategy by varying the related parameters in our model while keeping all other parameters as in the base case (Table 1), including the currently estimated value of dog vaccination rate (u) and human post vaccination rate (u).

Dog vaccination

Dog vaccination is one of the commonly practiced methods to control rabies in Nepal. The parameter u of our model can be used to incorporate the vaccination coverage of the dog population. With the vaccination rate of u, the target dog population, X, can be approximated using the solution ofi.e, . Then for a program that aims to vaccinate η% of the dog population in a period of t years, the vaccination rate u needs to be in such a way that , which gives .

Using this relation, we now compute rabies prevalence in dog and jackal populations as well as the total annual human rabies cases (Fig. 7) for the dog vaccination coverage, η, varying from 0 to 15% targeted to accomplish in a year. In the absence of dog vaccination (i.e, η = 0%), the prevalence among dog and jackal establishes at about 0.18% and 0.14%, respectively, causing about 102,400 annual human rabies cases. An increase in the dog vaccination coverage decreases the prevalence of rabies among both the dog and jackal population as well as the annual human rabies cases. For example, with the yearly dog vaccination coverage of 5%, the rabies prevalence among both dogs and jackals decreases to 0.04%, respectively. In this case, the total number of human rabies cases is about 28,190.

Fig. 7

Effect of dog vaccination program. The rabies prevalence among dog (a) and jackal (b) and the annual human rabies cases (c) under various monthly dog vaccination coverage.

Dog culling

In our model, the per capita rate of dog culling, ν, can be introduced through μ → μ + ν. As discussed in the case of dog vaccination above, the culling rate can be estimated as , where ξ is the percentage of dog culling/removal targeted to achieve in a period of t years.

In Fig. 8, we present the rabies prevalence in dog and jackal population as well as the total number of annual human rabies cases under various annual dog culling rates with all other parameters, including the dog vaccination rate (u) and human vaccination rates (u), taken from the base case scenario 1. The model predicts that without dog culling program, about 0.077% of the dog and 0.073% of the jackal remain prevalent with about 55,000 of the annual human rabies cases. An increase in the coverage of dog culling significantly decreases the prevalence of rabies animal populations as well as human cases and human deaths. For example, a culling program that covers 3% of dogs annually can bring the rabies prevalence among the dog population to almost 0.011% and among the jackal population to almost 0.010% in a long term. In this level of culling, the annual human rabies cases comes down to 7,983 which is about 85% lower than the annual human cases in the absence of culling.

Fig. 8

Effect of dog culling program. The rabies prevalence among dog (a) and jackal (b) and the annual human rabies cases (c) under various per year removal/culling coverage of infected dogs.

Dog sterilization

We now study the effects of dog sterilization on the rabies transmission dynamics. This strategy of controlling stray dogs is usually applied as an alternative to the dog culling strategy. Note that the sterilization does not remove the current dog population or block the infection, it rather stopsthe dog’s reproductivity. As a result, the effect of this program on reducing the dog population, including those infected with rabies, is normally seen late in the future.

To incorporate sterilization effect into our model, we take the recruitment rate Λ = μ × N, where N is the total number of dogs and μ represents the per-capita birth rate. We assume that the sterilization rate k reduces the per-capita birth rate giving the recruitment rate as Λ = (μ − k) × N. With this growth rate, the dog population at time t reaches . This implies that for a sterilization program that aims to reduce the dog recruitment by θ% in a period of t years, the sterilization rate should be . We define θ as the t-year sterilization strength.

In Fig. 9, we present the model prediction for the prevalence of rabies in the dog and jackal population as well as the annual human rabies cases for varying 5-year dog sterilization strength. As mentioned above, the baseline values for the dog vaccination rate (u) and the human vaccination rate (u) along with other parameters are used from Table 1. Without the sterilization effect, we observe that about 0.077% of the dog and 0.073% of the jackals remain prevalent with the annual human rabies cases of 55,000. An increase in the sterilization strength significantly decreases the prevalence of rabies animal populations as well as human cases. For example, 5-year sterilization strength of 40% decreases the rabies prevalence in dogs and jackals each to 0.021%, with the annual human rabies cases decreased to 15,420.

Fig. 9

Effect of dog sterilization program. The rabies prevalence among dog (a) and jackal (b) population and annual human rabies cases (c) under various 5-year dog sterilization strengths.

Jackal vaccination through bait

In many other countries, one of the ways that is frequently practiced to control rabies in wild animals is vaccination through baits, in which the food/bait containing the vaccine is spread throughout different locations so that the target species eat the food and become immune (North Carolina Department of Health and Human Services, 2018). To evaluate whether this potential program can be effective in the context of Nepal, we extended our model by introducing an additional compartment V representing the jackals which become immune due to the vaccination through baits. Denoting the vaccination rate u and implementing the similar method as in the dog vaccination, we get if the program aims to achieve immunity coverage in ζ% of jackal population in a period of t years.

Varying the yearly jackal vaccination coverage ζ from 0 to 70%, we found only a negligible impact of this program on human cases and prevalence in dogs, while some effect was observed in prevalence in jackals (Fig. 10). Specifically, the rabies prevalence in the dog population has no significant decrease (0.077% for 0% monthly coverage and 0.068% for 70% monthly coverage) whereas the prevalence in the jackal population is reduced from 0.0730% to 0.052% when the monthly jackal immunity coverage is increased from 0% to 70%. In this case, the reduction in the number of annual human rabies cases is small (only 5620 reduction).

Fig. 10

Effect of jackal vaccination program. The rabies prevalence among dog (a) and jackal (b) as well as the annual human rabies cases (c) for various yearly coverage of jackal vaccination through bait.

Interspecies and intraspecies transmissions: comparison with a wider parameter range

It is interesting to note that with the base case parameter for Nepal, the jackal intraspecies reproduction number is estimated to be , but the rabies persists at the steady state in jackal population (Fig. 4). Clearly, the persistence of rabies in jackal population is due to the interspecies transmission from the dog (Fig. 5). To study such possibility of persistence of rabies due to interspecies transmission, we extend our model to a wider parameter range and compare the potential role of interspecies and intraspecies transmission on rabies outbreak. For the purpose of demonstration, we consider ranges of parameters in such a way that and span over 0 to 1.5, i.e., and .

As derived above (Eq. (2)), the intraspecies basic reproduction numbers ( for dogs and for jackals) are related to overall basic reproduction number, , through parameters C and E, which contain interspecies transmission rates β and β. Using this relationship between intraspecies reproduction numbers () and the overall basic reproduction number () obtained from our model, we generate a 3-D diagram (Fig. 11a) showing as a function of and . R0 < 1 and R0 > 1 are separated by the horizontal plane at the level R0 = 1. Projecting this 3-D graph into the 2-D plane of () (Fig. 11a), we identified regions corresponding to (green region, no outbreak) and to (yellow region, rabies outbreak). We also identified the blue region where R0 > 1 but and .

Fig. 11

(a) The overall basic reproduction number, R0, as a function of individual basic reproduction numbers, and . The green region represents the area where all of the R0, and are less than 1. The yellow region represents area where one of the and , as well as R0 are greater than 1. The region in blue represents the area where each of and is less than 1 but the value of is greater than 1. (b) The bifurcation diagram of with the increasing interspecies transmission rates while all other parameters are taken in their baseline values. (c, d) The rabies prevalence among dog and jackal population with no interspecies transmission i.e, β = β = 0. (e, f) The rabies prevalence among dog and jackal population with interspecies transmission rates in their base line values i.e, β = β≠0. All the other parameters are taken as baseline parameter values from Table 1.

In the yellow region (Fig. 11a), at least one of the intraspecies basic reproduction numbers (, ) and the overall reproduction number () are greater than 1. This implies that with the parameters in this region the disease outbreak occurs (in both cases with species isolated or coupled). Similarly, we obtain opposite behavior for the green region, in which , , and are less than 1, and the outbreak is avoided in both cases with species isolated or coupled.

The blue region is particularly an interesting one, in which each of and is less than one but is greater than one implying that the disease outbreak occurs despite the absence of outbreak when each species is considered in isolation. In this case, the rabies outbreak does not occur among only dog or only jackal community if isolated from each other, but the combination of dogs and jackals together may lead to the disease outbreak. The size of the blue region depends on the value of the interspecies transmission rates, β, β; the higher the interspecies transmissions, the larger the size of the blue region (Fig. 11b). These results associated with the interspecies transmission underscore the importance of including multiple reservoirs with interspecies transmission in the model.

In order to better illustrate the important role of the interspecies transmission in the rabies dynamics, we now choose a set of parameters from the blue region in Fig. 11a, which correspond to the intraspecies basic reproduction numbers . In this case, if interspecies transmission was neglected, i.e., β = β = 0, the rabies got eliminated in each isolated dog and jackal populations (Fig. 11c and d), as expected. However, if the interspecies transmissions were introduced at their base values, i.e., β = 1.52 × 10−5, β = 1.90 × 10−5, the rabies persisted in both dog and jackal populations (Fig. 11e and f) despite and . This illustrates that although the intraspecies reproduction numbers are less than one asserting the condition for eradication of rabies in isolated populations, there is a possibility of overall reproduction number to be greater than one causing the disease persitence due to the interspecies transmission. Importantly, this implies that control strategies, focused on a single species, which aim to reduce the individual values of and to less than one, may not be enough to make the overall reproduction number less than one. Therefore, an ideal control strategy should target both interspecies and intraspecies transmissions for successful eradication of rabies.

Discussion

Rabies - one of the neglected infectious diseases - still poses a significant burden to public health in many countries, including Nepal. In Nepal, rabies is maintained by two epidemiological cycles: an urban-cycle involving stray dogs and a sylvetic cycle involving jackals (Devleesschauwer et al., 2016). While more than 96% of human rabies cases in Nepal are from the direct bites of dogs (Yadav, 2012), jackals mainly serve as the secondary source of infection to the urban cycle (Devleesschauwer et al., 2016). Therefore, contributions of the jackal population to the persistence of rabies and to the indirect transmission to humans through the intermediate dog population are the primary obstacle to the eradication of rabies in Nepal (The Himalayan Times, 2016, 2018). To describe this less apparent transmission mechanism, namely indirect transmission from jackals to humans via dogs, we developed a novel mathematical model in the form of system of ordinary differential equations, and used our model to describe transmission dynamics of rabies in Nepal as well as to evaluate potential control strategies.

The outcomes predicted by our model is consistent with the rabies outbreak data collected from Nepal for a decade long period from 2004 to 2013 (Devleesschauwer et al., 2016). We also estimated the basic reproduction number of rabies in Nepal to be , which is similar to the values found in rabies outbreak in other countries (Asamoah et al., 2017; Ega et al., 2015; Leung & Davis, 2017; Liu, 2013; Ruan, 2017; Zhang et al., 2011), where the models with only dog-human interaction were used. indicates that the rabies outbreak occurs in Nepal, consistent with the observed data, although some uncertainty remains in the data collected. While we acknowledge some uncertainty also on the estimated parameters, our model predicts that the current situation of rabies burden in Nepal could be a lot larger than that documented in the available reports.

Our model allowed us to draw several interesting conclusions regarding rabies epidemics. Results from our model simulations show that within 7 years of the outbreak, rabies prevalence in both the dog and jackal population reaches its peak and then declines until it reaches its steady state in about 25 years. The disease dynamics stabilizes at the positive steady state levels (0.08% and 0.07% prevalence in dogs and jackals, respectively). These results suggest that the rabies in Nepal is in endemic state, and as a result, human cases and potential rabies-related human deaths, though small in number, are persistently ongoing in Nepal. Because of a higher prevalence as well as a higher population size of dogs, at the steady state level, the contributions of interspecies transmission from dogs to jackals is higher than transmission from jackals to dogs.

We also used our model to evaluate various prevention and control strategies, and identified that the prevention and control programs focused on the strategies related to the dogs are more effective. In fact, currently practiced control strategies in Nepal also focus on the dog population, such as vaccination, sterilization, and culling of infected dogs. Though currently not practiced in Nepal, we also evaluated the vaccination in jackals through baits, which are used in many developed countries (North Carolina Department of Health and Human Services, 2018; United States Department of Agriculture Animal and Plant Health Inspection Service, 2018). According to our model simulations, this strategy appears to be less effective in the context of Nepal because most of the rabies transmission in humans and dogs are due to dog bites. Based on all of these results from our model, we suggest that while each of the dog-targeted strategies can be effective to control a rabies outbreak if implemented individually in a sufficiently high level, combined strategies may be needed to ultimately eradicate rabies in Nepal. While our model predicted that the strategies targeted at dog population is effective, it is important to note that the implementation of these strategies may have some obstacles and limitations. For example, dog culling program can be harder in some places because of the religious reasons; dogs are worshiped and celebrated for their friendship and loyalty in many Hindu communities. There may also be likely increase in puppy survival due to a relative increase in resources. It is worth noting that the World Health Organization (WHO) suggests that killing many dogs regardless of their infection status may not be helpful to eliminate dogs in a long run, but combining the dog vaccination with removal of unvaccinated stray dogs can be effective measures to control the rabies spread successfully (World Health Organization, 2020b).

The basic reproduction numbers derived based on our model implies that to the overall basic reproduction number, , there are significant contributions from interspecies transmissions, C (dog to jackal) and E (jackal to dog), in addition to the contributions from intraspecies transmissions, (intraspecific basic reproduction number in the isolated dog population) and (intraspecific basic reproduction number in the isolated jackal population). This relationship among reproduction numbers, along with the detailed bifurcation analysis performed over a wide parameter range (Fig. 11), shows that there can be a region in the parameter space, for which is greater than one due to the interspecies transmission even though individual values of and are less than one. In this case, prevention strategies, which have the potential to avoid rabies spread in dog and jackal populations if the transmission is restricted to a single species, may still fail to control rabies due to interspecies transmission. Therefore, in the countries like Nepal, where the urban cycle (stray dogs) and the sylvetic cycle (jackals) are in close proximity, the interspecies transmission can be an obstacle to the eradication of rabies and must be taken into account while devising prevention and control strategies.

We acknowledge several limitations of our study. Our parameter estimation is based on limited data sets, which lack detailed rabies cases, particularly in the jackal population. We did not consider the immunity loss of rabies vaccine, which may have some impact on the disease dynamics. Moreover, the previous reports (Devleesschauwer et al., 2016; World Health Organization, 2014) indicate that the data we used may be highly underestimated, especially because of the lack of mechanism to identify rabies in rural areas. Further studies with more detailed data may improve estimation of these parameters and our model predictions. The other limitations in our model include the lack of potential movement of animals from India, where the world’s largest rabies related deaths (more than 20,000 per year) and much higher jackal bite incidents occur. For simplicity, our model assumes that the same strain of rabies virus circulates in both jackals and dogs in Nepal. The model can similarly be extended to multiple strains, but requires more data with multiple strains for its validation.

In summary, the model developed in this study to describe rabies transmission in the context of Nepal underscores the potential importance of interspecies transmission from jackal to dog to human, which may cause persistent rabies epidemics in Nepal. We believe this model could be valuable to devise strategies to control rabies burden in Nepal if used in conjunction with detailed rabies epidemic data from the field.

Declaration of competing interest

Authors declare no conflict of interest.