Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 1, Number 4, December 2020, pp.355-372 https://doi.org/10.5206/mase/10847 CONTROLLING RABIES EPIDEMICS IN NEPAL WITH LIMITED RESOURCES: OPTIMAL CONTROL THEORY APPROACH BUDDHI PANTHA, HEM RAJ JOSHI, AND NAVEEN K. VAIDYA Abstract. In many developing countries, including Nepal, rabies epidemics constitute a serious public health concern, partly because of limited resources for proper implementation of control measures. In this study, we develop an extended model by incorporating various controls into the transmission dynamics model with both dog and jackal vectors. We apply the optimal control theory on the developed model system to identify optimal control strategy for mitigating rabies burden in Nepal with limited resources. Among the potential control strategies, human vaccination, dog vaccination, dog culling, dog sterilization, and jackal vaccination, considered in this study, our results show that a combination of dog vaccination and dog culling is the most effective strategy to control rabies in Nepal. Our optimal control solutions provide strategies for optimal implementation of these controls to suppress rabies prevalence in humans, dogs and jackals of Nepal using the minimum cost associated with controls. We found that given limited resources, implementing controls in a time-dependent manner with a higher level at the beginning of the outbreaks and reducing them during later part of the epidemics can provide maximum benefits. 1. Introduction Rabies, a viral zoonotic disease, remains an ongoing burden in many developing countries, including Nepal. Because of extremely high fatality rate (almost 100%) in rabid humans or animals having symp- toms such as violent movements, uncontrolled excitement, fear of water (hydrophobia), an inability to move parts of the body, confusion, and loss of consciousness [2, 6, 27, 31], this disease poses extreme threats of public health concerns. While multiple control strategies are available, including successful vaccine, limited resources and lack of proper allocation of resources often make developing countries fail to control rabies epidemics; it is also called a poor man’s disease because most of deaths (≥ 95%) occur in Asian and African countries [33]. Therefore, well-designed planning is necessary before implementing these control strategies to achieve optimal outcomes using limited resources. Mathematical modeling can provide an important means for identifying ideal planning strategies. Rabies virus is mainly transmitted through infected animal bites [8]. Among the potential animals, dogs are primarily considered as vector for the transmission of rabies. There are many existing models that describe transmission dynamics of the rabies with dogs as primary vectors [1, 2, 3, 7, 8, 9, 11, 12, 17, 18, 20, 22, 27, 30, 32]. However, in the context of some countries like Nepal, while dogs remain primary vectors in urban epidemiological cycle [26], jackals as secondary vectors in sylvatic epidemiological cycle Received by the editors 30 June 2020; accepted 1 December 2020; published online 10 December 2020. 2010 Mathematics Subject Classification. Primary 35K55; 92D30; Secondary 46E25, 20C20. Key words and phrases. dog and jackal vectors; mathematical model; optimal control theory; rabies in Nepal; resource limitation. Authors would like to thank Association of Nepalese Mathematicians in America (ANMA) for organizing a Workshop on Collaborative Research in Mathematical Sciences (May 25- 27, 2018), during which this work was initiated. Vaidya’s work was supported by NSF grants DMS-1951793, DMS-1616299, DMS-1836647, and DEB-2030479 from National Science Foundation of USA, and UGP award from San Diego State University. 355 356 B. PANTHA, H. JOSHI, AND N. VAIDYA also play an important role in the persistent ongoing rabies epidemics [21]. Such persistent low level of rabies infection in jackals were also found in Zimbabwe [24]. In a recent study, we developed a model that couple both dog and jackal along with human population to describe transmission dynamics of rabies in Nepal, and identified that consideration of both vectors is essential for successful mitigation of rabies in Nepal [21]. In particular, our model predicted that even though intraspecies transmission is prevented among each animal species (dogs and jackals), the rabies can still persist due to interspecies transmission [21]. Nepal has pledged to end rabies by 2030, but dog-bite cases are rising and as many as 26,312 peo- ple were administered with post-exposure anti-rabies vaccines at government health facilities in 2018 [14]. Despite continuous effort to control rabies, ongoing epidemic can partly be attributable to limited resources for implementing proper control and existence of secondary vectors. It is thus important to identify optimal strategy to implement control programs regarding underlying situation of two different vectors and limited resources. For such purposes, optimal control theory has been proved to be use- ful tool as in many previous epidemic controls, including rabies epidemics [5, 7, 13, 15, 16, 19, 20, 23, 25]. In this study, we introduce effects of various controls into the transmission dynamics model incor- porating both dog and jackal populations [21]. In particular, we focus on human vaccination, dog vaccination, dog culling, dog sterilization, and jackal vaccination through bait as well as combinations of them. Using the developed dynamical system model, we further formulate an optimal control problem to take resource limitation into account. Implementing optimal control theory and related numerical method, we compute the optimal control strategy for successful control of rabies epidemic in Nepal. 2. Mathematical Model with Controls Based on our previous model of rabies transmission dynamics in Nepal [21], we develop an extended model by incorporating various control measures to describe the dynamics of rabies transmission in multiple groups of species (jackals, dogs, and humans). Specifically, our model is a coupled system of differential equations that describe the rate of change of subpopulations of jackals (J), dogs (D), and humans (H) under various control programs implemented in the community. We use subscripts J, D and H in the variables and parameters to represent them corresponding to jackals, dogs, and hu- mans, respectively. The total population of each species is divided into four subpopulations: susceptible (SJ,SD,SH), exposed (EJ,ED,EH), vaccinated (VJ,VD,VH), and infected (IJ,ID,IH). The population in each group is recruited with rate Λ into susceptible class and die with natural death rate µ. Susceptible population get infected and enter into exposed class with rate β. Also, the exposed humans or animals transit into infected class with rate γ and the infected populations die due to rabies with rate δ. As observed in the context of Nepal, note that humans are infected through dog- bites only, while both dogs and jackals are infected through intra-species and inter-species transmissions. There are various potential prevention and control strategies that can be applied to break the Jackal- dog-human transmission sequence for rabies. In this study, we consider the five most frequently used strategies: human vaccination (post-exposure), dog vaccination (both pre and post-exposure), dog culling, dog sterilization, and jackal vaccination through bait. We use u1 to denote the rate at which the exposed humans get the rabies vaccine. The vaccination program for the rabies in dogs includes both pre- and post-exposure vaccination. As per guidelines in [4], if vaccinated dogs are exposed, revaccina- tion should be administered immediately. Therefore, we apply the dog vaccination to both susceptible CONTROLLING RABIES IN NEPAL 357 and exposed dogs at the same rate of u2, and we assume that the immunity is not lost for the period of dynamics considered in the study. The third rabies control strategy is dog culling, which we denote using the rate u3 implemented to cull dogs from all classes, regardless of their infection status. The bait vaccination for jackals is applied using foods containing the rabies vaccine that is spread in different locations so that the jackals consume the foods, reducing the rabies contraction [6, 28]. We denote the rate of jackal vaccination for both susceptible and exposed jackals by u4. The dog sterilization strategy is used to control reproductive rate of dogs, eventually reducing the recruitment rate. We denote the net effectiveness of dog sterilization by u5 so that the dog recruitment rate changes to (1 −u5)ΛD. As described above, the transmission dynamics of rabies under these five control strategies can be represented using the following system of differential equations. S′J = ΛJ −βJJSJIJ −βDJSJID − (µJ + u4)SJ E′J = βJJSJIJ + βDJSJID − (γJ + µJ + u4)EJ V ′J = u4(SJ + EJ) −µJVJ I′J = γJEJ − (µJ + δJ)IJ S′D = (1 −u5)ΛD −βJDSDIJ −βDDSDID − (µD + u2 + u3)SD E′D = βJDSDIJ + βDDSDID − (γD + µD + u2 + u3)ED V ′D = u2(ED + SD) − (µD + u3)VD (2.1) I′D = γDED − (µD + δD + u3)ID S′H = ΛH −βDHIDSH −µHSH E′H = βDHIDSH − (u1 + µH + γH)EH V ′H = u1EH −µHVH I′H = γHEH − (µH + δH)IH The model parameters related to the context of Nepal [21] are given in Table 1. 3. Impact of Controls on Rabies Epidemics in Nepal To identify the most impactful controls, we first evaluate the effects on constant level of control on preventing rabies epidemic and/or reducing rabies prevalence. In our previous study [21], we analyzed the impact of implementing the single control at a time, and found that the use of only one control is not able to mitigate the disease unless the level of control is significantly high. For example, rabies prevalence in dog and jackal can be reduced to low level only if the annual culling for the dog is more than 40% effective for 10 years or the effectiveness of dog vaccination is more than 60%. The requirement of unusually high level of control for a longer period of time implies that the multiple control strategies need to be applied simultaneously for the successful control of rabies. On the other hand, it is unlikely for developing countries like Nepal to implement many control programs at the same time due to resource limitation. Therefore, we mainly focus on combinations of two control programs. Here, we consider combinations of two different controls and evaluate which control combinations are more effective in reducing the basic reproduction number, R0, as well as reducing the rabies prevalence among dog and jackal populations. The most effective combination of two control measures are then considered to identify the optimal planning for implementing them in resource limited setting. 358 B. PANTHA, H. JOSHI, AND N. VAIDYA Table 1. The model variables and parameters. Model Variables and Initial Values Variable Description Initial values SJ Susceptible Jackals 73125 EJ Exposed Jackals 368 VJ Vaccinated Jackals 0 IJ Infected Jackals 73 SD Susceptible Dogs 15.898 × 105 ED Exposed Dogs 10 4 VD Vaccinated Dogs 4 × 105 ID Infected Dogs 200 SH Susceptible 25.265 × 106 EH Exposed 15534 IH Infected 1000 VH Vaccinated 14000 Model Parameters Params. Description Value µJ Jackal mortality rate 0.125 µD Dog mortality rate 0.2 µH Human mortality rate 0.0142 ΛJ Jackal recruitment rate µJ ×NJ (0) ΛD Dog recruitment rate µD ×ND(0) ΛH Human recruitment rate µH ×NH (0) δJ Jackal rabies related mortality rate 36.5 δD Dog rabies related mortality rate 36.5 δH Human rabies related mortality rate 36.5 γJ Jackal rate of moving from exposed to infected 6.64 γD Dog rate of moving from exposed to infected 2 γH Human rate of moving from exposed to infected 2 βJJ Transmission rate from Jackal to Jackal 3.79 × 10−5 βDJ Transmission rate from Dog to Jackal 1.90 × 10−5 βJD Transmission rate from Jackal to Dog 1.52 × 10−5 βDD Transmission rate from Dog to Dog 2.74 × 10−5 βDH Transmission rate from Dog to Human 1.71 × 10−6 uD Dog vaccination rate(susceptible and exposed) 0.03 uH Human vaccination rate (PEP) 2.05 3.1. Impact on the basic reproduction number. The basic reproduction number, R0, known as the expected number of secondary cases produced by a single (typical) infection in an entirely susceptible population [10], can be used to determine whether the outbreak occurs (R0 > 1) or infection dies out (R0 < 1) [10]. Applying the next generation matrix method [10] to our model, we consider the subsystem containing all of the equations except the equations for SH, SD and SJ. This subsystem is then linearized about the disease free equilibrium (DFE), given by( ΛJ µJ + u4 , 0, u4ΛJ µJ(µJ + u2) , 0, (1 −u5)ΛD (µD + u2 + u3) , 0, u2(1 −u5)ΛD (µD + u3)(µD + u2 + u3) , 0, ΛH µH , 0, 0, 0 ) . From the resulting equations, we obtain a matrix F containing infection terms and a matrix V containing transfer terms, as follows. F =   F11 F12 0F21 F22 0 0 F32 0   , and V =   V11 0 00 V22 0 0 0 V33   , CONTROLLING RABIES IN NEPAL 359 where F11 =   0 0 βJJ ΛJ µJ +u4 0 0 0 0 0 0   , F12 =   0 0 βDJ ΛJ µJ +u4 0 0 0 0 0 0   , F21 =   0 0 (1−u5)βJDΛD (µD+u2+u3) 0 0 0 0 0 0   , F22 =   0 0 (1−u5)βDDΛD (µD+u2+u3) 0 0 0 0 0 0   , F32 =   0 0 βDH ΛH µH 0 0 0 0 0 0   , V11 =   γJ + µJ + u4 0 0−u4 µJ 0 −γJ 0 µJ + δJ   , V22 =   γD + µD + u2 + u3 0 0−u2 µD + u3 0 −γD 0 µD + δD + u3   , and V33 =   u1 + µH + γH 0 0−u1 µH 0 −γH 0 µH + δH   . This implies FV−1 =   F11V−111 F12V−122 0F21V−111 F22V−122 0 0 F32V−122 0   . The basic reproduction number is then given by the spectral radius of the matrix FV−1. Therefore, R0 = ρ(FV−1) = 1 2 ( RJ0 + R D 0 + √ (RJ0 −RD0 )2 + 4CE ) , where RJ0 = βJJγJΛJ (µJ + u4)(δJ + µJ)(γJ + µJ + u4) , RD0 = (1 −u5)βDDγDΛD (δD + µD + u3)(µD + u2 + u3)(γD + µD + u2 + u3) , C = βDJγDΛJ (µJ + u4)(δD + µD + u3)(γD + µD + u2 + u3) , and (3.1) E = (1 −u5)βJDγJΛD (δJ + µJ)(γJ + µJ + u4)(µD + u2 + u3) . We now use the formula for R0 derived above and parameters given in Table 1 to compute the value of basic reproduction number for various control strategies with two controls taken at a time (Figure 1). Note that the human vaccination at its base level is always included in all strategies as it cannot be avoided in practice. We also identified the control levels for which R0 is less than 1, leading to eradication of rabies in Nepal. While each combination strategy has certain levels that can bring R0 below 1, the level required is quite high for some strategies compared to others. As presented in Figure 1 we observe that jackal vaccination is one of the least effective strategies on lowering the basic reproduction number; the increase in the coverage of jackal bait vaccination has negligible impact on the basic reproduction number (Figure 1a, 1c, 1e). Similarly, the dog sterilization does not show significant effectiveness either to reduce R0 (Figure 1a, 1b). The dog culling and dog vaccination on the other hand are highly impactful on bringing the value of R0 below 1 (Figure 1d). In particular, a combined 360 B. PANTHA, H. JOSHI, AND N. VAIDYA Figure 1. Heatmap of R0 for various combination of pair of controls. The dotted curve represents R0 = 1. strategy, including the dog vaccination and dog culling control measures, is the best approach to reduce R0, thereby preventing or eradicating the rabies in Nepal (Figure 1d). 3.2. Impact on long-term rabies prevalence. In this section, we present the effects of combined strategies with two control measures at a time on the long term prevalence of rabies among dogs and jackals. The parameter values and the initial values are used as given in Table 1. The Runge-Kutta method of fourth order is applied for the model simulations. Since there is no transmission of rabies from humans to animals, the human vaccination does not have any effect on the rabies prevalence in dog and jackal populations. Therefore, we keep the human vaccination at a constant level estimated using the data from Nepal [21] and focus on other four control measures. In the following simulations (Figure 2-7), we compute the prevalence of rabies in dog and jackal populations with the application of combined strategies for the period of 10 years. In our simulations, rate of each control strategy is estimated from the target percentage of the population that are to be covered by that strategy. For example, with the vaccination rate of u2, the target dog population, XD, can be approximated using the solution of dXD dt = −u2XD, i.e, XD(t) = XD(0)e−u2t. Then for a program that aims to vaccinate η2% of the dog population in a period of t years, the vaccination rate u2 needs to be in such a way that XD(t) = (1 −η2/100)XD(0) = XD(0)e−u2t, which gives u2 = − ln(1−η2/100) t . Similarly, the dog culling rate can be written as u3 = − ln(1−η3/100) t , where η3 is the % coverage for dog culling in a period of t years. First, we consider combined strategy with dog vaccination and dog culling control measures and com- pute the prevalence of rabies in dog and jackal populations. As expected, our model prediction shows that applying one control strategy at a lower level requires another control strategy at higher level to achieve sufficient decrease in the long-term prevalence of rabies among dogs and jackals. For example, annual dog vaccination coverage at a level of 5% for a period of 10 years requires more than 5% of CONTROLLING RABIES IN NEPAL 361 Figure 2. The rabies prevalence in (a) dog, and (b) jackal populations under various levels of annual coverage of dog vaccination and dog culling. annual dog culling for that period to lower the rabies prevalence among dogs below 0.006%, and similar dog vaccination and culling coverage is needed to achieve the same low level of rabies prevalence in jackal. Next, we consider dog vaccination and jackal bait vaccination together to compute the prevalence of rabies in dog and jackal populations. As shown in Figure 3, we observe that Jackal vaccination does not Figure 3. The rabies prevalence in (a) dog and (b) jackal populations under various level of annual coverage by dog vaccination and jackal vaccination. have significant impact in lowering the rabies prevalence in the dog population, i.e., the rabies preva- lence among dogs remains almost the same for any level of coverage by jackal vaccination. However, the rabies prevalence among jackals, is impacted by this combination of dog vaccination and jackal vaccination. For example, a combined program with 5% coverage of each of dog vaccination and jackal bait vaccination results in the rabies prevalence in jackal of about 0.033% while an increase in the jackal vaccination to 15% keeping the same 5% dog vaccination coverage decreases the rabies prevalence in jackal to 0.018%. The third strategy we consider is the combination of the dog vaccination and the dog sterilization. Note that sterilization does not control the rabies incidence directly, it rather controls the disease in a 362 B. PANTHA, H. JOSHI, AND N. VAIDYA Figure 4. The rabies prevalence among (a) dog and (b) jackal populations under various level of annual coverage by dog vaccination and dog sterilization. long run by decreasing the reproduction of the dog and eventually reducing the susceptible populations of dog. Therefore, for fare comparison as in our previous work [21], we measure the sterilization strength in a 5 year time frame as opposed to other strategies which are measured in 1 year time frame. In this strategy, with dog sterilization strength of 15% coverage over 5 years, at least 10% of annual coverage of dog vaccination is required for 10 years to keep the prevalence in dog population below 0.003% (Figure 4), while at least 11% dog vaccination coverage is required for 10 years to keep the rabies prevalence in jackal population below 0.003%. We observe that dog sterilization does not have significant impact in lowering the rabies prevalence in dog and jackal populations. The fourth strategy is the combination of the dog sterilization and the jackal vaccination. This strat- Figure 5. The rabies prevalence among (a) dog and (b) jackal populations under various level of annual coverage by jackal vaccination and dog sterilization. egy does not have significant impact on the rabies prevalence in dog population (Figure 5). However, it has considerable impact on the rabies prevalence in jackal population. For instance, dog sterilization level of 10% and jackal vaccination coverage of 15% for 10 years can bring the rabies prevalence below 0.050% in jackal population. For a higher dog sterilization coverage (for example, 17% for 10 years), the same level of rabies prevalence can be achieved with lower jackal vaccination level of 4%. CONTROLLING RABIES IN NEPAL 363 The fifth strategy considered in this study is the combination of dog sterilization and culling. In this Figure 6. The rabies prevalence among (a) dog and (b) jackal populations under various level of annual coverage by jackal vaccination and dog culling. combination of strategies, the dog culling has significant impact in reducing rabies prevalence in both dog and jackal populations (Figure 6a), but the dog sterilization has only a little impact in reducing the rabies prevalence in both populations (Figure 6b). For example, at 5% of dog culling for 10 years, a change of the level of dog sterilization from 5 to 15% causes the rabies prevalence in dog population to change from 0.0351% to 0.0122%. In this change, the prevalence in jackal populations also changes from 0.0289% to 0.0095%. The last strategy considered in this study is the combination of the jackal vaccination and the dog culling. As in the previous combination of strategies, dog culling has significant impact on reducing Figure 7. The rabies prevalence among (a) dog and (b) jackal populations under various level of annual coverage by jackal vaccination and dog culling. rabies prevalence in both populations but the jackal vaccination plays a noticeable role in reducing rabies prevalence in jackal population only. For example, at 5% of dog culling for 10 years, the rabies prevalence in jackal population changes from 0.45% to 0.011% on changing the level of jackal vaccina- tion from 0 to 25%, but the prevalence in dog populations does not change significantly (0.055% for 0% 364 B. PANTHA, H. JOSHI, AND N. VAIDYA coverage and 0.047% for 25% coverage of jackal vaccination). For comparison purpose, we consider the rabies prevalence among dogs and jackal populations under each combined strategy with 10% coverage level of each of two strategies included in the combination. For this level of (10%, 10%) coverage, our results show that the rabies prevalence for dog population remains 0.0001%, 0.0040%, 0.0036%, 0.1362%, 0.0023% and 0.0056% for combined strategies with dog vaccination and dog culling, jackal vaccination and dog vaccination, dog vaccination and dog steriliza- tion, dog sterilization and jackal vaccination, dog culling and dog sterilization, and jackal vaccination and dog culling, respectively. In the jackal population the corresponding prevalence is 0.0001%, 0.0026%, 0.0041%, 0.0652%, 0.0027% and 0.0025%, respectively. Therefore, dog vaccination and dog culling are the most effective measures to prevent and control rabies epidemics in Nepal. In the following section, we focus on the optimal strategy to implement the combination of dog vaccination and dog sterilization under resource limited settings. 4. Optimal Control of Rabies in Nepal In this section, we use our model to formulate optimal control problem, which allows us to identify optimal time-dependent strategies under limited resources to achieve maximum benefit from the control strategy implementation. While we formulate general optimal control problem including all possible control strategies, we particularly emphasize on two most effective strategies, namely the dog vaccination and the dog culling, identified in Section 3. 4.1. Formulation of optimal control problem. We consider controls as time-dependent functions, i.e, ui = ui(t) for i = 1, 2, ..., 5. To incorporate resource limitation while controlling the rabies burden, we set a goal of minimizing the number of exposed and infected humans and animals as well as min- imizing the anticipated cost for control implementation for a fixed period of time, say tf . To achieve this goal, we formulate an objective functional as follows. J(u1,u2,u3,u4,u5) = min (u1,u2,u3,u4,u5)∈U ∫ tf 0 [A1(EH + IH) + A2(ED + ID) + A3(EJ + IJ) + B1u1EH + B2u2(SD + ED) + B3u3ND + B4u4NJ + B5u5(ND − ID) +C1u 2 1 + C2u 2 2 + C3u 2 3 + C4u 2 4 + C5u 2 5 ] dt, (4.1) where Ai’s, Bi’s, and Ci’s are the positive constants, associated with weights corresponding to disease outcome and costs. Here, we assume that the admissible control set U is given by U ={(u1,u2,u3,u4,u5) ∈ R5 : 0 ≤ ui(t) ≤ bi for i = 1 . . . 5 and ui are Lebesgue measurable}, where bi’s are positive constants related to the availability of resources. We consider the formulated optimal control under dynamical system given by Eqs. (2.1). For the control problem formulated above, we can apply a result from Lukes [29] to prove the existence and uniqueness of solutions for the state system (2.1) with the given controls. The existence and uniqueness results for our optimal control problem can be summarized as in Theorem 4.1. Theorem 4.1. Given controls u = (u1,u2,u3,u4,u5) ∈ U, there exist non-negative bounded solutions (SJ,EJ,VJ,IJ,SD,ED,VD,ID,SH,EH,VH,IH) to the state system (2.1) in the finite interval [0,T] with given initial conditions. CONTROLLING RABIES IN NEPAL 365 The structure of system (2.1) gives the non-negativity and uniform boundedness of the state solutions. As stated in Theorem 4.2, we can also assert the existence of the optimal controls based on the results from [16]. Theorem 4.2. There exists an optimal control tuple u∗ = (u∗1,u ∗ 2,u ∗ 3,u ∗ 4,u ∗ 5) ∈ U with corresponding states (S∗J,E ∗ J,V ∗ J ,I ∗ J,S ∗ D,E ∗ D,V ∗ D,I ∗ D,S ∗ H,E ∗ H,V ∗ H,I ∗ H) that minimizes the objective functional J(u1,u2,u3,u4,u5). By using Pontryagin’s Maximum Principle as stated in [16, 23], we are able to derive necessary conditions for our optimal control and corresponding states. The Hamiltonian of the system is: H = A1(EH + IH) + A2(ED + ID) + A3(EJ + IJ) + B1u1EH + B2u2(SD + ED) + B3u3ND + B4u4NJ + B5u5(SD + ED + VD) + C1u 2 1 + C2u 2 2 + C3u 2 3 + C4u 2 4 + C5u 2 5 + λ1(ΛJ −βJJSJIJ −βDJSJID − (µJ + u4)SJ) + λ2(βJJSJIJ + βDJSJID − (γJ + µJ + u4)EJ) + λ3(u4(SJ + EJ) −µJVJ) + λ4(γJEJ − (µJ + δJ)IJ) + λ5((1 −u5)ΛD −βJDSDIJ −βDDSDID − (µD + u2 + u3)SD) + λ6(βJDSDIJ + βDDSDID − (γD + µD + u2 + u3)ED) + λ7(u2(ED + SD) − (µD + u3)VD) + λ8(γDED − (µD + δD + u3)ID) + λ9(ΛH −βDHIDSH −µHSH) + λ10(βDHIDSH − (u1 + µH + γH)EH) + λ11(u1EH −µHVH) + λ12(γHEH − (µH + δH)IH). For given optimal controls u1, ...,u5, there exist λi, i = 1, . . . , 12, with derivative λ ′ i(t) given by λ′i(t) = − ∂H ∂ (ith state variable) . Therefore, we obtain λ′1(t) = −(B4u4 −λ1(βJJIJ + βDJID + (µJ + u4)) + λ2(βJJIJ + βDJID) + λ3u4) , λ′2(t) = −(A3 + B4u4 −λ2(γJ + µJ + u4) + λ3u4 + λ4γJ) , λ′3(t) = −(B4u4 −λ3µJ) , λ′4(t) = −(A3 + B4u4 −λ1βJJSJ + λ2βJJSJ −λ4(µJ + δJ) −λ5βJDSD + λ6βJDSD) , λ′5(t) = −(B2u2 + B3u3 + B5u5 −λ5(βJDIJ + βDDID + µD + u2 + u3) , +λ6(βJDIJ + βDDID) + λ7u2) , λ′6(t) = −(A2 + B3u3 + B5u5 −λ6(γD + µD + u2 + u3) + λ7u2 + λ8γD) , λ′7(t) = −(B3u3 + B5u5 −λ7(µD + u3), ) λ′8(t) = −(A2 + B3u3 −λ1βDJSJ + λ2βDJSJ −λ5βDDSD , (4.2) +λ6βDDSD −λ8(µD + δD + u3) −λ9βDHSH + λ10βDHSH) , λ′9(t) = −(−λ9(βDHID + µH) + λ10βDHID) , λ′10(t) = −(A1 + B1u1 −λ10(u1 + µH + γH) + λ11u1 + λ12γH) , λ′11(t) = −(−λ11µH) , λ′12(t) = −(A1 −λ12(µH + δH)) , and the transversality conditions are λi(tf ) = 0, i = 1, . . . 12. 366 B. PANTHA, H. JOSHI, AND N. VAIDYA The optimal control solutions u∗i , i = 1, 2, ..., 5, can then be obtained by setting the derivative of Hamiltonian system with respect to each control to zero, i.e., ∂H/∂ui = 0, i = 1, 2, ..., 5, where ∂H ∂u1 = B1EH + 2C1u1 −λ10EH + λ11EH, ∂H ∂u2 = B2(SD + ED) + 2C2u2 −λ5SD −λ6ED + λ7(ED + SD), ∂H ∂u3 = B3(SD + ED + VD + ID) + 2C3u3 −λ5SD −λ6ED −λ7VD −λ8ID, ∂H ∂u4 = B4(SJ + EJ + VJ + IJ) + 2C4u4 −λ1SJ −λ2EJ + λ3(SJ + EJ), ∂H ∂u5 = B5(SD + ED + VD) + 2C5u5 −λ5ΛD, We obtain the optimal control solution as follows. u∗1 = min [ b1, max [ a1,− 1 2C1 (B1 −λ10 + λ11) EH ]] , u∗2 = min [ b2, max [ a2,− 1 2C2 (B2(SD + ED) −λ5SD −λ6ED + λ7(ED + SD)) ]] , u∗3 = min [ b3, max [ a3,− 1 2C3 (B3(SD + ED + VD + ID) −λ5SD −λ6ED −λ7VD −λ8ID) ]] , u∗4 = min [ b4, max [ a4,− 1 2C4 (B4(SJ + EJ + VJ + IJ) −λ1SJ −λ2EJ + λ3(SJ + EJ)) ]] , u∗5 = min [ b5, max [ a5,− 1 2C5 (B5(SD + ED + VD) −λ5ΛD) ]] . 4.2. Estimation of weight parameters and bounds for controls. It is important to determine the reasonable weight parameters Ai,Bi,Ci introduced into the objective functional, since the outcome of the minimization procedure may highly be impacted by the choice of these weights. Here, we follow a similar technique used in Mallela et al. [19] to make proper choice of these weight constants. To estimate the reasonable proportion of weights, we take A1 = 1, and compute other weight constants in such a way that the term corresponding to each weight is approximately the same as the term corresponding to A1. For example, we estimate A2 using∫ tf 0 A1(EH + IH)dt = ∫ tf 0 A2(ED + ID)dt, and obtain A2 = 1.0702. The similar technique allows us to obtain A3 = 95.8453, B1 = 0.5144, B2 = 0.0614, and B3 = 0.1259. To estimate Ci’s we use the average value, u av i , of minimum and maximum values of controls. The minimum value for both dog vaccination and culling rates are assumed to be zero for no vaccination and culling. For the upper bound, we assume that the available resource for dog vaccination corresponds to the maximum capacity of covering 40% of dogs in Nepal in a year. Using dXD dt = −u2XD, where XD is dog population remained to be vaccinated, we can compute that the maximum resource (i.e., 40% coverage) is equivalent to the dog vaccination rate of u2 = 0.5. Thus, we take the bounds for dog vaccination rate as 0 ≤ u2 ≤ 0.5. Next, as capturing the dogs, culling them and disposing them need more resources, culling process needs more manpower than vaccination process, and thus lower coverage for culling can be achieved with limited resources. Let us assume that the maximum available CONTROLLING RABIES IN NEPAL 367 resources for dog culling can cover 20% of dogs per year. With this assumption, the bounds for the dog culling rate is 0 ≤ u3 ≤ 0.22. Then we use A1 ∫ tf 0 (EH + IH)dt = Ci ∫ tf 0 (uavi ) 2dt, i = 1, 2, 3. to estimate C′is. The computation from our model solution provides C1 = 4.7562×10 4, C2 = 1.5991× 106, C3 = 8.2596 × 105. Since we intend to obtain optimal time-varying strategy with combined two most effective controls, u2(t) and u3(t), with underlying constant human post-infection vaccination (u1), the remaining two controls u4 and u5 are taken to be 0. As a result, the corresponding weight parameters vanish, i.e., B4 = B5 = C4 = C5 = 0. 4.3. Method for numerical computation. In this section, we briefly summarize the computational method used to obtain the optimal control solutions. Our technique is similar to the iterative algorithm introduced by Lenhart and Workman [16]. In particular, we use a backward-forward sweep iterative method with a fourth order Runge-Kutta scheme. Starting with initial guesses for the controls, the state equations (2.1) are solved forward in time. Then, the resulting state values are used to solve the adjoint equations (4.2) backward in time. The controls are then updated. This iteration process is repeated until the convergence is achieved. The convergence of iteration is defined as a condition, at which the value of variables in two successive iterations are negligibly close (i.e., their difference is smaller than a desired small number). The algorithm implemented in our study can be summarized as follows: Step 1: Input initial guess for controls over the interval [0, tf ]. Step 2: Using the initial values of state variables and the values of the controls, solve the system (2.1) forward in time, i.e., from t = 0 to t = tf . Step 3: Using the transversality conditions λi(tf ) = 0, i = 1, 2, ..., 12 as well as the values of the state variables from step 2 and the values of controls, solve the adjoint system (4.2) backward in time, i.e., from t = tf to t = 0. Step 4: Update the controls with the characterization u∗i , i = 1, 2, ...5 using the values of state and adjoint variables. Step 5: If the convergence is achieved, output the current values of solution, otherwise return to Step 2. 4.4. Solutions for the optimal control of rabies in Nepal. As discussed in Section (3), as well as in our previous work [21], the dog vaccination and the dog culling are two most effective intervention strategies to control the rabies epidemics in Nepal. Here, we present the numerical solutions of how these strategies can be implemented optimally given limited resources in the context of Nepal. First, we consider programs with only one of these two strategies and then consider a program with these two strategies combined. As mentioned earlier, note that there is always underlying post-infection human vaccination as this can not be avoided in practice. 4.4.1. Control program with dog vaccination only. We used the bounds for the vaccination rate 0 ≤ u2 ≤ 0.5 and the weight parameters from Section 4.2. All other model parameters are taken from Table 1. Our model simulations at the boundary levels of the dog vaccination control (Figure 8) show that without dog vaccination (u2 = 0) the prevalence of rabies in dog population increases continuously and reaches at a high level of 0.2543% as early as the ninth year, and in jackal population, the rabies prevalence increases and reaches 0.2162% in the eighth year. In human population, the rabies case increases continuously and reaches about 169,000 as soon as eighth year. On the other hand, with the 368 B. PANTHA, H. JOSHI, AND N. VAIDYA Figure 8. The rabies prevalence among (a) dog population, (b) jackal population, (c) num- ber of human rabies cases and (d) optimal dog vaccination profile, u2, in the control program with dog vaccination only. highest level of dog vaccination (u2 = 0.5), i.e., use of maximum resources, the prevalence of rabies in both dog and jackal populations reach approximately zero as early as fourth year and the number of rabies in humans is about one. Our optimal control solution indicates that the optimal strategy should be with the dog vaccination control profile u2(t) starting at rate 0.18 and tapering down to zero as shown in Figure 8d. With this profile, the value of objective functional can be brought to 36% less than without control and 64% less than with highest control (i.e., the value of J is 1.5310×106 without control, 6.2895 × 106 with highest control, and 2.2690 × 106 with optimal control), thereby utilizing the resources optimally. With this optimal vaccination strategy, the rabies prevalence in dog and jackal populations decreases slowly and maintains at about 0.0026%, as opposed to 0.2458% and 0.2027%, respectively, without vaccination. In this case, there are only about 2,510 human rabies infections at the end of tenth year. 4.4.2. Control program with dog culling only. The bounds for the dog culling rate, 0 ≤ u3 ≤ 0.22 and the weight parameters are taken from Section 4.2. All other model parameters are taken from Table 1. From the model simulations (Figure 9), we observe that without the dog culling program (u3 = 0), the Figure 9. The rabies prevalence among (a) dog population, (b) jackal population (c) number of human rabies cases and (d) optimal dog culling profile u3 in the control program with dog culling only. prevalence of rabies in dog population increases immediately after the outbreak begins. The prevalence then reaches as high as 0.2543% in the ninth year. Similarly, the rabies prevalence in jackal population peaks reaching to 0.2162% in eighth year and the number of human rabies reaches about 169,000. If the highest level of resource for culling is implemented (u3 = 0.22), the rabies prevalence in dog and jackal populations reaches negligible level (close to zero) at about 6 years. Also, the number of human rabies cases decreases rapidly and reaches about 6 at the end of tenth year. CONTROLLING RABIES IN NEPAL 369 In both cases, without dog culling (u3 = 0) and program with highest level of culling (u3 = 0.5), the objective functional remains higher at J = 3.5310 × 106 and J = 2.8437 × 106, respectively, asserting that neither of them is an optimal strategy. Our optimal control solution implies that with the dog culling strategy, u3(t), starting at the rate 0.184 and tapering down to zero, as shown in Figure 9d, can bring down the values of objective functional to 2.3017 × 106, which is 35% lower than without dog culling and 20% lower than the value in the highest culling rate. Under this optimal culling strategy, the peak rabies prevalence in dog and jackal populations are significantly lower than no culling strategy, while utilizing minimum resources. With this strategy, at the end of the study period of the tenth year, the rabies prevalence among dog and jackal populations remain 0.004% and 0.003%, respectively and the number of human rabies infections are about 2,608. 4.4.3. Control program with dog vaccination and dog culling combined. We now consider a combination of dog vaccination and dog culling strategies, and identify the optimal way of implementing them un- der resource limitation. As above, we take the range of u2 and u3 as 0 ≤ u2 ≤ 0.5, 0 ≤ u3 ≤ 0.22, respectively. As presented in Figure 10, under the combined program with the highest level of both Figure 10. The rabies prevalence among (a) dog population, (b) jackal population, (c) number of human rabies cases (d) optimal dog vaccination profile, u2, and (e) optimal dog culling profile, u3, in the control program with dog vaccination and dog culling combined. dog vaccination and dog culling (u2 = 0.5,u3 = 0.22), the prevalence of rabies in both dog and jackal populations as well as the number of human rabies cases approach to zero as early as in 3 years, while the long term prevalence remains 0.2543% and 0.2162%, respectively, in the absence of the program (u2 = 0,u3 = 0). In this combined approach, the optimal benefit can be achieved by implementing the dog vaccination and dog culling at the level of 0.11 and 0.09, respectively, at the beginning of the outbreak, and tapering down both to zero as shown in Figures 10(d,e). In this combined optimal strategy, the value of objective functional is J = 2.2532 × 106, which is about 37% lower than without control program (J = 3.5310×106) and about 67% lower than the highest level of both dog vaccination 370 B. PANTHA, H. JOSHI, AND N. VAIDYA and dog culling (J = 6.9613 × 106). With this strategy, the rabies prevalence stays approximately at the level of 0.00001% in both dog and jackal populations with only about 94 human rabies cases at the end of tenth year. As the more resources become available, we can utilize the weight constants of our optimal control model to represent the high resource scenario by assigning lower weight constants corresponding to the cost for controls (i.e., smaller Bi’s and Ci’s). For example, for 10-fold and 1000-fold higher resource availability, we take 10-fold and 1000-fold lower values of B2,C2,B3, and C3 than in the base case. For a higher resource availability (or lower B2,C2,B3,C3), the optimal control solution results in the optimal dog vaccination and dog culling strategies with a level higher than a case of limited resource (allowing the use of more resources) and tapering down (Figure 11). In this case, the value of the objective functional also comes out to be smaller (J = 2.2532×106 for the base case, J = 2.1062×106 for 10-fold higher resource, and J = 2.0724×106 for 1000-fold higher resource). For a higher resource availability, the number of human rabies cases as well as the prevalence of rabies in dog and jackal populations can be maintained at extremely low level with faster pace, particularly in the case of 1000-fold higher resource, compared to the case when resources are more limited (Figure 11). Figure 11. The rabies prevalence among (a) dog population, (b) jackal population, (c) number of human rabies cases (d) optimal dog vaccination profile, u2, and (e) optimal dog culling profile, u3, for 10-fold higher (red dashed curve) and 1000-fold higher(black dotted curve) resource availability (i.e., 10 times and 1000 times smaller values of B2,C2,B3,C3 than in the base case). 5. Conclusion Since the rabies is mostly problematic in developing counties of Asia and Africa, the control of rabies epidemics poses challenges due to limited resources available in these countries. The proper evalua- tion of control strategies and identifying optimal way of implementing such strategies are critical for CONTROLLING RABIES IN NEPAL 371 reducing rabies burden. In this study, we extended our basic rabies transmission dynamics model by adding the effects of commonly practiced control strategies, such as human post exposure vaccination, dog vaccination, dog culling, dog sterilization, and jackal bait vaccination. Furthermore, we formulated the optimal control model, which was then used to obtain the optimal time-varying strategy of imple- menting controls to mitigate rabies in Nepal under limited resources. Our model predicts that the dog vaccination and the dog culling are the most effective two control strategies to bring the basic reproduction number to low value, and also to reduce the human rabies cases and prevalence of rabies among dogs and jackals in Nepal. The optimal control formulation allowed us to identify time-dependent implementation of these two control strategies to achieve maximum benefit under limited resources. In general, applying higher level of controls at the beginning of outbreak and reducing during later part of the epidemic provide a maximum benefit, in both programs with single control strategy and two strategies combined. As revealed in our optimal control results, availability of more resources allows us to apply higher level of controls for longer period, resulting in lower level of rabies prevalence. We acknowledge that parameters used in our model are estimated from the literature or from the limited data set. Therefore, there may be some discrepancy between the model predictions and the actual prevalence of rabies in Nepal; more data sets from Nepal may help achieve better predictions of the model. The resource related parameters are particularly difficult to estimate, implying some uncertainty in quantitative optimal control results. However, the qualitative conclusion of results remain the same for a wider range of parameters, and therefore can be useful for implementing control strategies optimally. In summary, our model and optimal control theory provide framework to evaluate and implement control strategies under resource limited setting for mitigation of rabies epidemic burden in Nepal. References [1] E. Abo and E. Twizell, The Dynamics of a One-Dimensional Fox-Rabies Model, J. Appl. Math. & Computing, 22 (2006), 149–159. [2] J. Asamoah, F. Oduro, E. Bonyah, and B. Seldu, Modelling of Rabies Transmission Dynamics Using Optimal Control Analysis, Journal of Applied Mathematics, 2017(2017), Article ID 2451237, 23 pages. [3] R. Anderson, H. Jackson, R. May, and A. Smith, Population Dynamics of Fox Rabies in Europe, Nature, 289(1981), 765-771. [4] Center for Disease Control and Prevention: Rabies, Information for Specific Groups: Veterinarians, https://www.cdc.gov/rabies/index.html, accessed August, 2020. [5] H. Behncke, Optimal control of deterministic epidemics, Optimal Control Applications and Methods, 21(2000), 269-285. [6] Center for Disease Control and Prevention: Rabies, https://www.cdc.gov/rabies/index.html, accessed on December 2018. [7] T. Clayton, S. Duke-Sylvester, L. Gross, S. Lenhart, and L. Real, Optimal control of a rabies epidemic model with a birth pulse, Journal of Biological Dynamics, 4(2010), 43-58. [8] E. Demirci, A New Mathematical Approach for Rabies Endemy Applied Mathematical Sciences, 8 (2014), 59-67. [9] W. Ding, L. Gross, K. Langston, S. Lenhart, and L. Real, Rabies in raccoons: optimal control for a discrete time model on a spatial grid, Journal of Biological Dynamics, 1(2007), 379-393. [10] P. van den Driessche and J. Watmough, Reproduction number and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Bioscience, 180 (2002), 29-41. [11] T. Ega, L. Luboobi, and D. Kuznetsov, Modeling the Dynamics of Rabies Transmission with Vaccination and Stability Analysis, Applied and Computational Mathematics, 4(2015), 409-419. [12] Q. Hou, Z. Jin, and S. Ruan, Dynamics of rabies epidemics and the impact of control efforts in Guangdong Province, China, J Theor Biol., 300(2012), 39-47. 372 B. PANTHA, H. JOSHI, AND N. VAIDYA [13] H. R. Joshi, S. Lenhart, K. Albright and K. Gipson, Modeling the Effect of Information Campaigns On the HIV Epidemic In Uganda, Mathematical Biosciences and Engineering, 5(2008), 757-770. [14] The Kathmandu Post Nepal. https://kathmandupost.com/national/2019/02/04/nepal-has-pledged-to-end-rabies-by- 2030-but-dog-bite-cases-are-rising, accessed February 2019. [15] H. Laarabi, M. Rachik, O. E. Kahlaoui, and E. H. Labriji, Optimal Vaccination Strategies of an SIR Epidemic Model with a Saturated Treatment, Universal Journal of Applied Mathematics, 1(2013), 185-191. [16] S. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, 1st Edition, Chapman and Hall/CRC, 2007. [17] T. Leung and S. Davis, Rabies Vaccination Targets for Stray Dog Populations, Front. Vet. Sci.,4(2017): 52, 10 pages. [18] H. Liu, Spatial Spread of Rabies in Wildlife, Doctoral Dissertation. University of Arizona, 2013. [19] A. Mallela, S. Lenhart, and N. Vaidya, HIV-TB Co-infection Treatment: Modeling and Optimal Control Theory Perspectives, Journal of Computational and Applied Mathematics, 307(2016), 143-164. [20] R. Neilan, and S. Lenhart, Optimal vaccine distribution in a spatiotemporal epidemic model with an application to rabies and raccoons, J. Math. Anal. Appl., 378(2011), 603-619. [21] B. Pantha, S. Giri, H. Joshi, and N. Vaidya, Modeling Transmission Dynamics of Rabies in Nepal, Submitted. [22] V. Panjeti and L. Real, Mathematical Models for Rabies, Adv. Virus Res. 79(2011), 377-395. [23] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Inter Science Publishers, 1962. [24] C. Rhodes, R. Atkinson, R. Anderson, and D. Macdonald, Rabies in Zimbabwe: reservoir dogs and the implications for disease control, Phil. Trans. R. Soc. Lond. B., 353(1998), 999-1010. [25] B. Pantha, J. Day and S. Lenhart Optimal Control Applied in an Anthrax Epizootic Model., Journal of Biological System, 24(2016), 495-517. [26] S. Rimal, K. Ojha, W. Chaisowwong, Y. Shah, D. Pant, and A. Sirimalaisuwan, Detection of virus-neutralising antibodies and associated factors against rabies in the vaccinated household dogs of Kathmandu Valley, Nepal, PLoS One, 15 (2020): e0231967. [27] S. Ruan, Modeling the transmission dynamics and control of rabies in China, Mathematical Biosciences, 286 (2017), 65-93. [28] D. Slate, T. Algeo, K. Nelson, R. Chipman, D. Donovan, J. Blanton, M. Niezgoda, and C. Rupprecht, Oral rabies vaccination in North America: opportunities, complexities, and challenges, PLoS Negl. Trop. Dis. 3(2009): e549. [29] D. Lukes, Differential Equations: Classical to Controlled, Mathematics in Science and Engineering, Academic Press, 1982. [30] K. Tohma, M. Saito, C. Demetria, D. Manalo, B. Quiambao, H. Kamigaki, Molecular and mathematical modeling analyses of inter-island transmission of rabies into a previously rabies-free island in the Philippines, Infection, Genetics and Evolution, 38(2016), 22-28. [31] S. Yadav, Animal Rabies in Nepal and Raccoon Rabies in Albany County, New York, Report for Masters of Public Health, Kansas State University, 2012. [32] J. Zhang, Z. Jin, G. Sun, T. Zhou, and S. Ruan, Analysis of Rabies in China: Transmission Dynamics and Control, PLoS One, 6(2011): e20891. [33] World Health Organization/Nepal: http://www.searo.who.int/nepal/en/, Accessed in October 2018. Corresponding author, Department of Science and Mathematics, Abraham Baldwin Agricultural College, Tifton, GA, 31793, USA E-mail address: bpantha@abac.edu Department of Mathematics, Xavier University, Cincinnati, OH, USA E-mail address: joshi@xavier.edu Department of Mathematics and Statistics, Computational Science Research Center, & Viral Information Institute, San Diego State University, San Diego, CA, USA E-mail address: nvaidya@sdsu.edu