# Modelling of Rabies Transmission Dynamics Using Optimal Control Analysis

Authors: Joshua Kiddy K. Asamoah; Francis T. Oduro; Ebenezer Bonyah; Baba Seidu
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

We examine an optimal way of eradicating rabies transmission from dogs into the human population, using preexposure prophylaxis (vaccination) and postexposure prophylaxis (treatment) due to public education. We obtain the disease-free equilibrium, the endemic equilibrium, the stability, and the sensitivity analysis of the optimal control model. Using the Latin hypercube sampling (LHS), the forward-backward sweep scheme and the fourth-order Range-Kutta numerical method predict that the global alliance for rabies control’s aim of working to eliminate deaths from canine rabies by 2030 is attainable through mass vaccination of susceptible dogs and continuous use of pre- and postexposure prophylaxis in humans.


---

## Body

## 1. Introduction

Rabies is an infection that mostly affects the brain of an infected animal or individual, caused by viruses belonging to the genusLyssavirus of the family Rhabdoviridae and order Mononegavirales [1, 2]. This disease has become a global threat and it is also estimated that rabies occurs in more than 150 countries and territories [2]. Raccoons, skunks, bats, and foxes are the main animals that transmit the virus in the United States [2]. In Asia, Africa, and Latin America, it is known that dogs are the main source of transmission of the rabies virus into the human population [2]. When the rabies virus enters the human body or that of an animal, the infection (virus) moves rapidly along the neural pathways to the central nervous system; from there the virus continues to spread to other organs and causes injury by interrupting various nerves [2]. The symptoms of rabies are quite similar to those of encephalitis (see [3]). Due to movement of dogs in homes or the surroundings, the risk of not being infected by a rabid dog can never be guaranteed. Rabies is a major health problem in many populations dense with dogs, especially in areas where there are less or no preventive measures (vaccination and treatment) for dogs and humans. Treatment after exposure to the rabies virus is known as postexposure prophylaxis (PEP) and vaccination before exposure to the infection is known as preexposure prophylaxis.The study of optimal control analysis in maximizing or minimizing a said target was introduced by Pontryagin and his collaborators around 1950. They developed the key idea of introducing the adjoint function to a differential equation, by forming an objective functional [4], and since then there has been a considerable study of infectious disease using optimal control analysis (see [4–12]).Research published by Aubert [13], on the advancement of the expense of wildlife rabies in France, incorporated various variables. They follow immunization of domestic animals, the reinforcement of epidemiological reconnaissance system and the bolster given to indicative research laboratories, the costs connected with outbreaks of rabies, the clinical perception of those mammals which had bitten humans, the preventive immunization, and postexposure treatment of people. A significant percentage (72%) of the cost was the preventive immunization of local animals. In France, as in other European nations in which the red fox (Vulpes) is the species most affected, two primary procedures for controlling rabies were assessed in [13] at the repository level to be specific: fox termination and the oral immunization of foxes. The consolidated costs and advantages of both systems were looked at and included either the expenses of fox separation or the cost of oral immunization. The total yearly costs of both techniques stayed practically identical until the fourth year, after which the oral immunization methodology turned out to be more cost effective. This estimate was made in 1988 and readjusted in 1993 and affirmed by ex-postinvestigation five years later. Accordingly, it was presumed that fox termination brought about a transient diminishment in the event of the infection while oral immunization turned out to be equipped for wiping out rabies even in circumstances in which fox population was growing. Anderson and May [14] formulated a mathematical model based on each time step dynamic which was calculated independently in every cell. Later, Bohrer et al. [15] published a paper on the viability of different rabies spatial immunization designs in a simulated host population.The research presented by Bohrer [15] stated that, in desert environments, where host population size varies over time, nonuniform spreading of oral rabies vaccination may, under certain circumstances, be more effective than the commonly used uniform spread. The viability of a nonarbitrary spread of the immunization depends, to some extent, on the dispersal behavior of the carriers. The outcomes likewise exhibit that, in a warm domain in a few high-density regions encompassed by populations with densities below the critical threshold for the spread of the disease, the rabies infection can persist.Levin et al. [16] also presented a model for the immune responses to rabies virus in bats. Coyne et al. [17] proposed an SEIR model, which was also used in a study predicting the local dynamics of rabies among raccoons in the United States. Childs et al. [18] also researched rabies epidemics in raccoons with a seasonal birth pulse, using optimal control of an SEIRS model which describes the population dynamics. Hampson et al. [19] also noted that rabies epidemic cycles have a period of 3–6 years in dog populations in Africa, so they built a susceptible, exposed, infectious, and vaccinate model with an intervention response variable, which showed significant synchrony.Carroll et al. [20] also used compartmental models to describe rabies epidemiology in dog populations and explored three control methods: vaccination, vaccination pulse fertility control, and culling. An ordinary differential equation model was used to characterize the transmission dynamics of rabies between humans and dogs by [21, 22]. The work by Zinsstag et al. [23] further extended the existing models on rabies transmission between dogs to include dog-to-human transmission and concluded that human postexposure prophylaxis (PEP) with a dog vaccination campaign was the more cost effective in controlling the disease in the long run. Furthermore, Ding et al. [24] formulated an epidemic model for rabies in raccoons with discrete time and spatial features. Their goal was to analyze the strategies for optimal distribution of vaccine baits to minimize the spread of the disease and the cost of carrying out the control. Smith and Cheeseman [25] show that culling could be more effective than vaccination, given the same efficacy of control, but Tchuenche and Bauch suggest that culling could be counterproductive, for some parameter values (see [26]).The work in [27, 28] also presented a mathematical model of rabies transmission in dogs and from the dog population to the human population in China. Their study did not consider the use optimal control analysis to the study of the rabies virus in dogs and from the dog population to the human population. Furthermore, the insightful work of Wiraningsih et al. [29] studied the stability analysis of a rabies model with vaccination effect and culling in dogs, where they introduced postexposure prophylaxis to a rabies transmission model, but the paper did not consider the noneffectiveness of the pre- and postprophylaxis on the susceptible humans and exposed humans and that of the dog population and the use of optimal control analysis. Therefore, motivated by the research predictions of the global alliance of rabies control [30] and the work mention above, we seek to adjust the model presented in [27–29], by formulating an optimal control model, so as to ascertain an optimal way of controlling rabies transmission in dogs and from the dog population to the human population taking into account the noneffectiveness (failure) of vaccination and treatment.The paper is petition as follows. Section2 contains the model formulation, mathematical assumptions, the mathematical flowchart, and the model equations. Section 3 contains the model analysis, invariant region, equilibrium points, basic reproduction number R0, and the stability analysis of the equilibria. In Section 4 we present the parameter values leading to numerical values of the basic reproduction number R0, the herd immunity threshold and sensitivity analysis using Latin hypercube sampling (LHS), and some numerical plots. Section 5 contains the objective functional and the optimality system of the model. Finally, Sections 6 and 7 contain discussion and conclusion, respectively.


## 2. Model Formulation

We present two subpopulation transmission models of rabies virus in dogs and that of the human population (see Figure1), based on the work presented in [27–29]. The dog population has a total of four compartments. The compartments represent the susceptible dogs, SD(t), exposed dogs, ED, infected dogs, ID(t), and partially immune dogs, RD(t). Thus, the total dog population is ND(t)=SD(t)+ED(t)+ID(t)+RD(t). The human population also has four compartments representing susceptible humans, SH(t), exposed humans, EH(t), infected humans, IH(t), and partially immune humans, RH(t). Thus, the total human population is NH(t)=SH(t)+EH(t)+IH(t)+RH(t). It is assumed that there is no human to human transmission of the rabies virus in the human submodel (see [29]). In the dog submodel, it is assumed that there is a direct transmission of the rabies virus from one dog to the other and from the infected dog compartment to the susceptible human population. It is further assumed that the susceptible dog population, SD(t), is increased by recruitment at a rate AD and BH is the birth or immigration rate into the susceptible human population, SH(t). It is assumed that the transmission and contact rate of the rabid dog into the dog compartment is βDD. Suppose that νD represents the control strategy due to public education and vaccination in the dogs compartment; then the transmission dynamics become (1-νD)βDDSDID, where (1-νD) is the noneffectiveness (failure) of the vaccine. It is also assumed that the contact rate of infectious dogs to the human population is βDH. Similarly, administrating vaccination to the susceptible humans the progression rate of the susceptible humans to the exposed stage becomes (1-νH)βDHSHID, where νH is the preexposure prophylaxis (vaccination), (1-νH) represents the failure of the preexposure prophylaxis in the human compartment. Furthermore, administrating postexposure prophylaxis (treatment) to affected humans at the rate ρH decreases the progression rate of the rabies virus, at the exposed class to the infectious class as (1-ρH)δHγHEH, where (1-ρH) is the failure rate of the postexposure prophylaxis and δHγH represents the rate at which exposed humans progress to the infected compartment [27]. The rate of losing immunity in both compartments is represented by αD and αH, respectively.Figure 1
Optimal control model of rabies transmission dynamics.The exposed humans without clinical rabies that move back to the susceptible population are denoted by the rateδHεH. The natural death rate of dogs is mD, and mH denotes the mortality rate of humans (natural death rate), μD represents the death rate associated with rabies infection in dogs, and μH represents the disease induce death in humans. The rate at which exposed dogs die due to culling is CD, and δεD represents the rate at which exposed dogs without clinical rabies move back to the susceptible dog compartment. Subsequently, using the idea presented in [29], we assumed that the exposed dogs are treated or quarantined by their owners at the rate ρD; this implies that (1-ρD)δγDED is the progression rate of the exposed dogs to the infectious compartment, where (1-ρD) is the failure of the treatment or quarantined strategy, and δγDED denotes those exposed dogs that develop clinical rabies [27]. Figure 1 shows the mathematical dynamics of the rabies virus in both compartments.From Figure1 transmission flowchart and assumptions give the disease pathways as(1)dSDdt=AD-1-νDβDDSDID-mD+νDSD+δεDED+αDRD,dEDdt=1-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED,dIDdt=1-ρDδγDED-mD+μDID,dRDdt=νDSD+ρDED-mD+αDRD,dSHdt=BH-1-νHβDHSHID-mH+νHSH+δHεHEH+αHRH,dEHdt=1-νHβDHSHID-1-ρHδHγH+mH+ρH+δHεHEH,dIHdt=1-ρHδHγHEH-mH+μHIH,dRHdt=νHSH+ρHEH-mH+αHRH,withSD0>0,ED0≥0,ID0≥0,RD0≥0,SH0>0,EH0≥0,IH0>0,RH0>0.


## 3. Model Analysis

Model system (1) will be studied in a biological feasible region as outlined below. Model system (1) is basically divided into two regions; thus Ω=ΩD×ΩH.Lemma 1.
The solution setSD,ED,ID,RD,SH,EH,IH,RH∈R+8 of model system (1) is contained in the feasible region Ω.Proof.
SupposeSD,ED,ID,RD,SH,EH,IH,RH∈R+8 for all t>0. We want to show that the region Ω is positively invariant, so that it becomes sufficient to look at the dynamics of model system (1), given that(2)NDt=SDt+EDt+IDt+RDt,(3)NHt=SHt+EHt+IHt+RHt,where ND(t) is the total population of dogs at any time (t) and NH(t) is total population of humans at any time (t).
Equation (2) gives(4)dNDdt=AD-SD+ED+ID+RDmD-μDID-CDED,which yields(5)dNDdt=AD-mDND-μDID-CDED.
Similarly (3) gives(6)dNHdt=BH-mHNH-μHIH.
Now, assuming that there are no disease induced death rate and culling effect in the dogs’ compartment, it implies that (5) and (6) become(7)dNDdt=AD-mDND,dNHdt=BD-mHNH.Suppose dND/dt≤0, dNH/dt≤0, ND≤AD/mD, and NH≤BH/mH, and then imposing the theorem proposed in [32] on differential inequality results in 0≤ND≤AD/mD and 0≤NH≤BH/mH. Therefore (7) becomes(8)dNDdt≤AD-mDND,(9)dNHdt≤BD-mHNH.
Solve (8) and (9) using the integrating factor (IF) method. Thus dy/dt+p(t)y=Q, IF=e∫‍p(t)dt. After some algebraic manipulation the feasible solution of the dogs’ population in model system (1) is in the region(10)ΩD=SD,ED,ID,RD∈R+4,ND≤ADmD.
Similarly the human population follows suit, and from (9) this implies that the feasible solution of the human population of model system (1) is in the region(11)ΩH=SH,EH,IH,RH∈R+4,NH≤BHmH. Therefore, the feasible solutions are contained in Ω. Thus Ω=ΩD×ΩH. From the standard comparison theorem used on differential inequality in [33], it implies that(12)NDt≤ND0e-mDt+ADmD1-e-mDt,NHt≤NH0e-mHt+BHmH1-e-mHt.
Hence, the total dog population sizeND(t)→AD/mD as t→∞. Similarly, the total human population size NH(t)→BH/mH as t→∞. This means that the infected state variables ED,ID,EH,IH of the two populations tend to zero as time goes to infinity. Therefore, the region Ω is pulling (attracting) all the solutions in R+8. This gives the feasible solution set of model system (1) as(13)SDEDIDRDSHEHIHRH∈R+8∣SD>0ED≥0ID≥0ID≥0RD≥0SH>0EH≥0IH≥0RH≥0ND≤ADmDNH≤BHmH.Hence, (1) is mathematically well posed and epidemiologically meaningful.


### 3.1. Disease-Free EquilibriumE0

Suppose there is no infection of rabies in both compartments; thenED=0,ID=0,EH=0,IH=0. Incorporating this into (1) leads to (14)AD-mD+νDSD+αDRD=0,νDSD-mD+αDRD=0,BH-mH+νHSH+αHRH=0,νHSH-mH+αHRH=0.After some algebraic manipulation of (14), the disease-free equilibrium point becomes E0=SD0,ED0,ID0,RD0,SH0,EH0,IH0,RH0 with (15)E0=ADmD+αDmDmD+αD+νD,0,0,ADνDmDmD+αD+νD,BHmH+αHmHmH+αH+νH,0,0,BHνHmHmH+αH+νH.


### 3.2. Basic Reproduction NumberR0

Here, the basic reproduction number(R0) measures the average number of new infections produced by one infected dog in a completely susceptible (dog and human) population (see also [34]). Now taking ED,ID,EH, and IH as our infected compartments gives(16)f1=1-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED,f2=1-ρDδγDED-mD+μDID,f3=1-νHβDHSHID-1-ρHδHγH+mH+ρH+δHεHEH,f4=1-ρHδHγHEH-mH+μHIH,where f1=dED/dt,f2=dID/dt,f3=dEH/dt, and f4=dIH/dt.Now, using the next generation matrix operatorG=FV-1 and the Jacobian matrix(17)J=∂f1∂ED∂f1∂ID∂f1∂EH∂f1∂IH∂f2∂ED∂f2∂ID∂f2∂EH∂f2∂IH∂f3∂ED∂f3∂ID∂f3∂EH∂f3∂IH∂f4∂ED∂f4∂ID∂f4∂EH∂f4∂IH, as described in [34], results in(18)J=-1-ρDδγD+mD+ρD+δεD+CD1-νDβDDSD001-ρDδγD-mD+μD0001-νHβDHSH-1-ρHδHγH+mH+ρH+δHεH0001-ρHδHγH-mH+μH.Using the fact thatJ=F-V gives F and V evaluated at E0 as(19)FE0=01-νDβDDADmD+αDmDmD+νD+αD00000001-νHβDHmH+αHBHmHmH+νH+αH000000,VE0=1-ρDδγD+mD+ρD+δεD+CD000-1-ρDδγDmD+μD00001-ρHδHγH+mH+ρH+δHεH000-1-ρHδHγHmH+μH,where the element in matrix F constitutes the new infection terms, while that of matrix V constitutes the new transfer of infection terms from one compartment to another. Now, splitting matrix V into four 2×2 submatrices and finding its corresponding inverses result in G=FV-1, given by(20)G=1-ρD1-νDδγDβDDADmD+αD1-ρDδγD+mD+ρD+δεD+CDmD+μDmDmD+νD+αD1-νDβDDADmD+αDmDmD+νD+αD0000001-ρDδγD1-νHβDHBHmH+αH1-ρDδγD+mD+ρD+δεD+CDmD+μDmHmH+νH+αH1-νHβDHmH+αHBHmD+νDmHmH+νH+αH000000Letting(21)a=1-ρD1-νDδγDβDDADmD+αD1-ρDδγD+mD+ρD+δεD+CDmD+μDmDmD+νD+αD,b=1-νDβDDADmD+αDmDmD+νD+αD,c=1-ρD1-νHδγDβDHmH+αH1-ρDδγD+mD+ρD+δεD+CDmD+μDmHmH+νH+αH,d=1-νHβDHmH+αHBHmD+νDmHmH+νH+αHimplies(22)G=ab000000cd000000.Finding the matrix determinant of (22) and denoting it by D give the expression D=|G-Iλ|, where I is the identity matrix of a 4×4 matrix; thus(23)D=a-λb000-λ00cd-λ0000-λ=0.This gives a characteristic equation of the formλ3(a-λ)=0; solving the characteristic polynomial results in the following eigenvalues: λi=0,0,0,a. The basic reproduction number R0 is the spectral radius (largest eigenvalue) ρFV-1, also defined as the dominant eigenvalue of FV-1.Therefore,(24)R0=1-ρD1-νDδγDβDDADmD+αD1-ρDδγD+mD+ρD+δεD+CDmD+μDmDmD+νD+αD.Remark 2.
R
0 contains the secondary infection produced by the infectious compartment of dogs (in the presence of preexposure prophylaxis (vaccination), postexposure prophylaxis (treatment/quarantine), and culling of exposed dogs). When R0<1, the infection gradually leaves the dog compartment, but when R0>1, the rabies virus remains in the dog compartments for a longer time, thereby increasing the rate at which the susceptible dogs and humans get infected by a rabid dog.


### 3.3. Endemic EquilibriumE1

The endemic equilibrium is given as(25)SD∗=ADmD+αDmDmD+νD+αDR0,ED∗=mD+μD1-ρDδγDID∗,ID∗=1-ρDδγD+mD+ρD+δγDmD+μDmDmD+νD+αDR0-1mD+αD1-νDβDD1-ρDδγD+mD+CD+mD1-νDβDDρD,RD∗=ADνD1-νDβDD1-ρDδγDmD+αD+1-νDβDDρDmD+μDID∗mDmD+νD+αDR01-νDβDD1-ρDδγDmD+αD,SH∗=BHmH+αH+δHεH+αHρHEH∗1-νHmH+αHβDHID∗+mHmH+αH+νH,EH∗=1-νHBHmH+αHβDHID∗mH+αH1-νHβDHID∗1-ρHδHγH+mH+ρH+mH+νH1-ρHδHγH+mH+ρH+δHεH-1-νHβDHID∗αHρH,IH∗=1-ρHδHγHmH+μHEH∗,RH∗=BHνHmH+νH+νHδHεH+νHαHρH+ρH1-νHmH+αHβDHID∗+ρHmHmH+αH+νHEH∗1-νHmH+αH2βDHID∗+mH+αHmHmH+αH+νH.Note that ifR0=1, it results in the disease-free equilibrium; if R0>1, then there exists a unique endemic equilibrium; if R0<1, then there exist two endemic equilibriums.


### 3.4. Stability Analysis ofE0

Linearizing (1) at E0 and subtracting eigenvalue λ along the main diagonal yield(26)JE0=b1-λb7a1αD00000a2-λa3000000b10b2-λ00000νDρD0b3-λ000000a40b4-λb50αH00a500a6-λ0000000b9b6-λ00000νHρH0b8-λ,where(27)a1=-1-νDβDDmD+αDADmDmD+νD+αD,a2=-1-ρDδγD+mD+ρD+δεD+CD,a3=1-νDβDDmD+αDADmDmD+νD+αD,a4=-1-νHβDHmH+αHmHmH+νH+αH,a5=1-νHβDHmH+αHmHmH+νH+αH,a6=-1-ρHδHγH+mH+ρH+δHεH,b1=-mD+νD,b2=-mD+μD,b3=-mD+αD,b4=-νH+mH,b5=δHεH,b6=-mH+μH,b7=δεD,b8=-mH+αH,b9=1-ρHδHγH,b10=1-ρDδγD.Simplifying matrixJ(E0) gives(28)b6-λa6-λb4-λb8-λλ4+a11λ3+a12λ2+a13λ+a14=0,where(29)a11=-b2-a2-b1-b3,a12=νDαD+a2b3+a2b1+b2b3+b2b1+b3b1+b2a2-1-ρDδγDa3,a13=-a2νDαD-b2νDαD+a31-ρDδγDb2+a31-ρDδγDb1-a2b2b3-b2a2b1-a2b3b1-b2b3b1,a14=b1b2b3a2+1-ρDδγDa3b3b1+1-ρDδγDa3νD+νDαDb2a2.From (28) the four characteristic factors that are negative are (30)λ1=b6,λ2=a6,λ3=b4,λ4=b8,where a6=-((1-ρH)δHγH+mH+ρH+δH+δHεH), b6=-(mH+μH), b4=-(νH+mH), and b8=-(mH+αH). The other four characteristic factors can be obtained using the Routh-Hurwitz criterion. Routh-Hurwitz stability criterion is a test to ascertain the nature of the eigenvalues. If the roots of the polynomial are all positive, then the polynomial has a negative real part [35, 36]. The remaining four characteristic eigenvalues are obtained as follows:(31)λ4+a11λ3+a12λ2+a13λ+a14=0.Hence, simplifying the coefficient of the above characteristic polynomial in (31) yields (32)a11=1-ρDδγD+mD+ρD+δεD+CD+mD+μD+mD+αD+mD+νD,a12=νDαD+1-ρDδγD+mD+ρD+δεD+CDmD+αD+mD+νD+mD+μDmD+αD+mD+νD+mD+αDmD+νD+1-ρDδγD+mD+ρD+CDmD+μD1-R0,a13=1-ρDδγD+mD+ρD+δεD+CDνDαD+mD+μDνDαD+mD+μDmD+αDmD+νD+1-ρDδγD+mD+ρD+δεD+CDmD+αDmD+μD1-R0mD+μDmD+αD+1-ρDδγD+mD+ρD+δεD+CDmD+αDmD+νD1-R0m+α,a14=νDαDmD+νD1-ρDδγD+mD+ρD+δεD+CD+1-ρDδγD1-νDβDDmD+αDADνDmDmD+mD+νD+αD+mD+νDmD+μDmD+αD1-ρDδDγD+mD+ρD+δεD+CD1-R0.Therefore, from the Routh-Hurwitz criterion of order four, it implies that the conditions,a11>0, a12>0, a13>0, a14>0, and a11a12a13>a132+a112a14, are satisfied if R0<1. Hence, the disease-free equilibrium E0 is locally asymptotically stable when R0<1 (see [37]).


#### 3.4.1. Global Stability ofE0

Theorem 3.
The disease-free equilibriumE0 of model (1) is globally asymptotically stable if R0≤1 and unstable if R0>1.Proof.
LetV be a Lyapunov function with positive constants K1, K2, K3, and K4 such that(33)V=SD-SD0-SD0ln⁡SDSD0+K1ED+K2ID+RD-RD0-RD0ln⁡RDRD0+SH-SH0-SH0ln⁡SHSH0+K3EH+K4IH+RH-RH0-RH0ln⁡RHRH0.
Taken the derivative of the Lyapunov function with respect to time gives(34)dVdt=1-SD0SDdSDdt+K1dEDdt+K2dIDdt+1-RD0RDdRDdt+1-SH0SHdSHdt+K3dEHdt+K4dIHdt+1-RH0RHdRHdt.
Plugging (1) into (34) results in(35)dVdt=1-SD0SDAD-1-νDβDDSDID-mD+νDSD+δεDED+αRD+K11-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED+K21-ρDδγDED-mD+μDID+1-RD0RDνDSD+ρDED-mD+αDRD+1-SH0SHBH-1-νHβDHSHID-mH+νHSH+δHεHEH+αHRH+K31-νHβDHSHID-1-ρHδHγH+mH+ρH+δHεHEH+K41-ρHδHγHEH-mH+μHIH+1-RH0RHνHSH+ρHEH-mH+αHRH.
Now, after forming the Lyapunov functionV on the space of the eight state variables, thus (SD,ED,ID,RD,SH,EH,IH,RH), and introducing the idea from [37], it is clear that if ED(t), ID(t), EH(t), and IH(t) at the disease-free equilibrium are globally stable (thus, ED=0, ID=0, EH=0, and IH=0), then SD(t)→AD(mD+αD)/mD(mD+αD+νD), RD(t)→ADνD/mD(mD+αD+νD),  SH(t)→BH(mH+αH)/mH(mH+αH+nuH), and RH(t)→BHνH/mH(mH+αH+H+νH) as t→∞.
Therefore, it can be assumed that(36)SD≤SD0=ADmD+αDmDmD+αD+νD,RD≤RD0=ADνDmDmD+αD+νD,SH≤SH0=BHmH+αHmHmH+αH+nuH,RH≤RH0=BHνHmHmH+αH+H+νH,(see [38]) and replacing it into (35) yields(37)dVdt≤K11-νDβDDADmD+αDmDmD+αD+νDID-1-ρDδγD+mD+ρD+δεD+CDED+K21-ρDδγDED-mD+μDID+K31-νHβDHBHmH+αHmHmH+αH+νHId-1-ρHδHγH+mH+ρH+δHεHEH+K41-ρHδHγHEH-mH+μHIH,
This implies that(38)dVdt≤K11-νDβDDADmD+αDmDmD+αD+νD-K2mD+μD+K31-νHβDHBHmH+αHmHmH+αH+νHID+K21-ρDδγD-K11-ρDδγD+mD+ρD+δεD+CDED+K41-ρHδHγH-K31-ρHδHγH+mH+δHεHEH-K4mH+μH.
Equating the coefficient ofID,  ED,  IH, and EH in (38) to zero gives(39)K4=K3=0,K2=1-ρDδγD+mD+ρD+δεD+CD,K1=1-ρDδγD, and we obtain (40)dVdt≤1-ρDδγD+mD+ρD+δεD+CDmD+μDR0-1ID,≤0,ifR0≤1.
AdditionallydV/dt=0 if and only if ID=0. Therefore, for ED=ID=EH=IH=0 it shows that SDt→ADmD+αD/mDmD+αD+νD, RDt→ADνD/mDmD+αD+νD, SH(t)→BH(mH+αH)/mH(mH+αH+νH), and RH(t)→BHνH/mH(mH+αH+νH) as t→∞. Hence, the largest compact invariant set in (SD,ED,ID,RD,SH,EH,IH,RH)∈Ω:dV/dt≤0 is the singleton set E0. Therefore, from La Salle’s invariance principle, we conclude that E0 is globally asymptotically stable in Ω if R0≤1 (see also [38, 39]).


### 3.5. Global Stability of Endemic EquilibriumE1

Theorem 4.
The endemic equilibriumE1 of model (1) is globally asymptotically stable whenever R0>1.Proof.
SupposeR0>1; then the existence of the endemic equilibrium point is assured. Using the common quadratic Lyapunov function(41)Vx1,x2,…,xn=∑i=1nci2xi-xi∗2, as illustrated in [40], we consider a Lyapunov function with the following candidate:(42)VSD,ED,ID,RD,SH,EH,IH,RH=12SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗2+12SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗2.
Now, differentiating (42) along the solution curve of (1) gives(43)dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗dSD+ED+ID+RDdt+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗dSH+EH+IH+RHdt.
From (1) it implies that d(SD+ED+ID+RD)/dt=AD-mD(SD+ED+ID+RD)-CDED-μDID and d(SH+EH+IH+RH)/dt=B-m(SH+EH+IH+RH)-μHIH, which when plugged into (43) gives(44)dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗AD-mDSD+ED+ID+RD-CDED-μDID+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗BH-mSH+EH+IH+RH-μHIH.
Now assuming(45)AD=mDSD∗+ED∗+ID∗+RD∗+CDED∗+μDID∗,BH=mHSH∗+EH∗+IH∗+RH∗+μHIH∗ and substituting it into (44), we have(46)dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗mDSD∗+ED∗+ID∗+RD∗+CDED∗+μDID∗-mDSD+ED+ID+RD-CDED-μDID+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗mHSH∗+EH∗+IH∗+RH∗+μHIH∗-mSH+EH+IH+RH-μHIH,dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗-mDSD-SD∗-mDED-ED∗-mDID-ID∗-mDRD-RD∗-CDED-ED∗-μDID-ID∗+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗-mHSH-SH∗-mHEH-EH∗-mHIH-IH∗-mHRH-RH∗-μHIH-IH∗.
This also implies that(47)dVdt=-mDSD-SD∗2-CD+mDED-ED∗2-mD+μDID-ID∗2-mDRD-RD∗2-2mD+CDSD-SD∗ED-ED∗-2mD+μDSD-SD∗ID-ID∗-2mD+μD+CDED-ED∗ID-ID∗-2mDRD-RD∗ID-ID∗-2mD+μD+CDRD-RD∗ID-ID∗-mHSH-SH∗2-mHEH-EH∗2-mH-μHIH-IH∗2-mHRH-RH∗2-2mHSH-SH∗EH-EH∗-2mH-μHSH-SH∗IH-IH∗-2mH+μHEH-EH∗IH-IH∗-mHIH-IH∗RH-RH∗+SH-SH∗RH-RH∗.
This shows thatdV/dt is negative and dV/dt=0, if and only if SD=SD∗, ED=ED∗, ID=ID∗, RD=RD∗, SH=SH∗, EH=EH∗, IH=IH∗, RH=RH∗. Additionally every solution of (1) with the initial conditions approaches E1 as t→∞ (see [38, 39]); therefore, the largest compact invariant set in (SD,ED,ID,RD,SH,EH,IH,RH)∈Ω:dV/dt≤0 is the singleton set E1. Therefore, from Lasalle’s invariant principle [41], it implies that the endemic equilibrium E1 is globally asymptotically stable in Ω whenever R0>1.


## 3.1. Disease-Free EquilibriumE0

Suppose there is no infection of rabies in both compartments; thenED=0,ID=0,EH=0,IH=0. Incorporating this into (1) leads to (14)AD-mD+νDSD+αDRD=0,νDSD-mD+αDRD=0,BH-mH+νHSH+αHRH=0,νHSH-mH+αHRH=0.After some algebraic manipulation of (14), the disease-free equilibrium point becomes E0=SD0,ED0,ID0,RD0,SH0,EH0,IH0,RH0 with (15)E0=ADmD+αDmDmD+αD+νD,0,0,ADνDmDmD+αD+νD,BHmH+αHmHmH+αH+νH,0,0,BHνHmHmH+αH+νH.


## 3.2. Basic Reproduction NumberR0

Here, the basic reproduction number(R0) measures the average number of new infections produced by one infected dog in a completely susceptible (dog and human) population (see also [34]). Now taking ED,ID,EH, and IH as our infected compartments gives(16)f1=1-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED,f2=1-ρDδγDED-mD+μDID,f3=1-νHβDHSHID-1-ρHδHγH+mH+ρH+δHεHEH,f4=1-ρHδHγHEH-mH+μHIH,where f1=dED/dt,f2=dID/dt,f3=dEH/dt, and f4=dIH/dt.Now, using the next generation matrix operatorG=FV-1 and the Jacobian matrix(17)J=∂f1∂ED∂f1∂ID∂f1∂EH∂f1∂IH∂f2∂ED∂f2∂ID∂f2∂EH∂f2∂IH∂f3∂ED∂f3∂ID∂f3∂EH∂f3∂IH∂f4∂ED∂f4∂ID∂f4∂EH∂f4∂IH, as described in [34], results in(18)J=-1-ρDδγD+mD+ρD+δεD+CD1-νDβDDSD001-ρDδγD-mD+μD0001-νHβDHSH-1-ρHδHγH+mH+ρH+δHεH0001-ρHδHγH-mH+μH.Using the fact thatJ=F-V gives F and V evaluated at E0 as(19)FE0=01-νDβDDADmD+αDmDmD+νD+αD00000001-νHβDHmH+αHBHmHmH+νH+αH000000,VE0=1-ρDδγD+mD+ρD+δεD+CD000-1-ρDδγDmD+μD00001-ρHδHγH+mH+ρH+δHεH000-1-ρHδHγHmH+μH,where the element in matrix F constitutes the new infection terms, while that of matrix V constitutes the new transfer of infection terms from one compartment to another. Now, splitting matrix V into four 2×2 submatrices and finding its corresponding inverses result in G=FV-1, given by(20)G=1-ρD1-νDδγDβDDADmD+αD1-ρDδγD+mD+ρD+δεD+CDmD+μDmDmD+νD+αD1-νDβDDADmD+αDmDmD+νD+αD0000001-ρDδγD1-νHβDHBHmH+αH1-ρDδγD+mD+ρD+δεD+CDmD+μDmHmH+νH+αH1-νHβDHmH+αHBHmD+νDmHmH+νH+αH000000Letting(21)a=1-ρD1-νDδγDβDDADmD+αD1-ρDδγD+mD+ρD+δεD+CDmD+μDmDmD+νD+αD,b=1-νDβDDADmD+αDmDmD+νD+αD,c=1-ρD1-νHδγDβDHmH+αH1-ρDδγD+mD+ρD+δεD+CDmD+μDmHmH+νH+αH,d=1-νHβDHmH+αHBHmD+νDmHmH+νH+αHimplies(22)G=ab000000cd000000.Finding the matrix determinant of (22) and denoting it by D give the expression D=|G-Iλ|, where I is the identity matrix of a 4×4 matrix; thus(23)D=a-λb000-λ00cd-λ0000-λ=0.This gives a characteristic equation of the formλ3(a-λ)=0; solving the characteristic polynomial results in the following eigenvalues: λi=0,0,0,a. The basic reproduction number R0 is the spectral radius (largest eigenvalue) ρFV-1, also defined as the dominant eigenvalue of FV-1.Therefore,(24)R0=1-ρD1-νDδγDβDDADmD+αD1-ρDδγD+mD+ρD+δεD+CDmD+μDmDmD+νD+αD.Remark 2.
R
0 contains the secondary infection produced by the infectious compartment of dogs (in the presence of preexposure prophylaxis (vaccination), postexposure prophylaxis (treatment/quarantine), and culling of exposed dogs). When R0<1, the infection gradually leaves the dog compartment, but when R0>1, the rabies virus remains in the dog compartments for a longer time, thereby increasing the rate at which the susceptible dogs and humans get infected by a rabid dog.


## 3.3. Endemic EquilibriumE1

The endemic equilibrium is given as(25)SD∗=ADmD+αDmDmD+νD+αDR0,ED∗=mD+μD1-ρDδγDID∗,ID∗=1-ρDδγD+mD+ρD+δγDmD+μDmDmD+νD+αDR0-1mD+αD1-νDβDD1-ρDδγD+mD+CD+mD1-νDβDDρD,RD∗=ADνD1-νDβDD1-ρDδγDmD+αD+1-νDβDDρDmD+μDID∗mDmD+νD+αDR01-νDβDD1-ρDδγDmD+αD,SH∗=BHmH+αH+δHεH+αHρHEH∗1-νHmH+αHβDHID∗+mHmH+αH+νH,EH∗=1-νHBHmH+αHβDHID∗mH+αH1-νHβDHID∗1-ρHδHγH+mH+ρH+mH+νH1-ρHδHγH+mH+ρH+δHεH-1-νHβDHID∗αHρH,IH∗=1-ρHδHγHmH+μHEH∗,RH∗=BHνHmH+νH+νHδHεH+νHαHρH+ρH1-νHmH+αHβDHID∗+ρHmHmH+αH+νHEH∗1-νHmH+αH2βDHID∗+mH+αHmHmH+αH+νH.Note that ifR0=1, it results in the disease-free equilibrium; if R0>1, then there exists a unique endemic equilibrium; if R0<1, then there exist two endemic equilibriums.


## 3.4. Stability Analysis ofE0

Linearizing (1) at E0 and subtracting eigenvalue λ along the main diagonal yield(26)JE0=b1-λb7a1αD00000a2-λa3000000b10b2-λ00000νDρD0b3-λ000000a40b4-λb50αH00a500a6-λ0000000b9b6-λ00000νHρH0b8-λ,where(27)a1=-1-νDβDDmD+αDADmDmD+νD+αD,a2=-1-ρDδγD+mD+ρD+δεD+CD,a3=1-νDβDDmD+αDADmDmD+νD+αD,a4=-1-νHβDHmH+αHmHmH+νH+αH,a5=1-νHβDHmH+αHmHmH+νH+αH,a6=-1-ρHδHγH+mH+ρH+δHεH,b1=-mD+νD,b2=-mD+μD,b3=-mD+αD,b4=-νH+mH,b5=δHεH,b6=-mH+μH,b7=δεD,b8=-mH+αH,b9=1-ρHδHγH,b10=1-ρDδγD.Simplifying matrixJ(E0) gives(28)b6-λa6-λb4-λb8-λλ4+a11λ3+a12λ2+a13λ+a14=0,where(29)a11=-b2-a2-b1-b3,a12=νDαD+a2b3+a2b1+b2b3+b2b1+b3b1+b2a2-1-ρDδγDa3,a13=-a2νDαD-b2νDαD+a31-ρDδγDb2+a31-ρDδγDb1-a2b2b3-b2a2b1-a2b3b1-b2b3b1,a14=b1b2b3a2+1-ρDδγDa3b3b1+1-ρDδγDa3νD+νDαDb2a2.From (28) the four characteristic factors that are negative are (30)λ1=b6,λ2=a6,λ3=b4,λ4=b8,where a6=-((1-ρH)δHγH+mH+ρH+δH+δHεH), b6=-(mH+μH), b4=-(νH+mH), and b8=-(mH+αH). The other four characteristic factors can be obtained using the Routh-Hurwitz criterion. Routh-Hurwitz stability criterion is a test to ascertain the nature of the eigenvalues. If the roots of the polynomial are all positive, then the polynomial has a negative real part [35, 36]. The remaining four characteristic eigenvalues are obtained as follows:(31)λ4+a11λ3+a12λ2+a13λ+a14=0.Hence, simplifying the coefficient of the above characteristic polynomial in (31) yields (32)a11=1-ρDδγD+mD+ρD+δεD+CD+mD+μD+mD+αD+mD+νD,a12=νDαD+1-ρDδγD+mD+ρD+δεD+CDmD+αD+mD+νD+mD+μDmD+αD+mD+νD+mD+αDmD+νD+1-ρDδγD+mD+ρD+CDmD+μD1-R0,a13=1-ρDδγD+mD+ρD+δεD+CDνDαD+mD+μDνDαD+mD+μDmD+αDmD+νD+1-ρDδγD+mD+ρD+δεD+CDmD+αDmD+μD1-R0mD+μDmD+αD+1-ρDδγD+mD+ρD+δεD+CDmD+αDmD+νD1-R0m+α,a14=νDαDmD+νD1-ρDδγD+mD+ρD+δεD+CD+1-ρDδγD1-νDβDDmD+αDADνDmDmD+mD+νD+αD+mD+νDmD+μDmD+αD1-ρDδDγD+mD+ρD+δεD+CD1-R0.Therefore, from the Routh-Hurwitz criterion of order four, it implies that the conditions,a11>0, a12>0, a13>0, a14>0, and a11a12a13>a132+a112a14, are satisfied if R0<1. Hence, the disease-free equilibrium E0 is locally asymptotically stable when R0<1 (see [37]).


### 3.4.1. Global Stability ofE0

Theorem 3.
The disease-free equilibriumE0 of model (1) is globally asymptotically stable if R0≤1 and unstable if R0>1.Proof.
LetV be a Lyapunov function with positive constants K1, K2, K3, and K4 such that(33)V=SD-SD0-SD0ln⁡SDSD0+K1ED+K2ID+RD-RD0-RD0ln⁡RDRD0+SH-SH0-SH0ln⁡SHSH0+K3EH+K4IH+RH-RH0-RH0ln⁡RHRH0.
Taken the derivative of the Lyapunov function with respect to time gives(34)dVdt=1-SD0SDdSDdt+K1dEDdt+K2dIDdt+1-RD0RDdRDdt+1-SH0SHdSHdt+K3dEHdt+K4dIHdt+1-RH0RHdRHdt.
Plugging (1) into (34) results in(35)dVdt=1-SD0SDAD-1-νDβDDSDID-mD+νDSD+δεDED+αRD+K11-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED+K21-ρDδγDED-mD+μDID+1-RD0RDνDSD+ρDED-mD+αDRD+1-SH0SHBH-1-νHβDHSHID-mH+νHSH+δHεHEH+αHRH+K31-νHβDHSHID-1-ρHδHγH+mH+ρH+δHεHEH+K41-ρHδHγHEH-mH+μHIH+1-RH0RHνHSH+ρHEH-mH+αHRH.
Now, after forming the Lyapunov functionV on the space of the eight state variables, thus (SD,ED,ID,RD,SH,EH,IH,RH), and introducing the idea from [37], it is clear that if ED(t), ID(t), EH(t), and IH(t) at the disease-free equilibrium are globally stable (thus, ED=0, ID=0, EH=0, and IH=0), then SD(t)→AD(mD+αD)/mD(mD+αD+νD), RD(t)→ADνD/mD(mD+αD+νD),  SH(t)→BH(mH+αH)/mH(mH+αH+nuH), and RH(t)→BHνH/mH(mH+αH+H+νH) as t→∞.
Therefore, it can be assumed that(36)SD≤SD0=ADmD+αDmDmD+αD+νD,RD≤RD0=ADνDmDmD+αD+νD,SH≤SH0=BHmH+αHmHmH+αH+nuH,RH≤RH0=BHνHmHmH+αH+H+νH,(see [38]) and replacing it into (35) yields(37)dVdt≤K11-νDβDDADmD+αDmDmD+αD+νDID-1-ρDδγD+mD+ρD+δεD+CDED+K21-ρDδγDED-mD+μDID+K31-νHβDHBHmH+αHmHmH+αH+νHId-1-ρHδHγH+mH+ρH+δHεHEH+K41-ρHδHγHEH-mH+μHIH,
This implies that(38)dVdt≤K11-νDβDDADmD+αDmDmD+αD+νD-K2mD+μD+K31-νHβDHBHmH+αHmHmH+αH+νHID+K21-ρDδγD-K11-ρDδγD+mD+ρD+δεD+CDED+K41-ρHδHγH-K31-ρHδHγH+mH+δHεHEH-K4mH+μH.
Equating the coefficient ofID,  ED,  IH, and EH in (38) to zero gives(39)K4=K3=0,K2=1-ρDδγD+mD+ρD+δεD+CD,K1=1-ρDδγD, and we obtain (40)dVdt≤1-ρDδγD+mD+ρD+δεD+CDmD+μDR0-1ID,≤0,ifR0≤1.
AdditionallydV/dt=0 if and only if ID=0. Therefore, for ED=ID=EH=IH=0 it shows that SDt→ADmD+αD/mDmD+αD+νD, RDt→ADνD/mDmD+αD+νD, SH(t)→BH(mH+αH)/mH(mH+αH+νH), and RH(t)→BHνH/mH(mH+αH+νH) as t→∞. Hence, the largest compact invariant set in (SD,ED,ID,RD,SH,EH,IH,RH)∈Ω:dV/dt≤0 is the singleton set E0. Therefore, from La Salle’s invariance principle, we conclude that E0 is globally asymptotically stable in Ω if R0≤1 (see also [38, 39]).


## 3.4.1. Global Stability ofE0

Theorem 3.
The disease-free equilibriumE0 of model (1) is globally asymptotically stable if R0≤1 and unstable if R0>1.Proof.
LetV be a Lyapunov function with positive constants K1, K2, K3, and K4 such that(33)V=SD-SD0-SD0ln⁡SDSD0+K1ED+K2ID+RD-RD0-RD0ln⁡RDRD0+SH-SH0-SH0ln⁡SHSH0+K3EH+K4IH+RH-RH0-RH0ln⁡RHRH0.
Taken the derivative of the Lyapunov function with respect to time gives(34)dVdt=1-SD0SDdSDdt+K1dEDdt+K2dIDdt+1-RD0RDdRDdt+1-SH0SHdSHdt+K3dEHdt+K4dIHdt+1-RH0RHdRHdt.
Plugging (1) into (34) results in(35)dVdt=1-SD0SDAD-1-νDβDDSDID-mD+νDSD+δεDED+αRD+K11-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED+K21-ρDδγDED-mD+μDID+1-RD0RDνDSD+ρDED-mD+αDRD+1-SH0SHBH-1-νHβDHSHID-mH+νHSH+δHεHEH+αHRH+K31-νHβDHSHID-1-ρHδHγH+mH+ρH+δHεHEH+K41-ρHδHγHEH-mH+μHIH+1-RH0RHνHSH+ρHEH-mH+αHRH.
Now, after forming the Lyapunov functionV on the space of the eight state variables, thus (SD,ED,ID,RD,SH,EH,IH,RH), and introducing the idea from [37], it is clear that if ED(t), ID(t), EH(t), and IH(t) at the disease-free equilibrium are globally stable (thus, ED=0, ID=0, EH=0, and IH=0), then SD(t)→AD(mD+αD)/mD(mD+αD+νD), RD(t)→ADνD/mD(mD+αD+νD),  SH(t)→BH(mH+αH)/mH(mH+αH+nuH), and RH(t)→BHνH/mH(mH+αH+H+νH) as t→∞.
Therefore, it can be assumed that(36)SD≤SD0=ADmD+αDmDmD+αD+νD,RD≤RD0=ADνDmDmD+αD+νD,SH≤SH0=BHmH+αHmHmH+αH+nuH,RH≤RH0=BHνHmHmH+αH+H+νH,(see [38]) and replacing it into (35) yields(37)dVdt≤K11-νDβDDADmD+αDmDmD+αD+νDID-1-ρDδγD+mD+ρD+δεD+CDED+K21-ρDδγDED-mD+μDID+K31-νHβDHBHmH+αHmHmH+αH+νHId-1-ρHδHγH+mH+ρH+δHεHEH+K41-ρHδHγHEH-mH+μHIH,
This implies that(38)dVdt≤K11-νDβDDADmD+αDmDmD+αD+νD-K2mD+μD+K31-νHβDHBHmH+αHmHmH+αH+νHID+K21-ρDδγD-K11-ρDδγD+mD+ρD+δεD+CDED+K41-ρHδHγH-K31-ρHδHγH+mH+δHεHEH-K4mH+μH.
Equating the coefficient ofID,  ED,  IH, and EH in (38) to zero gives(39)K4=K3=0,K2=1-ρDδγD+mD+ρD+δεD+CD,K1=1-ρDδγD, and we obtain (40)dVdt≤1-ρDδγD+mD+ρD+δεD+CDmD+μDR0-1ID,≤0,ifR0≤1.
AdditionallydV/dt=0 if and only if ID=0. Therefore, for ED=ID=EH=IH=0 it shows that SDt→ADmD+αD/mDmD+αD+νD, RDt→ADνD/mDmD+αD+νD, SH(t)→BH(mH+αH)/mH(mH+αH+νH), and RH(t)→BHνH/mH(mH+αH+νH) as t→∞. Hence, the largest compact invariant set in (SD,ED,ID,RD,SH,EH,IH,RH)∈Ω:dV/dt≤0 is the singleton set E0. Therefore, from La Salle’s invariance principle, we conclude that E0 is globally asymptotically stable in Ω if R0≤1 (see also [38, 39]).


## 3.5. Global Stability of Endemic EquilibriumE1

Theorem 4.
The endemic equilibriumE1 of model (1) is globally asymptotically stable whenever R0>1.Proof.
SupposeR0>1; then the existence of the endemic equilibrium point is assured. Using the common quadratic Lyapunov function(41)Vx1,x2,…,xn=∑i=1nci2xi-xi∗2, as illustrated in [40], we consider a Lyapunov function with the following candidate:(42)VSD,ED,ID,RD,SH,EH,IH,RH=12SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗2+12SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗2.
Now, differentiating (42) along the solution curve of (1) gives(43)dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗dSD+ED+ID+RDdt+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗dSH+EH+IH+RHdt.
From (1) it implies that d(SD+ED+ID+RD)/dt=AD-mD(SD+ED+ID+RD)-CDED-μDID and d(SH+EH+IH+RH)/dt=B-m(SH+EH+IH+RH)-μHIH, which when plugged into (43) gives(44)dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗AD-mDSD+ED+ID+RD-CDED-μDID+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗BH-mSH+EH+IH+RH-μHIH.
Now assuming(45)AD=mDSD∗+ED∗+ID∗+RD∗+CDED∗+μDID∗,BH=mHSH∗+EH∗+IH∗+RH∗+μHIH∗ and substituting it into (44), we have(46)dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗mDSD∗+ED∗+ID∗+RD∗+CDED∗+μDID∗-mDSD+ED+ID+RD-CDED-μDID+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗mHSH∗+EH∗+IH∗+RH∗+μHIH∗-mSH+EH+IH+RH-μHIH,dVdt=SD-SD∗+ED-ED∗+ID-ID∗+RD-RD∗-mDSD-SD∗-mDED-ED∗-mDID-ID∗-mDRD-RD∗-CDED-ED∗-μDID-ID∗+SH-SH∗+EH-EH∗+IH-IH∗+RH-RH∗-mHSH-SH∗-mHEH-EH∗-mHIH-IH∗-mHRH-RH∗-μHIH-IH∗.
This also implies that(47)dVdt=-mDSD-SD∗2-CD+mDED-ED∗2-mD+μDID-ID∗2-mDRD-RD∗2-2mD+CDSD-SD∗ED-ED∗-2mD+μDSD-SD∗ID-ID∗-2mD+μD+CDED-ED∗ID-ID∗-2mDRD-RD∗ID-ID∗-2mD+μD+CDRD-RD∗ID-ID∗-mHSH-SH∗2-mHEH-EH∗2-mH-μHIH-IH∗2-mHRH-RH∗2-2mHSH-SH∗EH-EH∗-2mH-μHSH-SH∗IH-IH∗-2mH+μHEH-EH∗IH-IH∗-mHIH-IH∗RH-RH∗+SH-SH∗RH-RH∗.
This shows thatdV/dt is negative and dV/dt=0, if and only if SD=SD∗, ED=ED∗, ID=ID∗, RD=RD∗, SH=SH∗, EH=EH∗, IH=IH∗, RH=RH∗. Additionally every solution of (1) with the initial conditions approaches E1 as t→∞ (see [38, 39]); therefore, the largest compact invariant set in (SD,ED,ID,RD,SH,EH,IH,RH)∈Ω:dV/dt≤0 is the singleton set E1. Therefore, from Lasalle’s invariant principle [41], it implies that the endemic equilibrium E1 is globally asymptotically stable in Ω whenever R0>1.


## 4. Numerical Analysis

Considering the parameter values in Table1, we will ascertain the numerical importance of our analysis.Table 1
Parameter values.
Parameter
Description
Standard value
Source
A
D
Recruitment rate of dogs
3
×
1
0
6
y
-
1
[27]
α
D
Loss of immunity in dogs
1
y
-
1
[27]
C
D
Death rate of dogs due to culling
0.3
y
-
1
Assumed
m
D
Natural death rate of dogs
0.056
y
-
1
[27]
μ
D
Disease induced mortality in dogs
1
y
-
1
[27]
ν
D
Preexposure prophylaxis for dogs
0.25
y
-
1
Assumed
ρ
D
Postexposure prophylaxis for dogs
0.2
y
-
1
[27]
β
D
D
Transmission rate in dogs
1.58
×
1
0
-
7
y
-
1
[27]
γ
D
Latency period in dogs
2.37
/
6
y
-
1
[27]
δ
ε
D
Rate of no clinical rabies
0.4
y
-
1
[27]
B
H
Birth rate (humans)
0.0314
y
-
1
[31]
β
D
H
Transmission rate (dog-humans)
2.29
×
1
0
-
12
y
-
1
[27]
α
H
Loss of immunity (humans)
1
y
-
1
[27]
m
H
Natural death rate (humans)
0.0074
y
-
1
[31]
μ
H
Disease induced mortality (humans)
1
y
-
1
[27]
ν
H
Preexposure prophylaxis for humans
0.54
y
-
1
Assumed
ρ
H
Postexposure prophylaxis for humans
0.1
y
-
1
[27]
γ
H
Latency rate (humans)
1
/
6
y
-
1
[27]
γ
H
ε
H
Rate of no clinical rabies (humans)
2.4
y
-
1
[27]


### 4.1. Different Scenarios of the Basic Reproduction NumberR0

We shall denoteR0 without pre- and postexposure prophylaxis (treatment) as R0∗ and R0 without preexposure prophylaxis and culling as R0∗∗ and the R0 without postexposure prophylaxis (treatment) and culling as R0∗∗∗. Therefore, using the parameter values in Table 1, R0∗, R0∗∗, and R0∗∗∗ are given as follows: (48)R0∗=βDDADδγDδγD+mD+δεD+CD×mD+μDmD,R0∗=3.027,R0∗∗=1-ρDδγDβDDAD1-ρDδγD+mD+ρD+δεDmD+μDmD,R0∗∗=2.181,R0∗∗∗=1-νDδγDβDDADmD+αDδγD+mD+δεDmD+μDmDmD+νD+αD,R0∗∗∗=1.914.Therefore, from the above calculations it indicates that the best way in reducing or minimizing the rabies virus in the dogs compartment is to use more of preexposure prophylaxis (vaccination).


#### 4.1.1. Herd Immunity ThresholdH1

Therefore, from the above numerical values, we are motivated to know the number of humans or dogs that should be vaccinated whenR0∗=3.027. (49)H1≔1-1R0∗=0.66.This shows that ifR0∗=3.027, then 66% of individuals and dogs should receive vaccination.


### 4.2. Sensitivity Analysis

To determine parameters that contribute most to the rabies transmission, we used two sensitivity analysis approach: the normalised forward sensitivity index as presented in [37] and the Latin hypercube sampling as described in [42]. To determine the dependence of parameters in R0, using a sampling size, n=1000, the partial rank correction coefficients (PRCC) value of the ten parameters in R0 are shown in Figure 2(a). The longer the bar in Figure 2(a) suggests that the statistical influence of those parameters to changes in R0 is high. Also, using the normalised forward sensitivity index gives the following values and the nature of their signs in Table 2, based on the parameter value given in Table 1. The plus sign or minus sign signifies that the influence is positive or negative, respectively [42],(50)ΓR0βDD=∂R0∂βDβDDR0=1,ΓR0AD=∂R0∂ADADR0=1,ΓR0μD=∂R0∂μDμDR0=-μDmD+μD=-0.95,ΓR0δεD=∂R0∂δεDδεDR0=δεD1-ρDδγD-δεD-CD-mD-ρD=-1.61,ΓR0CD=∂R0∂CDCDR0=CD1-ρDδγD-δεD-CD-mD-ρD=-0.45,ΓR0αD=∂R0∂αDαDR0=0.28,ΓR0mD=∂R0∂mDmDR0=-1.64,ΓR0δγD=∂R0∂δγDδγDR0=1.33,ΓR0ρD=∂R0∂ρDρDR0=-0.5,ΓR0νD=∂R0∂νDνDR0=-0.52.Table 2
Sensitivity signs ofR0 to the parameters in (24).
Parameter
Description
Sensitivity sign
βDD
Transmission rate of dogs
+ve
AD
Recruitment rate of dogs
+ve
μD
Disease induce death rate of dogs
−ve
δεD
Rate of no clinical rabies
−ve
CD
Culling of exposed dogs
−ve
αD
Loss of immunity in dogs
+ve
mD
Natural death rate of dogs
−ve
δγD
Rate at which exposed dogs become infective (infective rate)
+ve
ρD
Postexposure prophylaxis (treatment/quarantined)
−ve
νD
Preexposure prophylaxis (vaccination)
−veFigure 2
The graphical representation of some parameters inR0 and the effect of varying some initial state values on the model.
(a)
LHS plot for the parameters inR0
(b)
Effect of varying recruitment rate on the infected humans
(c)
Effect of an increase inR0 on the infected humans
(d)
Effect of varying the initial infected dog population size on the infected humansTherefore, from Table2 it shows that an addition or a reduction in the values of βDD,  αD,  δγD, and AD will have an increase or a decrease in the spread of the rabies virus. For example, ΓR0βD=1 indicates that increasing or reducing the transmission rate by 5% may increase or reduce the number of secondary infection by 5%. The negative sign in Table 2 will have a reduction in the basic reproduction number, R0, when the values of those parameters are increased, and a reduction in the values of ρD, νD, μD, mD, and δεD will lead to an increase in the number of secondary infections.The Latin hypercube sampling (LHS) in Figure2(a) shows that μD, CD, αD, and δγD have a minimal influence on the rate at which the rabies virus is spread. The Latin hypercube sampling (LHS) plots for the ten parameters in R0 show that culling of exposed dogs does not actually minimize the spread of rabies as compared to vaccination of susceptible dogs. Figure 2(a) also shows that the most influential parameter in spreading the infection is βDD followed by AD. Figure 2(c) shows that an increase in the basic reproduction number will contribute to a high level of secondary infection in the human population. Similarly, Figure 2(a) shows that vaccination of dogs νD is the most effective way of controlling the rabies virus in the dog population as compared to the treatment/quarantine of exposed dogs, ρD. Figure 3(a) gives the contour nature of νD and ρD, which shows a more saturated effect on the basic reproduction number. Figure 3(b) shows that βDD and  αD have a positive relation with the basic reproduction number R0. Therefore, an increase in βDD and αD will have a direct increase in the spread of the rabies virus. Figure 2(b) indicates that with a high number of recruitment of dogs into the susceptible dog’s compartment will have a corresponding high increase in the number of infected humans. Figure 2(d) demonstrates that a high number of infected dogs in the compartment will lead to an increase in the number of infected humans. Figure 3(c) shows that a high increase in the number of disease induce death rate and natural death rate will have a negative reflection on R0; biologically, we would not recommend this approach in minimizing the spread of the disease, since an increase in both μD and mD may result in a high rate of the disease in the human population, even though μD and mD naturally reduce the number of susceptible and infected dogs in the population. Finally, Figure 3(d) shows the 3D plot of Figure 3(a).Figure 3
The graphical representation of some parameters inR0.
(a)
The contour plot ofνD and ρD to R0
(b)
The 3D plot ofR0 to αD and βDD
(c)
The 3D plot ofR0 to mD and μD
(d)
The 3D plot ofR0 to ρD and νD


## 4.1. Different Scenarios of the Basic Reproduction NumberR0

We shall denoteR0 without pre- and postexposure prophylaxis (treatment) as R0∗ and R0 without preexposure prophylaxis and culling as R0∗∗ and the R0 without postexposure prophylaxis (treatment) and culling as R0∗∗∗. Therefore, using the parameter values in Table 1, R0∗, R0∗∗, and R0∗∗∗ are given as follows: (48)R0∗=βDDADδγDδγD+mD+δεD+CD×mD+μDmD,R0∗=3.027,R0∗∗=1-ρDδγDβDDAD1-ρDδγD+mD+ρD+δεDmD+μDmD,R0∗∗=2.181,R0∗∗∗=1-νDδγDβDDADmD+αDδγD+mD+δεDmD+μDmDmD+νD+αD,R0∗∗∗=1.914.Therefore, from the above calculations it indicates that the best way in reducing or minimizing the rabies virus in the dogs compartment is to use more of preexposure prophylaxis (vaccination).


### 4.1.1. Herd Immunity ThresholdH1

Therefore, from the above numerical values, we are motivated to know the number of humans or dogs that should be vaccinated whenR0∗=3.027. (49)H1≔1-1R0∗=0.66.This shows that ifR0∗=3.027, then 66% of individuals and dogs should receive vaccination.


## 4.1.1. Herd Immunity ThresholdH1

Therefore, from the above numerical values, we are motivated to know the number of humans or dogs that should be vaccinated whenR0∗=3.027. (49)H1≔1-1R0∗=0.66.This shows that ifR0∗=3.027, then 66% of individuals and dogs should receive vaccination.


## 4.2. Sensitivity Analysis

To determine parameters that contribute most to the rabies transmission, we used two sensitivity analysis approach: the normalised forward sensitivity index as presented in [37] and the Latin hypercube sampling as described in [42]. To determine the dependence of parameters in R0, using a sampling size, n=1000, the partial rank correction coefficients (PRCC) value of the ten parameters in R0 are shown in Figure 2(a). The longer the bar in Figure 2(a) suggests that the statistical influence of those parameters to changes in R0 is high. Also, using the normalised forward sensitivity index gives the following values and the nature of their signs in Table 2, based on the parameter value given in Table 1. The plus sign or minus sign signifies that the influence is positive or negative, respectively [42],(50)ΓR0βDD=∂R0∂βDβDDR0=1,ΓR0AD=∂R0∂ADADR0=1,ΓR0μD=∂R0∂μDμDR0=-μDmD+μD=-0.95,ΓR0δεD=∂R0∂δεDδεDR0=δεD1-ρDδγD-δεD-CD-mD-ρD=-1.61,ΓR0CD=∂R0∂CDCDR0=CD1-ρDδγD-δεD-CD-mD-ρD=-0.45,ΓR0αD=∂R0∂αDαDR0=0.28,ΓR0mD=∂R0∂mDmDR0=-1.64,ΓR0δγD=∂R0∂δγDδγDR0=1.33,ΓR0ρD=∂R0∂ρDρDR0=-0.5,ΓR0νD=∂R0∂νDνDR0=-0.52.Table 2
Sensitivity signs ofR0 to the parameters in (24).
Parameter
Description
Sensitivity sign
βDD
Transmission rate of dogs
+ve
AD
Recruitment rate of dogs
+ve
μD
Disease induce death rate of dogs
−ve
δεD
Rate of no clinical rabies
−ve
CD
Culling of exposed dogs
−ve
αD
Loss of immunity in dogs
+ve
mD
Natural death rate of dogs
−ve
δγD
Rate at which exposed dogs become infective (infective rate)
+ve
ρD
Postexposure prophylaxis (treatment/quarantined)
−ve
νD
Preexposure prophylaxis (vaccination)
−veFigure 2
The graphical representation of some parameters inR0 and the effect of varying some initial state values on the model.
(a)
LHS plot for the parameters inR0
(b)
Effect of varying recruitment rate on the infected humans
(c)
Effect of an increase inR0 on the infected humans
(d)
Effect of varying the initial infected dog population size on the infected humansTherefore, from Table2 it shows that an addition or a reduction in the values of βDD,  αD,  δγD, and AD will have an increase or a decrease in the spread of the rabies virus. For example, ΓR0βD=1 indicates that increasing or reducing the transmission rate by 5% may increase or reduce the number of secondary infection by 5%. The negative sign in Table 2 will have a reduction in the basic reproduction number, R0, when the values of those parameters are increased, and a reduction in the values of ρD, νD, μD, mD, and δεD will lead to an increase in the number of secondary infections.The Latin hypercube sampling (LHS) in Figure2(a) shows that μD, CD, αD, and δγD have a minimal influence on the rate at which the rabies virus is spread. The Latin hypercube sampling (LHS) plots for the ten parameters in R0 show that culling of exposed dogs does not actually minimize the spread of rabies as compared to vaccination of susceptible dogs. Figure 2(a) also shows that the most influential parameter in spreading the infection is βDD followed by AD. Figure 2(c) shows that an increase in the basic reproduction number will contribute to a high level of secondary infection in the human population. Similarly, Figure 2(a) shows that vaccination of dogs νD is the most effective way of controlling the rabies virus in the dog population as compared to the treatment/quarantine of exposed dogs, ρD. Figure 3(a) gives the contour nature of νD and ρD, which shows a more saturated effect on the basic reproduction number. Figure 3(b) shows that βDD and  αD have a positive relation with the basic reproduction number R0. Therefore, an increase in βDD and αD will have a direct increase in the spread of the rabies virus. Figure 2(b) indicates that with a high number of recruitment of dogs into the susceptible dog’s compartment will have a corresponding high increase in the number of infected humans. Figure 2(d) demonstrates that a high number of infected dogs in the compartment will lead to an increase in the number of infected humans. Figure 3(c) shows that a high increase in the number of disease induce death rate and natural death rate will have a negative reflection on R0; biologically, we would not recommend this approach in minimizing the spread of the disease, since an increase in both μD and mD may result in a high rate of the disease in the human population, even though μD and mD naturally reduce the number of susceptible and infected dogs in the population. Finally, Figure 3(d) shows the 3D plot of Figure 3(a).Figure 3
The graphical representation of some parameters inR0.
(a)
The contour plot ofνD and ρD to R0
(b)
The 3D plot ofR0 to αD and βDD
(c)
The 3D plot ofR0 to mD and μD
(d)
The 3D plot ofR0 to ρD and νD


## 5. Objective Functional

Given thaty(t)∈Y∈Rn is a state variable of model system (1) and u(t)∈U∈Rn are the control variables at any time (t) with t(0)≤t≤t(f), then an optimal control problem consists of finding a piecewise continuous control u(t) and its corresponding state y(t). This optimizes the cost functional J[y(t),u(t)] using Pontryagin’s maximum principle [43]. Therefore we set the following likelihood control strategies:(1)
u1=νD is the control effort aimed at increasing the immunity of susceptible dogs (preexposed prophylaxis).(2)
u2=ρD is the control effort aimed at treating the exposed dogs (postexposed prophylaxis).(3)
u3=νH is the control effort aimed at increasing the immunity of susceptible humans (preexposure prophylaxis).(4)
u4=ρH is the control effort aimed at treating the exposed humans (postexposed prophylaxis).Our goal is to seek optimal controls such asνD∗, ρD∗, νH∗, and ρH∗ that minimize the objective functional:(51)J=min∫t0tfA1ED+A2EH+A3ID+A4IH+B12νD2+B22ρD2+B32νH2+B42ρH2dt.Therefore, (51) is subject to (52)dSddt=AD-1-νDβDDSDID-mD+νDSD+δεDED+αDRD,dEddt=1-νDβDDSDID-1-ρDδγD+mD+ρD+δεD+CDED,dIDdt=1-ρDδγDED-mD+μDID,dRDdt=νDSD+ρDED-mD+αDRD,dSHdt=BH-1-νHβDHSHId-mH+νHSH+δHεHEH+αHRH,dEHdt=1-νHβDHSHId-1-ρHδHγH+mH+ρH+δHεHEH,dIHdt=1-ρHδHγHEH-mH+μHIH,dRHdt=νHSH+ρHEH-mH+αHRH,SD>0,ED≥0,ID≥0,RD≥0,SH>0,EH≥0,IH≥0,RH≥0.From (51) the quantities A1 and A2 denote the weight constants of the exposed classes and A3 and A4 are the weight of the infectious classes, respectively. B1,B2,B3,B4 are the weight constants for the dog and human controls. B1νD2,B2ρD2,B3νH2,B4ρH2 describe the cost associated with rabies vaccination and treatment. The square of the control variables shows the severity of the side effects of the vaccination and treatment. Employing Pontryagin’s maximum principle, we form the Hamiltonian equation with state variables SD=SD∗,ED=ED∗,ID=ID∗,RD∗ and SH=SH∗,EH=EH∗,IH=IH∗,RH∗ as (53)H=A1ED∗+A2EH∗+A3ID∗+A4IH∗+B12νD2+B22ρD2+B32νH2+B42ρH2+λ1AD-1-νDβDDSD∗ID∗-mD+νDSD∗+δεED∗+αDRD∗+λ21-νDβDDSD∗ID∗-1-ρDδγD+mD+ρD+δεD+CDED∗+λ31-ρDδγDED∗-mD+μDID∗+λ4νDSD∗+ρDED∗-mD+αDRD∗+λ5BH-1-νHβDHSH∗ID∗-mH+νHSH∗+δHεHEH∗+αHRH∗+λ61-νHβDHSH∗ID∗-1-ρHδHγH+mH+ρH+δHεHEH∗+λ71-ρHδHγHEH∗-mH+μHIH∗+λ8νHSH∗+ρHEH∗-mH+αHRH∗.Considering the existence of adjoint functionsλi, i=1,2,…,8, satisfying (54)dλ1dt=-∂H∂SD∗=λ11-νDβDDID∗+mD+νD-λ21-νDβDDID∗-λ4νD,dλ2dt=-∂H∂ED∗=λ21-ρDδγD+mD+ρD+δεD+CD-λ1δεD-λ31-ρDδγD-λ4ρD-A1,dλ3dt=-∂H∂ID∗=λ3mD+μD+λ11-νDβDDSD∗+λ51-νHβDHSH∗-λ21-νDβDSD∗-λ61-νHβDHSH∗-A3,dλ4dt=-∂H∂RD∗=λ4mD+αD-λ1αD,dλ5dt=-∂H∂SH∗=λ51-νHβdHID∗+mH+νH-λ61-νHβDHID∗-λ8νH,dλ6dt=-∂H∂EH∗=λ61-ρHδHγH+mH+ρH+δHεH-λ5δHεH-λ71-ρHδHγH-λ8ρH-A2,dλ7dt=-∂H∂IH∗=λ7mH+μH-A4,dλ8dt=-∂H∂RH∗=λ8mH+αH-λ5αH, with transversality condition λi(tf)=0 for i=1,…,8 for the control set ui, hence we have(55)∂H∂ui=0,wherei=1,2,3,4,∂H∂νDνD=νD∗≔B1νD∗-λ1SD∗+λ4SD∗+λ1βDDSD∗ID∗-λ2βDSD∗ID∗=0,νD∗=λ1SD∗-λ4SD∗+λ2-λ1βDDID∗SD∗B1,∂H∂ρDρD=ρD∗≔B2ρD∗-λ1ED∗+λ4ED∗+λ2EDδγDED∗-λ3δγDED∗=0,ρD∗=λ2ED∗-λ4ED∗+λ3-λ2δγDED∗B2,∂H∂νHνH=νH∗≔B3νH∗-λ5SH+λ8SH+λ5βDHSHID-λ6βDHSH∗ID∗=0,νH∗=λ5SH∗-λ8SH∗+λ6-λ5βDHSH∗ID∗B3,∂H∂ρHρH=ρH∗≔B4ρ∗-λ6EH∗+λ8EH∗+λ6δHγHEH∗-λ7δHγHEH∗=0,ρH∗=λ6EH∗-λ8EH∗+λ7-λ6δHγHEH∗B4.Now, using an appropriate variation argument and taking the bounds into account, the optimal control strategies are given as(56)νD∗=min⁡max⁡0,λ1-λ4SD∗+λ2-λ1βDDID∗SD∗B1,νDmax,(57)ρD∗=min⁡max⁡0,λ2-λ4ED∗+λ3-λ2δγDED∗B2,ρDmax,(58)νH∗=min⁡max⁡0,λ5-λ8SH∗+λ6-λ5βDHSH∗ID∗B3,νHmax,(59)ρH∗=min⁡max⁡0,λ6-λ8EH∗+λ7-λ6δHγHEH∗B4,ρHmax.Optimality System. Substituting the representation of the optimal vaccination and treatment control with corresponding adjoint function, we have the optimality system as(60)dSDdt=AD-1-minmax0,λ1-λ4SD∗+λ2-λ1βDDID∗SD∗B1,νDmaxβDDSDID-mDSD-minmax0,λ1-λ4SD∗+λ2-λ1βDDID∗SD∗B1,νDmaxSD+δεDED+αDRD,dEDdt=1-minmax0,λ1-λ4SD∗+λ2-λ1βDDID∗SD∗B1,νDmaxβDDSDID-1-minmax0,λ2-λ4ED∗+λ3-λ2δγDED∗B2,ρmaxδγD+mD+δεD+CDED-minmax0,λ2-λ4ED∗+λ3-λ2δγDED∗B2,ρDmaxED,dIDdt=δγDED-mD+μDID,dRDdt=minmax0,λ1-λ4SD∗+λ2-λ1βDDID∗SD∗B1,νDmaxSD-mD+αDRD+minmax0,λ2-λ4ED∗+λ3-λ2δγDED∗B2,ρDmaxED,dSHdt=BH-1-minmax0,λ5-λ8SH∗+λ6-λ5βDHSH∗ID∗B3,νHmaxβDHSHID-mHSH-minmax0,λ5-λ8SH∗+λ6-λ5βDHSH∗ID∗B3,νHmaxSH+δHεHEH+αHRH,dEHdt=1-minmax0,λ6-λ8EH∗+λ7-λ6δHγHEH∗B4,ρHmaxβDHSHID-δHγH+mH+δHεHEH-minmax0,λ6-λ8EH∗+λ7-λ6δHγHEH∗B4,ρHmaxEH,dIHdt=δHγHEH-mH+μHIH,dRHdt=minmax0,λ5-λ8SH∗+λ6-λ5βDHSH∗ID∗B3,νHmaxSH-mH+αHRH+minmax0,λ6-λ8EH∗+λ7-λ6δHγHEH∗B4,ρHmaxEH,dλ1dt,dλ2dt,dλ3dt,dλ4dt,dλ5dt,dλ6dt,dλ7dt,dλ8dt,with,λitf=0,i=1,2,3,4,5,6,7,8.


### 5.1. Numerical Simulations of the Optimality System

To determine the control strategiesνD,ρD,νH, and ρH, as given in the objective functional, we began an iteration of the model until convergence is achieved. The results of the simulation of the control strategies are displayed below. We consider equal weights of (A1=1,A2=1,A3=1,A4=1) for both exposed and infected classes. We varied the cost associated with the objective functional, which indicate that, with low cost of vaccination, the rate at which individuals will seek for vaccination of their susceptible dogs will increase, and this could result in low transmission of rabies in a heterogeneous population. We consider the various cost of preexposure prophylaxis and postexposure prophylaxis to be (B1=1,B2=4,B3=1,B4=4). We found that the optimal time in controlling the infection using preexposure prophylaxis in dogs is much better than using postexposure prophylaxis in dogs, as shown by the trajectories of the red line and blue line in Figure 4, respectively. The blue line in Figure 4 indicates that applying postexposure prophylaxis will considerably take a longer time in controlling of rabies in dogs. The green line in Figure 4 signifies that preexposure prophylaxis in humans increases the immunity levels of humans and hence reduces the rate at which individuals move to the infected stage. Figures 5 and 6 show the effect of using only one control strategy on the model. Therefore, Figure 5(a) shows that applying only postexposure prophylaxis (treatment or quarantine) of dogs has a low positive impact on the model. Figure 5(b) shows that sticking to the use of pre- and postexposure prophylaxis in human without administering pre- and postexposure prophylaxis in the dog population will result in a high of the rabies infection in the human population. Figure 6(a) also shows that combining pre- and postexposure prophylaxis (vaccination and treatment/quarantine) in the dog compartment will reduce the spread of the rabies virus, thereby reducing the using of pre- and postexposure prophylaxis (vaccination and treatment) in humans. Figure 6(b) indicates that a rapid use of pre- and postexposure prophylaxis in the human population will reduce the number of rabies deaths in the human population. Figure 7 shows the simulation effects of applying both controls on the model. Figure 7(a) shows that, with the use of the optimal control strategies, the rate of the infection in the susceptible dogs will reduce significantly. Figures 7(b) and 7(c) show that there is a proportional decrease in the number of exposed and infected dogs when the control measures are applied. Similarly, Figures 7(e) and 7(f) show a significant decrease in the number of infected and exposed humans when the control measures are applied. Figure 7(d) shows that there is a proportional increase in the number of recovered dogs when the control measures are applied. Finally, Figures 8(a)–8(h) show the simulation effect of corresponding adjoint functions.Figure 4
The simulation effect of the controls.Figure 5
The trajectories of the model with and without pre- and postexposure prophylaxis on exposed humans and that of the exposed dogs.
(a)
A plot ofνD=ρD≠0 and νH=ρH=0
(b)
A plot ofνD=ρD=0 and νH=ρH≠0Figure 6
The trajectories of the model with and without pre- and postexposure prophylaxis on infected humans and that of the infected dogs.
(a)
A plot ofνD=ρD≠0 and νH=ρH=0
(b)
A plot ofν=ρ=0 and νH=ρH≠0Figure 7
The trajectories of the model with and without optimal control on individual compartments.
(a)
Susceptible dogs with and without control
(b)
Exposed dogs with and without control
(c)
Infected dogs with and without control
(d)
Recovered dogs with and without control
(e)
Infected humans with and without control
(f)
Exposed humans with and without controlFigure 8
The trajectories of the model with and without optimal control on individual compartments and corresponding adjoint function.
(a)
The cost functionλ1 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(b)
The cost functionλ2 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(c)
The cost functionλ3 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(d)
The cost functionλ4 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(e)
The cost functionλ5 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(f)
The cost functionλ6 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(g)
The cost functionλ7 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(h)
The cost functionλ8 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4


## 5.1. Numerical Simulations of the Optimality System

To determine the control strategiesνD,ρD,νH, and ρH, as given in the objective functional, we began an iteration of the model until convergence is achieved. The results of the simulation of the control strategies are displayed below. We consider equal weights of (A1=1,A2=1,A3=1,A4=1) for both exposed and infected classes. We varied the cost associated with the objective functional, which indicate that, with low cost of vaccination, the rate at which individuals will seek for vaccination of their susceptible dogs will increase, and this could result in low transmission of rabies in a heterogeneous population. We consider the various cost of preexposure prophylaxis and postexposure prophylaxis to be (B1=1,B2=4,B3=1,B4=4). We found that the optimal time in controlling the infection using preexposure prophylaxis in dogs is much better than using postexposure prophylaxis in dogs, as shown by the trajectories of the red line and blue line in Figure 4, respectively. The blue line in Figure 4 indicates that applying postexposure prophylaxis will considerably take a longer time in controlling of rabies in dogs. The green line in Figure 4 signifies that preexposure prophylaxis in humans increases the immunity levels of humans and hence reduces the rate at which individuals move to the infected stage. Figures 5 and 6 show the effect of using only one control strategy on the model. Therefore, Figure 5(a) shows that applying only postexposure prophylaxis (treatment or quarantine) of dogs has a low positive impact on the model. Figure 5(b) shows that sticking to the use of pre- and postexposure prophylaxis in human without administering pre- and postexposure prophylaxis in the dog population will result in a high of the rabies infection in the human population. Figure 6(a) also shows that combining pre- and postexposure prophylaxis (vaccination and treatment/quarantine) in the dog compartment will reduce the spread of the rabies virus, thereby reducing the using of pre- and postexposure prophylaxis (vaccination and treatment) in humans. Figure 6(b) indicates that a rapid use of pre- and postexposure prophylaxis in the human population will reduce the number of rabies deaths in the human population. Figure 7 shows the simulation effects of applying both controls on the model. Figure 7(a) shows that, with the use of the optimal control strategies, the rate of the infection in the susceptible dogs will reduce significantly. Figures 7(b) and 7(c) show that there is a proportional decrease in the number of exposed and infected dogs when the control measures are applied. Similarly, Figures 7(e) and 7(f) show a significant decrease in the number of infected and exposed humans when the control measures are applied. Figure 7(d) shows that there is a proportional increase in the number of recovered dogs when the control measures are applied. Finally, Figures 8(a)–8(h) show the simulation effect of corresponding adjoint functions.Figure 4
The simulation effect of the controls.Figure 5
The trajectories of the model with and without pre- and postexposure prophylaxis on exposed humans and that of the exposed dogs.
(a)
A plot ofνD=ρD≠0 and νH=ρH=0
(b)
A plot ofνD=ρD=0 and νH=ρH≠0Figure 6
The trajectories of the model with and without pre- and postexposure prophylaxis on infected humans and that of the infected dogs.
(a)
A plot ofνD=ρD≠0 and νH=ρH=0
(b)
A plot ofν=ρ=0 and νH=ρH≠0Figure 7
The trajectories of the model with and without optimal control on individual compartments.
(a)
Susceptible dogs with and without control
(b)
Exposed dogs with and without control
(c)
Infected dogs with and without control
(d)
Recovered dogs with and without control
(e)
Infected humans with and without control
(f)
Exposed humans with and without controlFigure 8
The trajectories of the model with and without optimal control on individual compartments and corresponding adjoint function.
(a)
The cost functionλ1 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(b)
The cost functionλ2 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(c)
The cost functionλ3 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(d)
The cost functionλ4 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(e)
The cost functionλ5 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(f)
The cost functionλ6 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(g)
The cost functionλ7 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4
(h)
The cost functionλ8 for A1=A2=A3=A4=1 and B1=1,B2=4,B3=1,B4=4


## 6. Discussion

The numerical simulations of the resulting optimality system show that, during the case where it is more expensive to vaccinate than treatment, more resources should be invested in treating affected individuals until the disease prevalence begins to fall. This option, however, does not reduce the number of individuals expose to the disease quickly enough, thus resulting in an overall increase in the infected human population. On the other hand, if it is more expensive to treat than to vaccinate, then more susceptible dogs should be vaccinated, so as to lower the rate at which newborn dogs get infected. Nevertheless, in the case where both measures are equally expensive, the simulation shows that the optimal way to drive the epidemic towards eradication within any specified period is to use more preexposure prophylaxis in both compartments.


## 7. Conclusion

We studied an optimal control model of rabies transmission dynamics in dogs and the best way of reducing death rate of rabies in humans. The stability analysis shows that the disease-free equilibrium is locally and globally asymptotically stable. We also obtained an optimal control solution for the model which predicts that the optimal way of eliminating deaths from canine rabies as projected by the global alliance for rabies control [30] is using more of preexposure prophylaxis in both dogs and humans and public education; however, the results show that the effective and optimal consideration of preexposure prophylaxis and postexposure prophylaxis in humans without an optimal use of vaccination in the dog population is not beneficial if total elimination of the disease is desirable in Africa and Asia. Any combination strategy which involves vaccination in the dogs’ population gives a better result and hence it may be beneficial in eliminating the disease in Asia, Africa, and Latin America.



---

*Source: 2451237-2017-07-16.xml*