Dog Rabies in Dhaka, Bangladesh, and Implications for Control

Controlling rabies among free-roaming street dogs has been a huge challenge in many parts of the world. Vaccination is a commonly used strategy to control rabies, however, sufficient vaccination coverage is very challenging when it comes to street dogs. Also, dog rabies data is scarce, making it difficult to develop proper strategies. In this study, we use a logistic growth incorporated epidemic model to understand the prevalence of rabies in the dog population of Dhaka, Bangladesh. The study shows that, the basic reproduction number for dog rabies in Dhaka lies between 1.1 to 1.249 and the environmental carrying capacity lies approximately between 58,110 to 194,739. Considering the vaccination and neuter programs administered in the last decade, we attempt to explain rabies transmission among dogs in this population. We found that the high basic reproduction number is associated with high environmental carrying capacity and vice versa. Further, we compare different type of control strategies, viz., constant vaccination, pulse vaccination, and optimal vaccination strategies. In the case of high environmental carrying capacity, vaccination, and neuter strategy is not sufficient for controlling rabies in street dogs, whereas carrying capacity control through waste management coupled with vaccination and neuter is more effective.


Introduction
Rabies is a zoonotic viral infection that is transmitted through the bites and scratches of foxes, dogs, bats, raccoons, cats, etc. It infects the central nervous system and ultimately leads to a fatal infection of the brain and spinal cord in almost all the cases if not treated before clinical manifestations occur. Rabies takes a heavy toll on human lives every year around the world, mainly in Asia and Africa, and more than 99% of such fatalities are attributed to dog-mediated rabies [1]. According to a data driven empirical study conducted in 2008, rabies kills about 2100 people every year in Bangladesh [2], which was 1.35 times higher than the World Health Organization (WHO) prediction.
Rabies is associated with certain dog-keeping practices and responsible ownership which vary with the cultural and economic conditions of a society [3][4][5][6]. Wandeler classified dogs into three categories: transmission of rabies have been proposed in recent years [12,14,16,[31][32][33][34][35][36][37]. These studies have focused on the transmission of rabies within dogs and to humans and suggested different control strategies to limit this transmission. Zhang et al. [31] studied the transmission of rabies within domestic dogs and from domestic dogs to humans in China, and the authors suggested that reducing the dog birth rate and increasing dog immunization coverage may be effective ways to control rabies. Another mathematical model focusing on domestic and free-roaming dogs suggested that rabies incidence can be reduced by increasing the vaccination rate of domestic dogs and reducing the free-roaming dog population [14]. In another recent sophisticated study conducted by Shigui Ruan [12], the author introduced a rabies transmission model with increasing complexity.
The data driven study of epidemics and population dynamics is based on time series data of infection prevalence. However, no such data are available for rabies in FRSDs. So, using mathematical modeling tools, we inferred key parameters of the dog population and its rabies dynamics.
In this study, we investigated the transmission dynamics of rabies within dogs in the presence of different vaccination and sterilization strategies to investigate the effectiveness of vaccination in stopping the endemic spread. Besides the CNVR program, we also considered food waste management as another control measure to reduce the carrying capacity. We then analyzed the model to investigate whether these combined control measures could be effective in stopping the spread or significantly reducing the incidence of rabies that is mainly caused by FRSDs. To the best of the authors' knowledge, these combined control measures have not been considered in any previously proposed mathematical models addressing rabies transmission in dogs, especially in free-roaming street dogs.

Basic Model
Dogs in low and middle income countries are usually free roaming street dogs living on food waste, which is a determinant of the carrying capacity of the dog population [7,13,15]. In our model, we divide dogs into three compartments: susceptible (S); exposed (E); infectious (I). We assume the logistic growth of dog population relates the growth rate with carrying capacity. We assume that only susceptible dogs give birth to susceptible newborns at the rate r. The carrying capacity is D. The growth term is modeled in such a way that the number of newborns depends on susceptible dogs only, while all dogs contribute to the limitations in available resources. The per-capita transition rate, β, is a product of the contact rate and transmission probability at each contact. In addition, homogeneous mixing of dogs is considered, so the probability that a susceptible dog will meet an infected dog is considered to be I(t) N(t) . Therefore, the transition rate of dogs from the susceptible to exposed compartment is βS(t) N(t) . Dogs remain exposed for a period of 1 σ years on average. Following [31], we assume, γ fraction of the exposed dogs (E) become infectious at the rate γσ, while the rest do not become infectious and instead return to the susceptible compartment at the rate (1 − γ)σ. The parameters m and µ represent the normal and disease related death rates, respectively. With these assumptions, the rabies transmission model is given by where N(t) = S(t) + E(t) + I(t). The model (1) has a disease-free equilibrium(DFE), E 0 = {S 0 , E 0 , I 0 } = (D, 0, 0). The basic reproduction number, R 0 , is the average number of new rabid dogs due to a single rabid dog introduced in DFE. We use the next generation matrix method and follow the notations developed by Van den Driessche and Watmough [32] in calculating the basic reproduction number.
F E 0 V −1 E 0 termed as next generation matrix [32] gives the expected number of secondary infected dogs in each compartment produced by a single rabid dog. The spectral radius of F E 0 V −1 E 0 is the basic reproduction number R 0 = βσγ (µ+m)(m+σ) , which is reasonable from the perspective of an epidemic due to the following.
×probability of contact with susceptible dog × probability of developing infection.
The eigenvalues of the Jacobian of the system (1) at DFE E 0 are: −r and here, the first value is clearly negative, while the other two are negative if R 0 < 1 and at least one is positive if R 0 > 1. Therefore, E 0 is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1.
In addition, the model has another equilibrium , which is biologically relevant only when R 0 > 1. Here, The expression of N * shows that when R 0 = 1, the population reaches the carrying capacity; however, if R 0 > 1, then the short lifespan of infected dogs hinders the total population to reach the carrying capacity. If R 0 > 1 + r(µ+m+σγ) (µ+m)(m+σγ) , then the dog population will become extinct. A similar indication can be observed from the stability of E * .
A stability investigation of E * requires tedious mathematical analysis; therefore, we used the numerical software MATCONT [38] to investigate the stability of this equilibrium. Due to their importance in epidemic outbreaks, β and γ are chosen as bifurcation parameters.
The bifurcation analysis provides an integrated view of the dynamics of the model (1). It reveals that E * is stable and rabies persists for R 0 > 1 (see Figure 1). The figure shows several ranges of values for β and γ for a given r. The endemic equilibrium E * that is stable in region II bifurcates on the Hopf curve emitting a stable limit cycle through supercritical Hopf bifurcation, resulting in the formation of region III. As a result, the system exhibits periodic oscillations when the values of β and γ are chosen from this region. The size of the limit cycle that emerged from the Hopf bifurcation increases in the direction of increasing γ, which finally disappears via a saddle homoclinic bifurcation on the homoclinic curve and thus, region IV is born. As a result, the values of β and γ chosen from region IV lead to total extinction of the population, suggesting that the values of β and γ corresponding to this region are unreasonably large compared to r.

Controlling Rabies
The above analysis shows that the dynamics of model (1) are largely controlled by the parameter r in association with β and γ. Therefore, one can speculate that the spread of rabies can be hindered by controlling the birth rate, or by securing the newborns, for example through vaccination. We experiment on three different vaccination strategies: One is continuous vaccination, in which the susceptible population is vaccinated at a constant rate. The next one is pulse vaccination, which refers to vaccinating a fraction of the entire susceptible population in a single pulse applied every θ years. And, the last one is time dependent vaccination accompanied with carrying capacity control.

Continuous Vaccination Model
Here, we assume that dogs are vaccinated and neutered at the rate k and move to the class R. As there are vaccines that need booster dose after 2-3 years, while the street dog life span is only 3 years, we assume vaccine provides permanent immunity.
The basic reproduction number reduces to, From the expression for R 0 , we find that it increases with β and decreases with k. For R 0 < 1, we can easily find that, k > β σγm (µ+m)(m+σ) − m i.e., the threshold for k is a linear function of β. This shows that vaccination coverage should be high in the case of a high transmission rate.

Pulse Vaccination Model
Here, we present the pulse vaccination model derived from model (1). According to this strategy, a fraction k = kθ of the entire susceptible dog population is vaccinated in a single pulse every θ years [39], where the value of k is such that 0 < kθ < 1.

Optimal Control Model
In this section, we intend to explore strategies apart from vaccination that can be used independently or in combination with vaccination to reduce virus transmission among the dog population effectively. We consider vaccination and moderated growth rate through waste management as potentially effective two control measures. As we intend to limit the growth rate by reducing the carrying capacity, we aim to enhance the second-order term in the logistic growth model, since a decrease in carrying capacity can alternatively be achieved by increasing this term. We assume u 1 can increase the second-order term ρ times at maximum, where u 1 ∈ [0, 1] and ρ are constants, that depends on how much stress can be applied to suppress growth. In our simulation, we assume ρ = 10, which means the impact of limiting of resources (second order term in the logistic growth) can be increased a maximum of 10 times by implementing carrying capacity control.
We assume the vaccine is provided at the rate u 2 , where u 2 ∈ [0, 1]. The model including the above two controls reduces to the following.
Though our objective is to reduce rabid dogs, it is neither possible to reduce and maintain the growth rate at zero nor to vaccinate all dogs. So, we construct a cost functional as follows to balance between the two competing aims of reducing rabid dogs and implementing control measures.
Here, T is the final time; C is a constant representing the cost associated with each infected dog and plausible human infections resulted from the infected dog; B 1 and B 2 are the costs associated with the implementation of the controls u 1 and u 2 , respectively. Our aim is to identify the time dependent piece-wise continuous controls u * 1 (t) and u * 2 (t) that minimize the cost J(u 1 (t), u 2 (t)) subject to the system (8). Two time-dependent controls govern the dynamics of the system so as to minimize the numbers of infected dogs as well as the cost in order to implement the controls.
Mathematically speaking, the aim is to determine the optimal controls u * 1 (t) and u * 2 (t) such that Such an optimal control pair (u * 1 (t), u * 2 (t)) exists (Please see the Appendix B for details) and the necessary conditions to be satisfied by the optimal control are derived from Pontryagin's Maximum Principle [40]. We construct the Hamiltonian, H, given by, Here, λ 1 , λ 2 , λ 3 , and λ 4 are called adjoint variables and given by the following system of ordinary differential equations.
with transversality conditions, Finally, the optimal controls u * 1 (t), u * 2 (t) are characterized by ∂H ∂u 1 = 0 and ∂H ∂u 2 = 0, which gives In order to numerically solve the optimal control system, we use the Forward-Backward Sweep method [41]. We assume C = 1, B 1 = B 2 = 10. We assume B 1 = B 2 to understand the relative importance of controlling carrying capacity and vaccination. The controls cannot be 100% effective, i.e., as we can neither reduce the growth rate to zero through waste food management, nor vaccinate all of the dogs at once. We assume that u 1 does not exceed 0.8. u 2 refers to the continuous vaccination rate, which has been reported to be 0.09 in [31] and 0.24 in [42]. As this study deals with street dogs, we assume it does not exceed 0.3. Here, we present the details of a simulation to observe the optimal control strategies for a period of 5 years.

Parameterization
The population dynamics of street dogs is very different from that of domestic well-supervised dogs. Yearly net per-capita growth is assumed to be 0.0027 × 365 − 0.33 = 0.6555 [16]. The life-span of a street dog is about 3 years [43], which is much shorter than that of a domestic dog [31]. So we assume m = 0.33 per year. We assume a 2 month incubation period [31], i.e., σ = 6. The probability of clinical outcome γ = 0.4. Infectious lifetime is assumed to be 5.7 days [43], which gives a disease related death rate µ = 1/(5.7/365) ≈ 64. Street dog rabies related parameters are described in Table 1. We now infer the carrying capacity, D, and transmission rate, β, as well as R 0 in Dhaka. From the data, we observe that the dog population almost doubled just in 3 years from 2013 to 2016, i.e., according to our model 1 < R 0 < 1.249. We are now especially interested in which value of R 0 and D lead to that particular growth over 3 years. Assuming the population was at equilibrium in 2013, we simulated our model for 1 < R 0 < 1.249, with different values of D. It should be noted that the initial conditions in 2013 were adapted according to the CNVR campaign outcome before the simulation [21]. The whole process is summarized in the following.
Step 1: Find all possible values of endemic equilibrium for N * = 18585 and 1 < R 0 < 1.249 to be used as initial condition.
Step 3: Run the model for each R 0 and associated initial condition found in Step 2 to find the value of D that produces 37, 009 dogs in 2016 [24].
The above algorithm gives us plausible R 0 − D pairs to explain the rabies epidemic and dog population dynamics in Dhaka from 2013 to 2016, which is portrayed in Figure 2.

Results and Discussion
The model considered in this study is concerned with rabies transmission among FRSDs. The analyses show that a rabies outbreak is not possible for R 0 < 1. For R 0 > 1, an epidemic takes off. We consider reducing the dog carrying capacity as a means of reproduction control in addition to vaccination as possible effective control measures. The effect of carrying capacity on the dog population has already been established [7,15]. Reducing shelters and food sources through waste management could be a cheap way to reduce dog turn over [7,15,44]. Here, we examine theoretically whether this method is effective in controlling rabies transmission.
For each R 0 − D pair shown in Figure 2 we find a distribution of dog population for four different compartments for the year 2016. In further analysis, we use this as the initial condition and simulate model from 2016. It is worth mentioning here that there was another CNVR program in 2016, in which 26% of the dogs in the North and 13% of dogs in the South of Dhaka were vaccinated [2,17]. We update the initial conditions considering the vaccination that was administered in 2016. We then used these updated initial conditions to simulate the model with constant vaccination, the pulse vaccinating model, and the model that introduce optimal vaccination along with carrying capacity control. We simulated these models for a period of 5 years, where the value of R 0 was chosen as 1.1 and 1.23 along with two types of vaccination coverage (0 < u 2 < 0.1 and 0 < u 2 < 0.3). In order to make a comparison, we simulated the model for no vaccination as well. These results are presented in Figure 3A-D. The red dotted and solid lines in Figure 3 show the number of infected dog population that is obtained from the constant vaccination model with and without vaccination. The blue lines show the number of infected dog population for pulse vaccination at three different rates, and the green lines show this dog population with optimal controls (viz. carrying capacity control, u 1 and vaccination, u 2 ). We simulated the optimal control model under three different scenarios: when both u 1 and u 2 are implemented, when only u 1 is implemented, and finally when only u 2 is implemented. The corresponding values of the time dependent controls are shown on the right panel. In our simulations, we refer to R 0 = 1.1 as low transmission setting and to R 0 = 1.23 as high transmission setting. Similarly, we refer to k = 0.1 as low vaccination rate and to k = 0.3 as high vaccination rate.
We notice that the implications of different controls used in this study may provide a significantly better dog rabies scenario compared to no control at all; however, when we look further into the details, we see that one produces better results in reducing the number of rabid dogs than the other. Note that their performance also depends on the transmission settings as well as on the vaccination rates.
We find that both the constant and pulse vaccination corresponding to low vaccination rate reduce the number of rabid dogs greatly; however, the optimal control model where the vaccination is accompanied by carrying capacity control produces even a better result. This is true both in low and high transmission settings (see Figure 3A,B). It is worth mentioning that the carrying control alone (when only u 1 is implemented) is shown to be more promising compared to constant and pulse vaccinations in high transmission settings, as opposed to low transmission. This indicates that the carrying capacity control alone may not be a reliable strategy to control rabies in general, and signifies the simultaneous implementation of vaccination and carrying capacity control. In contrast, the pulse vaccination is shown to be more promising compared to other control strategies in the case of a high vaccination rate (see Figure 3C,D). It is true both in low and high transmission settings; however, it should be noted that vaccinating FRSDs at the rate of 30% every year is obviously a quite challenging task.

Conclusions
We studied a dog population model of rabies transmission dynamics to investigate the role of vaccination and dog population management in controlling rabies. The analysis of the basic model showed that the birth rate of dogs significantly affects the dog dynamics and rabies transmission in FRSDs. This quantity is directly associated with an upper limit of the transmission coefficient, in other words an upper limit of the basic reproduction number at which the dog population gets extinct; therefore, it was primarily expected that reducing this quantity may help to control or reduce the spread of rabies. The idea was validated by the results of the advanced model, where we found that carrying capacity reduction significantly reduces the number of rabid dogs; however, this control measure alone was not effective enough to control rabies, and it showed different performances under different transmission settings and vaccination rates as follows.
Corresponding to a low vaccination rate, the vaccination was found to dominate other control measures in low transmission setting, whereas the carrying capacity reduction was found to play a dominant role in the case of high transmission setting; therefore, the carrying capacity control was expected to be a good supplementary control measure. It leads us to conclude that simultaneous implementation of vaccination and carrying capacity control could be a reliable approach to control rabies in FRSDs. This notion is supported by the results corresponding to the low vaccination rate, as shown in Figure 3A,B. In contrast, an alternative scenario was observed in the case of a high vaccination rate, where the pulse vaccination was found to dominate other control measures both in low and high transmission settings; however, vaccination at the proposed rate of 30% every year is obviously a quite challenging task for FRSDs [17], as the maximum vaccination rate in Dhaka was 26% in 2016 [2]. Therefore, we recommend the simultaneous use of carrying capacity control and vaccination as control measures.
In conclusion, the results of this study demonstrate that vaccines alone should not be considered as the sole instrument to control rabies in FRSDs. The carrying capacity control obtained through food waste management should also be considered as an additional control measure. Food waste control requires mass education and awareness, and may not be equally effective in all communities as some FRSDs might have an alternative food source [6,45,46]. Community specific studies considering specific cultural practices towards FRSDs and available food sources may help to identify more appropriate control measures. The results of this study should encourage similar efforts to raise community awareness against dumping food waste in the undesignated area, leaving food garbage opens for dogs to access or feeding street dogs. It is also important to implement the suggested guidelines as a means to realize appropriate control measures for communities. Finally, the model parameters were estimated on a trial basis due to a shortage of data. These parameters could be estimated more accurately when there are sufficient data and hence, leading to better accuracy of the predicted results.  and
Proof. First, the system (8) is bounded in the positive orthant and for each control pairs there exists solutions for the corresponding state variable. Thus the set of controls and the state variables is nonempty. By construction, the control set U is convex and closed. Next, the integrand of the objective functional is linear in state variables and quadratic in controls u 1 , u 2 , and thus is convex. Further, we choose c 1 = min{B 1 , B 2 }, c 2 2 = min{C 1 E(t), C 2 I(t)} and ρ = 2 then, the integrand of the objective functional satisfies the following.
These assures the existence of optimal control and completes the proof [40].