# Mathematical Analysis of Rabies Infection

Authors: C. S. Bornaa; Baba Seidu; M. I. Daabo
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

A mathematical model is proposed to study the dynamics of the transmission of rabies, incorporating predation of dogs by humans. The model is shown to have a unique disease-free equilibrium which is globally asymptotically stable wheneverℛ0≤1. Local sensitivity analysis suggests that the disease can be controlled through reducing contact with infected dogs, increasing immunization of dogs, screening recruited dogs, culling of infected dogs, and use of dog meat as a delicacy.


---

## Body

## 1. Introduction

Rabies is a zoonotic viral disease that is usually transmitted by an infected animal (such as foxes, raccoons, cats, coyotes, bats, skunks, and dogs) through bites or scratch that introduces the virus into the blood of another animal or a human [1, 2]. Once the virus enters the body, it either migrates directly to the brain via the peripheral nervous system to replicate there or stays in the muscles to replicate before migrating to the brain via the neuromuscular junctions. The infection becomes fatal upon reaching the brain since it produces acute inflammations of the brain, leading to coma and eventually death. The infection usually goes through five distinct phases of incubation (3 to 12 weeks), prodome (2 to 10 days, with worsening symptoms over time), acute neurologic period, coma, and death. Early symptoms of rabies infection are general body weakness, fever, and headache, which are also peculiar with common flu and other viral diseases [3]. Rabies is preventable and treatment is possible at early stages of infection. Antirabies vaccination of all pets and domestic animals and humans, screening of imported animals, and awareness creation are some preexposure prophylaxis strategies that governments can embark on to prevent rabies. Individuals can also prevent the disease by vaccinating their pets and keeping away from possible infection sources. Postexposure prophylaxis strategies including washing the bitten area with soapy water for some time and administering a series of shots of rabies immune globulin at the early stages are often effective.All mammals are possible reservoirs of the rabies virus, but small rodents such as rats and rabbits are less likely to cause transfers. A survey conducted in Ghana by Addy [4] revealed that, of 1,514 exposed domestic animals examined, dogs, cats, and cattle are responsible for 98.7%, 0.07%, and 0.70% infections, respectively. Rural areas of Southeast Asia and Africa bear the brunt of rabies with tens of thousands contracting the disease annually. Most of these infections are caused by dog bites due to the fact that dogs are closest in contact with humans since they are used as pets, as guards, for hunting purposes, and even as delicacies.The impact of mathematics continues to increase over the years as more evidence is found to support the fact that mathematical modeling helps to increase our understanding of the dynamics of infectious diseases. This has led to increased use of mathematics to model many infectious disease (see [5–18], and references therein).Zoonotic diseases have received their fair share of attention by mathematicians. Specifically, Wang and Lou [19] developed an ODE model to study the impact of combinations of culling and vaccination on the control of rabies transmission and observed that vaccination alone was a better strategy whilst culling was worst in the control of rabies transmission. The work of Wang and Lou [19] is however contradicted by that of Carroll et al. [20] that suggested that culling is more effective than vaccination, arguing that adding immunocontraceptives on animal populations could help improve the effectiveness of vaccination. To analyze the strategies for optimal distribution of vaccine baits in order to minimize the spread of the disease and the cost of carrying out the control, Ding et al. [21] formulated a model to describe rabies in raccoons with discrete time and spatial features. Hou et al. [12] suggested that the control of rabies needs to include awareness creation about rabies, increased domestic dog vaccination, and reduction in the number of stray dogs. Asamoah et al. [5] also suggested that vaccination of pets and use of preexposure and postexposure prophylaxis could help in the control of rabies spread. In China and some parts of Africa (for example, the Upper East and West regions of Ghana), dog meat is a delicacy. It is noteworthy that while so much research has been carried out on the transmission and control of rabies, the impact of predation of dogs by humans has not been studied in the transmission of rabies. This current work seeks to incorporate consumption of dog meat by humans into a dog-human rabies model.The rest of the paper is organized as follows. In Section2, the mathematical model under consideration is derived, while Section 3 discusses some basic qualitative properties of the model. In Section 4, numerical simulation of the model is carried out to study the impact of various factors on the transmission of rabies. The findings and conclusions are presented in Section 5.


## 2. Model Formulation

We consider the transmission of rabies from a dog population to a dog-predator human population. Each of the dog and human populations is divided into three compartments of Susceptibles (S), Exposed (E), and Infected (I). We use subscripts d and h to denote dogs and humans, respectively, so that Susceptible, Exposed, and Infected dogs are represented by Sd, Ed, and Id, respectively, while Susceptible, Exposed, and Infected humans are represented by Sh, Eh, and Ih, respectively. The total dog and human populations are therefore given by Wdt=Sdt+Edt+Idt and Wht=Sht+Sht+Iht. Human-to-human transmission of the rabies virus is assumed to be nonexistent or insignificant so that only direct dog-to-dog and dog-to-human transmission is considered. We assume that all recruitments of dogs and humans are made into the susceptible classes Sdt and Sht, at constant rates Pd and Ph, respectively. Dogs and humans contract the disease through effective contact with infected dogs, with transmission probabilities βd and βh, respectively. Rabies transmission is assumed to be reduced through vaccination of susceptible dogs and humans at rates δd and δh, respectively. Upon effective contact with infected dogs, dogs and humans are moved to the exposed classes before being moved (if postexposure prophylaxis is ineffective) to the infected classes at rates βDD and βH for dogs and humans, respectively. A successful treatment at rate, θ, of exposed and infected humans is assumed while successful culling at rate Cd is applied to exposed and infected dogs. Dogs and humans are assumed to suffer rabies-induced death at rates σd and σh, respectively. Parameters μd and μh are taken to be the natural death rates in dogs and humans, respectively. We assume that dogs suffer predation from humans at rate γ1. Since dogs in the infected class are assumed to be culled, we exclude the infected dog population from the predation.With these assumptions, the mathematical model describing the dynamics of dog-human rabies transmission is represented by the set of differential equations in equation (1). Table 1 shows the descriptions of the parameters used in the model.Table 1
Description of parameters of the model.
Parameter
Description
P
d
Recruitment into the susceptible dog population (birth and immigration)
P
h
Recruitment into the susceptible human population (birth and immigration)
γ
1
Rate of human predation of dogs
γ
≪
γ
1
Rate at which consumed dogs are converted into susceptible humans
β
d
Rate of disease transmission within the dog population
β
D
D
<
σ
d
Rate at which exposed dogs develop clinical rabies
β
h
Rate of dog-to-human disease transmission
β
H
<
σ
h
Rate at which exposed humans develop clinical rabies
μ
d
<
β
D
D
Natural death rate of dogs
μ
h
Natural death rate of humans
C
d
Death rate of dogs due to culling
σ
d
Disease induce death rate within the dog population
σ
h
Disease-induced death rate within the human population
δ
d
Rate of immunization of dogs
δ
h
Rate of immunization of humans
θ
Rate of treatment of infected human population(1)
d
S
d
d
t
=
P
d
−
β
d
S
d
E
d
+
I
d
1
−
δ
d
−
μ
d
+
γ
1
S
h
+
E
h
S
d
,
d
E
d
d
t
=
β
d
S
d
E
d
+
I
d
1
−
δ
d
−
μ
d
+
γ
1
S
h
+
E
h
+
C
d
+
β
D
D
E
d
,
d
I
d
d
t
=
β
D
D
E
d
−
μ
d
+
γ
1
S
h
+
E
h
+
σ
d
+
C
d
I
d
,
d
S
h
d
t
=
P
h
+
γ
S
h
+
E
h
S
d
+
E
d
+
I
d
+
θ
I
h
+
E
h
−
β
h
I
d
+
E
d
1
−
δ
h
+
μ
h
S
h
,
d
E
h
d
t
=
β
h
S
h
I
d
+
E
d
1
−
δ
h
−
θ
+
μ
h
+
β
H
E
h
,
d
I
h
d
t
=
β
H
E
h
−
μ
h
+
σ
h
+
θ
I
h
,
with initial conditions
S
d
0
≥
0
,
E
d
0
≥
0
,
I
d
0
≥
0
,
S
h
0
≥
0
,
E
h
0
≥
0
,
I
h
0
≥
0.For convenience, we make the following substitutions(2)ψ1=βd1−δd,ψ2=βh1−δh,ψ3=μd+Cd+βDDψ4=μd+σd+CCd,ψ5=θ+μh+βH,ψ6=μh+σh+θ.In the next section, we discuss some basic properties of the model.


## 3. Basic Properties of the Model

### 3.1. Positivity of Solutions

Model (1) is an epidemiological model, and hence, it is necessary that the associated population sizes be positive. Model (1) should be considered in a feasible region where such property (nonnegative) is preserved. This is provided in Theorem 1.Theorem 1.
If positive conditions initial conditions are provided for (1), then all its solutions remain positive for t>0.Proof.
From the first equation, if at some pointt>0, we have Sdt=0 and then, we have dSd/dt=Pd>0. This shows that Sdt>0.
Similar arguments can be used to show that the other state variables have nonnegative solutions for allt>0, hence completing the proof.Theorem1 indicates that the model is epidemiologically meaningful and mathematically well posed.


### 3.2. Boundedness of Solution of the Model

Theorem 2 (Boundedness of the Model).
All solutions of system (1) starting in ℝ+6 are bounded.Proof.
ConsiderWdt=Sdt+Edt+Idt and Wht=Sht+Eht+Iht.
Then,(3)dWddt=Pd−μd+γ1Sh+EhSd−μd+γ1Sh+Eh+CdEd−μd+γ1Sh+Eh+σdId=Pd−μd+γ1Sh+EhSd+Ed+Id−CdEd−σdId,dWhdt=Ph+γSh+EhSd+Ed+Id−μhSh−μhEh−Ihμh+σh=Ph+γSh+EhSd+Ed+Id−Sh+Eh+Ihμh−Ihσh.
From the above equations, we have(4)dWddt≤Pd−γ1Sh+EhSd+Ed+Id−μdWd,dWhdt≤Ph+γSh+EhSd+Ed+Id−μhWh.
If we setW=Wh+Wd and P=Pd+Ph, then,(5)dWdt≤P−minμd,μhW,so that(6)Wt≤W0e−Pt+Pminμd,μh1−e−Pt.
The feasible region for the model is thus given by(7)Ω=Sd,Ed,Id,Sh,Eh,Ih∈ℝ6:0≤W≤Pminμd,μh.
System (1) can therefore be conveniently studied in Ω.


### 3.3. Equilibrium Points of the Model

The model exhibits a biologically reasonable disease-free equilibriumε0=Sd0,0,0,Sh0,0,0, where Sd0=Pd/μd+γ1Sh0 and Sh0 satisfy the following quadratic equation:(8)γ1μhSh02−Pdγ+Phγ1−μhμdSh0−μdPh=0.Hence,(9)Sh0=Pdγ+Phγ1−μhμd+Pdγ+Phγ1−μhμd2+4γ1μhμdPh2γ1μh,Sd0=2μhPdPdγ+Phγ1+μhμd+Pdγ+Phγ1−μhμd2+4γ1μhμdPh.Remark 1.
Model (1) has a unique realistic DFE.Proof.
Since the coefficient ofSh02 and the constant term in (8) have different signs, Descartes’ rule of sign shows that equation (8) has only one positive solution. This concludes the proof of the remark.
The basic reproduction ratio,ℛ0, is the average number of secondary infections that are caused by one infectious individual, which is introduced into an initially disease-free population, during its infectious period. Using the next-generation method of Diekmann et al. [22], we compute ℛ0 as(10)ℛ0=Sd0ψ1γ1Sh0+βDD+ψ4γ1Sh0+ψ4γ1Sh0+ψ3.
In the presence of Infectives, model (1) is said to be exhibiting an endemic equilibrium point ε∗=Sd∗,Ed∗,Id∗,Sh∗,Eh∗,Ih∗, where Sd∗,Ed∗,Id∗, and Ih∗ are given by(11)Sd∗=ψ3+γ1Sh∗+Eh∗ψ4+γ1Sh∗+Eh∗ψ1ψ4+γ1Sh∗+Eh∗+βDD,Ed∗=ψ5Eh∗ψ4+γ1Sh∗+Eh∗ψ2Sh∗ψ4+γ1Sh∗+Eh∗+βDD,Id∗=βDDψ5Eh∗ψ2Sh∗ψ4+γ1Sh∗+Eh∗+βDDIh∗=βHEh∗ψ6,and Sh∗ and Eh∗ satisfy the following set of equations (see Appendix for derivation):(12)ShgEhT+b1TEhT=0,ShℋEhT+b2TEhT=0.
Here,(13)Sh=Sh∗4,Sh∗3,Sh∗2,Sh∗,Eh=Eh∗3,Eh∗2,Eh∗,1,g=g11000G213g1100g312g21+γ12ψ5ψ13g110g41g31+γ1ψ3+ψ4ψ1ψ5g21+2γ12ψ5ψ1G11,b1=0ψ5ψ1ψ3ψ4γ1ψ1ψ5ψ3+ψ4γ12ψ5ψ1,ℋ=ℋ11000ℋ213ℋ1100ℋ31ℋ21+ε13ℋ110ℋ41ℋ31+ε3ℋ21+2ε1+ε2ℋ11,b2=00γψ6ψ5ψ1ψ4+βDDγψ6ψ5ψ1γ1,where(14)g11=ψ2γ13,g31=ψ2γ1ψ3+ψ4μd−Pdψ1+ψ3ψ4,ℋ11=γψ6ψ2γ12,ℋ41=Phψ1ψ4+βDDψ2ψ6,ℋ31=ψ2ψ6Phψ1γ1+γψ3ψ4−ψ1ψ2ψ6μhψ4+βDD,ε1=γψ6ψ2ψ3+ψ2ψ4+ψ1ψ5+θψ1ψ2ψ6+βH−ψ1ψ2ψ6ψ5γ1,ε3=ψ1ψ4+βDDθψ2ψ6+βH+ψ6ψ5γ−ψ2+ψ2ψ6μh,g21=γ12ψ2ψ3+ψ4+μd,g41=ψ2ψ3ψ4μd−Pdψ1ψ4+βDD,ℋ21=ψ2ψ6γ1γψ3+ψ4−μhψ1,where Sh∗ and Eh∗ can be found by solving equation (12). Even though the solution can be done numerically, the nature of the problem suggests that an algebraic solution may be obtainable. The problem falls in the category of Open Problem 1.Open Problem 1.
Letb=b1,b2,…,bn,c=c1,c2,…,cn, x=xn,xn−1,…x,y=yn−1,xn−2,…y,1, and A and B be n−dimensional lower triangular square matrices. Is there an algebraic technique that be used to determine x and y satisfying the following equations?(15)xAyT+byT=0,xByT+cyT=0.
In the next section, the local stability of the disease-free equilibrium is discussed.


### 3.4. Local Stability of the Disease-Free Equilibrium Point

We use the Lyapunov indirect method to study the local stability of the disease-free equilibrium point of the model. An equilibrium point is said to be locally asymptotically stable if all eigenvalues of the Jacobian model evaluated at the given equilibrium point have negative real parts.The Jacobian matrix of the model evaluated atε0 is given by(16)Jε0=−μd−γ1Sh0−ψ1Sd0−ψ1Sd0−γ1Sd0−γ1Sd000ψ1Sd0−ψ3−γ1Sh0ψ1Sd00000βDD−ψ4−γ1Sh0000γSh0γ−ψ2Sh0γ−ψ2Sh0γSd0−μhγSd0+θθ0ψ2Sh0ψ2Sh00−ψ500000βH−ψ6.Two of the eigenvalues (−ψ5 and −ψ6) of Jε0 are clearly negative and the remaining eigenvalues are the zeros of P1λP2λ=0, where(17)P1λ=λ2+2γ1Sh0+ψ3+ψ4−ψ1Sd0λ+ψ4+γ1Sh0ψ3+γ1Sh01−ℛ0=0,P2λ=λ2+μh−γSd0+γ1Sh0+μdλ+μhγ1Sh0+μhμd−γSd0μd=0.The second equation of (8) has both zeros being negative if (remembering thatSd0=Pd/μd+γ1Sh0)(18)μhγ1Sh0+μd−γPdμdμd+γ1Sh0>0,γ1Sh0+μd+μh−Pdγμd+γ1Sh0>0.These conditions are easy to establish using (8).Also, the first equation in (17) has zeros with negative real parts if(19)ψ4+γ1Sh0ψ3+γ1Sh01−ℛ0>0,2γ1Sh0+ψ3+ψ4−ψ1Sd0>0.The first condition is true wheneverℛ0<1 whilst the second condition can be written as ψ4+ψ3+2γ1Sh0ψ4+γ1Sh0+βDD/ψ4+γ1Sh0ψ3+γ1Sh0−ℛ0>0, which holds if ℛ0<1. The following result is therefore established.Theorem 3.
The disease-free equilibrium state,ε0, of model (1) is locally asymptotically stable whenever ℛ0<1 and unstable otherwise.


### 3.5. Local Stability Analysis of Endemic Equilibrium

The stability of the endemic equilibriums of the model is similarly established by finding the eigenvalues of the Jacobian matrix of the model evaluated at the positive solutions of (12). The Jacobian model evaluated at ε∗ if it exists is given by(20)Jε∗=J11−Sd∗ψ1−Sd∗ψ1−γ1Sd∗−γ1Sd∗0ψ1EIdJ22Sd∗ψ1−γ1Ed∗−γ1Ed∗00βDDJ33−γ1Id∗−γ1Id∗0γSEhγSEh−ψ2Sh∗γSEh−ψ2Sh∗J44γSEId+θθ0ψ2Sh∗ψ2Sh∗ψ2EId−ψ500000βH−ψ6,where(21)J11=−ψ1EId+SEhγ1+μd,J33=−SEhγ1+ψ4,EId=Ed∗+Id∗,SEId=Sd∗+Ed∗+Id∗,J22=Sdψ1−SEhγ1−ψ3,J44=γSEId−ψ2EId−μh,SEh=Sh∗+Eh∗.The Jacobian matrixJε∗ has a characteristic polynomial given by(22)∑n=06Λnλn=0,where(23)Λ6=1,Λ5=ψ6+ψ5+γEhSEIdSh+PdSd+θEh+Ih+PhSh+Sdψ1IdEd+βDDEdId,Λ4=ψ6Sd+PhShSdPd+ψ6Sh+ψ5ShPh+ψ5+γSEIdEhSh+θEh+IhSh+βDDEdId+Sdψ1IdEdPdSd+ψ6+θIhSh+βDDEdId+Sdψ1IdEdψ5+γγ1SEIdSEh+EIdψ12Sd+γSEIdEhSh+PhSh+θEh+IhShSdψ1IdEd+βDDEdId,Λ3=γSEhγ1+EIdψ12SdEdβDDEd+ψ5PhSh+βDDEdId+Sdψ1IdEd+PdSd+γSEIdEhSh+θEhShβDDEdId+Sdψ1IdEd+PdSdψ6+γγ1SEIdSEh+EIdψ12Sd+PhSh+θIhShβDDEdId+Sdψ1IdEd+PdSd+βDDEdPdSdId+ψ1IdPdEdψ6+ψ5+Sdψ1SEhId+SEhEIdPdSd+Ed+SdβDDEdId+Id+SdSdψ1IdEdSEhγγ1+Sdψ12γSEIdEhSh+PhSh+θEh+IhSh+ψ12SdβDDEdIdEId+βDDEdPdSdId+ψ1IdPdEdγSEIdEhSh+PhSh+θEh+IhSh+γEIdSEIdSEhγ1ψ2.Clearly,Λn>0,∀n=3,4,5,6, and Λ0, Λ1, and Λ2 can similarly be shown to be positive. By Descartes’ rule of signs, the following result is thus established.Theorem 4.
The endemic equilibriumε∗ (whenever it exists) is locally asymptotically stable.


### 3.6. Global Stability of Disease-Free Equilibrium Point

To study the global stability ofℰ0, we employ the technique of Castillo-Chavez et al. [23] as follows:IfXG=Sd,Sh and ZG=Ed,Id,Eh,Ih, then model (1) can be written as(24)dXGdt=FXG,ZG,dZGdt=GXG,ZG,GXG,0=0.By the theorem of Castillo-Chavez et al. [23], ε0 is globally asymptotically stable if the following conditions are satisfied:GS1:ε0 is locally asymptotically stableGS2:XG∗ is globally stable for dXG/dt=GX,0GS3:GXG,ZG=AGZG−G^XG,ZG, G^XG,ZG≥0, where AG=DZGXG∗,0 is an M-MatrixThe first condition GS1 has been established forℛ0<1 in Theorem 3.For the second condition, we note thatFXG,0=Pd−μdSd−γ1ShSdPh+γShSd−μhSh is a limiting function of dXG/dt=FXG,ZG. Therefore, the second condition is satisfied.As for the third condition, we observe thatAG is given by(25)AG=Sd0ψ1−ψ3−γ1Sh0Sd0ψ100βDD−ψ4+γ1Sh000Sh0ψ2Sh0ψ2−ψ5000βH−ψ6,which is an M−matrix (i.e., all nondiagonal entries are nonnegative).Furthermore, the matrixG^XG,ZG is given by(26)G^XG,ZG=ψ1Sd0−SdEd+Id+Edγ1Sh+Eh−Sh0Idγ1Sh+Eh−Sh0ψ2Sh0−ShEd+Id0.FromdSd/dt, we have(27)Sdt=ϑt0Sdt0+Pdtϑt,where(28)ϑt=exp∫μd+γ1Sh+Eh+Ed+Idψ1dt.Hence,Sd0≥Sdt. Similarly, Sht≥Sh0 and hence, G^≥0, satisfies condition GS3. Now, since all three conditions are satisfied, the following result is established.Theorem 5.
The disease-free equilibrium point,ε0, is globally asymptotically stable whenever ℛ0≤1.For simulation purposes, we use the following parameter values:(29)Pd=0.45,Ph=0.1,βd=0.062,δd=0.5βDD=0.6,σd=0.8,βh=0.0229,δh=0.005,βH=0.0001,σh=0.0001,γ1=0.01,γ=0.0001,μd=0.06,Cd=0.05,θ=0.9,μh=0.04.


### 3.7. Sensitivity Analysis

To determine the impact of model parameters on the spread of rabies, we employed global sensitivity analysis technique of Marino et al. [24] to calculate the partial rank correlation coefficients of the model parameters on which the basic reproduction number depends. The result of sensitivity analysis informs which model parameters are most important in determining the output variable (in this case, ℛ0) and hence should be given much attention in terms of measurement for accuracy. The PRCC of results are shown in Figure 1.Figure 1
Plot of PRCC ofℛ0 with respect to model parameters.Local sensitivity analysis discusses the impact of small changes in parameter values on the spread of the disease. To study the sensitivity of the disease spread to parameter changes, the basic reproduction numberℛ0 is used, since it determines the persistence or possible eradication of the disease. The normalized forward sensitivity index is used to perform the analysis. It is defined as follows.Letℛ0 be a differentiable function xi. Then, the normalized forward sensitivity index of R0 relative to xi is given by(30)ηℛ0xi=∂ℛ0∂xi×xiℛ0.Local sensitivity indexes are presented in Table2. This index measures the relative percentage change in ℛ0 due to a percentage change in xi. The sensitivity indexes of the parameters determining ℛ0 are calculated and presented in Table 2. A positive sensitivity index of a parameter implies that an increase (decrease) in the value of the parameter will result in an increase (decrease) in ℛ0. On the other hand, a negative sensitivity index implies that an increase (decrease) in the value of the parameter will result in a decrease (an increase) in ℛ0.Table 2
Sensitivity indexes ofℛ0 with respect to model parameters.
C
d
P
d
P
h
β
D
D
β
d
δ
d
μ
d
μ
h
σ
d
γ
γ
1
−0.009
−2.010
−0.0143
−0.012
1.000
−1.000
1.056
3.024
−0.011
−3.967
0.943


### 3.8. Bifurcation Analysis

To study the stability of endemic equilibrium points, the center manifold theory as described in Theorem 4.1 of [25] is used as an alternative to the Lyapunov indirect method.Choosingβd=βd∗=ψ3+γ1Sh0ψ4+γ1Sh0/1−δdSd0γ1Sh0+σd+ψ3 as a bifurcation parameter leads the Jacobian model to have a zero simple eigenvalue. The right and left eigenvectors associated with this zero eigenvalue are given by w=w2η1,1,η3,η4,η5,η6T and v=0,v2,ζ3v2,0,0,0 with η3ζ3+1w2v2=1, where(31)η1=Sd0φ1γ1−φ2γ+μhφ2Sd0μd−μhμd+γ1Sh0,η3=βDDψ4+γ1Sh0,η4=φ2γSh0−φ1μd+γ1Sh0Sd0μd−μhμd+γ1Sh0,η5=ψ2Sh0γ1Sh0+ψ4+βDDψ4+γ1Sh0ψ5,η6=βHψ2Sh0γ1Sh0+ψ4+βDDψ4+γ1Sh0ψ5ψ6,ζ3=γ1Sh0+ψ3γ1Sh0+ψ4+βDD,φ1=Sd0ψ11+η3+γ1η5,φ2=Sh01+η3γ−ψ2+γSd0+θη5+θη6.The bifurcation coefficients are thus given by(32)a=2w2p∗−1,b=Sd01−δd1+η31+η3ζ3,where p∗=βd∗η11−δd1+η3/γ1η4+η51+η3ζ3.Clearly,b is always positive and the sign of a is dependent on p∗. The following result therefore easily follows [25].Theorem 6.
Model (1) undergoes a backward bifurcation at ℛ0=1 when p∗>1 and a forward bifurcation when p∗<1 in which case the endemic equilibrium will be locally asymptotically stable.We note that sinceℰ0 is GAS whenever ℛ0≤1, then backward bifurcation may only occur whenever ℛ approaches unity from above.


## 3.1. Positivity of Solutions

Model (1) is an epidemiological model, and hence, it is necessary that the associated population sizes be positive. Model (1) should be considered in a feasible region where such property (nonnegative) is preserved. This is provided in Theorem 1.Theorem 1.
If positive conditions initial conditions are provided for (1), then all its solutions remain positive for t>0.Proof.
From the first equation, if at some pointt>0, we have Sdt=0 and then, we have dSd/dt=Pd>0. This shows that Sdt>0.
Similar arguments can be used to show that the other state variables have nonnegative solutions for allt>0, hence completing the proof.Theorem1 indicates that the model is epidemiologically meaningful and mathematically well posed.


## 3.2. Boundedness of Solution of the Model

Theorem 2 (Boundedness of the Model).
All solutions of system (1) starting in ℝ+6 are bounded.Proof.
ConsiderWdt=Sdt+Edt+Idt and Wht=Sht+Eht+Iht.
Then,(3)dWddt=Pd−μd+γ1Sh+EhSd−μd+γ1Sh+Eh+CdEd−μd+γ1Sh+Eh+σdId=Pd−μd+γ1Sh+EhSd+Ed+Id−CdEd−σdId,dWhdt=Ph+γSh+EhSd+Ed+Id−μhSh−μhEh−Ihμh+σh=Ph+γSh+EhSd+Ed+Id−Sh+Eh+Ihμh−Ihσh.
From the above equations, we have(4)dWddt≤Pd−γ1Sh+EhSd+Ed+Id−μdWd,dWhdt≤Ph+γSh+EhSd+Ed+Id−μhWh.
If we setW=Wh+Wd and P=Pd+Ph, then,(5)dWdt≤P−minμd,μhW,so that(6)Wt≤W0e−Pt+Pminμd,μh1−e−Pt.
The feasible region for the model is thus given by(7)Ω=Sd,Ed,Id,Sh,Eh,Ih∈ℝ6:0≤W≤Pminμd,μh.
System (1) can therefore be conveniently studied in Ω.


## 3.3. Equilibrium Points of the Model

The model exhibits a biologically reasonable disease-free equilibriumε0=Sd0,0,0,Sh0,0,0, where Sd0=Pd/μd+γ1Sh0 and Sh0 satisfy the following quadratic equation:(8)γ1μhSh02−Pdγ+Phγ1−μhμdSh0−μdPh=0.Hence,(9)Sh0=Pdγ+Phγ1−μhμd+Pdγ+Phγ1−μhμd2+4γ1μhμdPh2γ1μh,Sd0=2μhPdPdγ+Phγ1+μhμd+Pdγ+Phγ1−μhμd2+4γ1μhμdPh.Remark 1.
Model (1) has a unique realistic DFE.Proof.
Since the coefficient ofSh02 and the constant term in (8) have different signs, Descartes’ rule of sign shows that equation (8) has only one positive solution. This concludes the proof of the remark.
The basic reproduction ratio,ℛ0, is the average number of secondary infections that are caused by one infectious individual, which is introduced into an initially disease-free population, during its infectious period. Using the next-generation method of Diekmann et al. [22], we compute ℛ0 as(10)ℛ0=Sd0ψ1γ1Sh0+βDD+ψ4γ1Sh0+ψ4γ1Sh0+ψ3.
In the presence of Infectives, model (1) is said to be exhibiting an endemic equilibrium point ε∗=Sd∗,Ed∗,Id∗,Sh∗,Eh∗,Ih∗, where Sd∗,Ed∗,Id∗, and Ih∗ are given by(11)Sd∗=ψ3+γ1Sh∗+Eh∗ψ4+γ1Sh∗+Eh∗ψ1ψ4+γ1Sh∗+Eh∗+βDD,Ed∗=ψ5Eh∗ψ4+γ1Sh∗+Eh∗ψ2Sh∗ψ4+γ1Sh∗+Eh∗+βDD,Id∗=βDDψ5Eh∗ψ2Sh∗ψ4+γ1Sh∗+Eh∗+βDDIh∗=βHEh∗ψ6,and Sh∗ and Eh∗ satisfy the following set of equations (see Appendix for derivation):(12)ShgEhT+b1TEhT=0,ShℋEhT+b2TEhT=0.
Here,(13)Sh=Sh∗4,Sh∗3,Sh∗2,Sh∗,Eh=Eh∗3,Eh∗2,Eh∗,1,g=g11000G213g1100g312g21+γ12ψ5ψ13g110g41g31+γ1ψ3+ψ4ψ1ψ5g21+2γ12ψ5ψ1G11,b1=0ψ5ψ1ψ3ψ4γ1ψ1ψ5ψ3+ψ4γ12ψ5ψ1,ℋ=ℋ11000ℋ213ℋ1100ℋ31ℋ21+ε13ℋ110ℋ41ℋ31+ε3ℋ21+2ε1+ε2ℋ11,b2=00γψ6ψ5ψ1ψ4+βDDγψ6ψ5ψ1γ1,where(14)g11=ψ2γ13,g31=ψ2γ1ψ3+ψ4μd−Pdψ1+ψ3ψ4,ℋ11=γψ6ψ2γ12,ℋ41=Phψ1ψ4+βDDψ2ψ6,ℋ31=ψ2ψ6Phψ1γ1+γψ3ψ4−ψ1ψ2ψ6μhψ4+βDD,ε1=γψ6ψ2ψ3+ψ2ψ4+ψ1ψ5+θψ1ψ2ψ6+βH−ψ1ψ2ψ6ψ5γ1,ε3=ψ1ψ4+βDDθψ2ψ6+βH+ψ6ψ5γ−ψ2+ψ2ψ6μh,g21=γ12ψ2ψ3+ψ4+μd,g41=ψ2ψ3ψ4μd−Pdψ1ψ4+βDD,ℋ21=ψ2ψ6γ1γψ3+ψ4−μhψ1,where Sh∗ and Eh∗ can be found by solving equation (12). Even though the solution can be done numerically, the nature of the problem suggests that an algebraic solution may be obtainable. The problem falls in the category of Open Problem 1.Open Problem 1.
Letb=b1,b2,…,bn,c=c1,c2,…,cn, x=xn,xn−1,…x,y=yn−1,xn−2,…y,1, and A and B be n−dimensional lower triangular square matrices. Is there an algebraic technique that be used to determine x and y satisfying the following equations?(15)xAyT+byT=0,xByT+cyT=0.
In the next section, the local stability of the disease-free equilibrium is discussed.


## 3.4. Local Stability of the Disease-Free Equilibrium Point

We use the Lyapunov indirect method to study the local stability of the disease-free equilibrium point of the model. An equilibrium point is said to be locally asymptotically stable if all eigenvalues of the Jacobian model evaluated at the given equilibrium point have negative real parts.The Jacobian matrix of the model evaluated atε0 is given by(16)Jε0=−μd−γ1Sh0−ψ1Sd0−ψ1Sd0−γ1Sd0−γ1Sd000ψ1Sd0−ψ3−γ1Sh0ψ1Sd00000βDD−ψ4−γ1Sh0000γSh0γ−ψ2Sh0γ−ψ2Sh0γSd0−μhγSd0+θθ0ψ2Sh0ψ2Sh00−ψ500000βH−ψ6.Two of the eigenvalues (−ψ5 and −ψ6) of Jε0 are clearly negative and the remaining eigenvalues are the zeros of P1λP2λ=0, where(17)P1λ=λ2+2γ1Sh0+ψ3+ψ4−ψ1Sd0λ+ψ4+γ1Sh0ψ3+γ1Sh01−ℛ0=0,P2λ=λ2+μh−γSd0+γ1Sh0+μdλ+μhγ1Sh0+μhμd−γSd0μd=0.The second equation of (8) has both zeros being negative if (remembering thatSd0=Pd/μd+γ1Sh0)(18)μhγ1Sh0+μd−γPdμdμd+γ1Sh0>0,γ1Sh0+μd+μh−Pdγμd+γ1Sh0>0.These conditions are easy to establish using (8).Also, the first equation in (17) has zeros with negative real parts if(19)ψ4+γ1Sh0ψ3+γ1Sh01−ℛ0>0,2γ1Sh0+ψ3+ψ4−ψ1Sd0>0.The first condition is true wheneverℛ0<1 whilst the second condition can be written as ψ4+ψ3+2γ1Sh0ψ4+γ1Sh0+βDD/ψ4+γ1Sh0ψ3+γ1Sh0−ℛ0>0, which holds if ℛ0<1. The following result is therefore established.Theorem 3.
The disease-free equilibrium state,ε0, of model (1) is locally asymptotically stable whenever ℛ0<1 and unstable otherwise.


## 3.5. Local Stability Analysis of Endemic Equilibrium

The stability of the endemic equilibriums of the model is similarly established by finding the eigenvalues of the Jacobian matrix of the model evaluated at the positive solutions of (12). The Jacobian model evaluated at ε∗ if it exists is given by(20)Jε∗=J11−Sd∗ψ1−Sd∗ψ1−γ1Sd∗−γ1Sd∗0ψ1EIdJ22Sd∗ψ1−γ1Ed∗−γ1Ed∗00βDDJ33−γ1Id∗−γ1Id∗0γSEhγSEh−ψ2Sh∗γSEh−ψ2Sh∗J44γSEId+θθ0ψ2Sh∗ψ2Sh∗ψ2EId−ψ500000βH−ψ6,where(21)J11=−ψ1EId+SEhγ1+μd,J33=−SEhγ1+ψ4,EId=Ed∗+Id∗,SEId=Sd∗+Ed∗+Id∗,J22=Sdψ1−SEhγ1−ψ3,J44=γSEId−ψ2EId−μh,SEh=Sh∗+Eh∗.The Jacobian matrixJε∗ has a characteristic polynomial given by(22)∑n=06Λnλn=0,where(23)Λ6=1,Λ5=ψ6+ψ5+γEhSEIdSh+PdSd+θEh+Ih+PhSh+Sdψ1IdEd+βDDEdId,Λ4=ψ6Sd+PhShSdPd+ψ6Sh+ψ5ShPh+ψ5+γSEIdEhSh+θEh+IhSh+βDDEdId+Sdψ1IdEdPdSd+ψ6+θIhSh+βDDEdId+Sdψ1IdEdψ5+γγ1SEIdSEh+EIdψ12Sd+γSEIdEhSh+PhSh+θEh+IhShSdψ1IdEd+βDDEdId,Λ3=γSEhγ1+EIdψ12SdEdβDDEd+ψ5PhSh+βDDEdId+Sdψ1IdEd+PdSd+γSEIdEhSh+θEhShβDDEdId+Sdψ1IdEd+PdSdψ6+γγ1SEIdSEh+EIdψ12Sd+PhSh+θIhShβDDEdId+Sdψ1IdEd+PdSd+βDDEdPdSdId+ψ1IdPdEdψ6+ψ5+Sdψ1SEhId+SEhEIdPdSd+Ed+SdβDDEdId+Id+SdSdψ1IdEdSEhγγ1+Sdψ12γSEIdEhSh+PhSh+θEh+IhSh+ψ12SdβDDEdIdEId+βDDEdPdSdId+ψ1IdPdEdγSEIdEhSh+PhSh+θEh+IhSh+γEIdSEIdSEhγ1ψ2.Clearly,Λn>0,∀n=3,4,5,6, and Λ0, Λ1, and Λ2 can similarly be shown to be positive. By Descartes’ rule of signs, the following result is thus established.Theorem 4.
The endemic equilibriumε∗ (whenever it exists) is locally asymptotically stable.


## 3.6. Global Stability of Disease-Free Equilibrium Point

To study the global stability ofℰ0, we employ the technique of Castillo-Chavez et al. [23] as follows:IfXG=Sd,Sh and ZG=Ed,Id,Eh,Ih, then model (1) can be written as(24)dXGdt=FXG,ZG,dZGdt=GXG,ZG,GXG,0=0.By the theorem of Castillo-Chavez et al. [23], ε0 is globally asymptotically stable if the following conditions are satisfied:GS1:ε0 is locally asymptotically stableGS2:XG∗ is globally stable for dXG/dt=GX,0GS3:GXG,ZG=AGZG−G^XG,ZG, G^XG,ZG≥0, where AG=DZGXG∗,0 is an M-MatrixThe first condition GS1 has been established forℛ0<1 in Theorem 3.For the second condition, we note thatFXG,0=Pd−μdSd−γ1ShSdPh+γShSd−μhSh is a limiting function of dXG/dt=FXG,ZG. Therefore, the second condition is satisfied.As for the third condition, we observe thatAG is given by(25)AG=Sd0ψ1−ψ3−γ1Sh0Sd0ψ100βDD−ψ4+γ1Sh000Sh0ψ2Sh0ψ2−ψ5000βH−ψ6,which is an M−matrix (i.e., all nondiagonal entries are nonnegative).Furthermore, the matrixG^XG,ZG is given by(26)G^XG,ZG=ψ1Sd0−SdEd+Id+Edγ1Sh+Eh−Sh0Idγ1Sh+Eh−Sh0ψ2Sh0−ShEd+Id0.FromdSd/dt, we have(27)Sdt=ϑt0Sdt0+Pdtϑt,where(28)ϑt=exp∫μd+γ1Sh+Eh+Ed+Idψ1dt.Hence,Sd0≥Sdt. Similarly, Sht≥Sh0 and hence, G^≥0, satisfies condition GS3. Now, since all three conditions are satisfied, the following result is established.Theorem 5.
The disease-free equilibrium point,ε0, is globally asymptotically stable whenever ℛ0≤1.For simulation purposes, we use the following parameter values:(29)Pd=0.45,Ph=0.1,βd=0.062,δd=0.5βDD=0.6,σd=0.8,βh=0.0229,δh=0.005,βH=0.0001,σh=0.0001,γ1=0.01,γ=0.0001,μd=0.06,Cd=0.05,θ=0.9,μh=0.04.


## 3.7. Sensitivity Analysis

To determine the impact of model parameters on the spread of rabies, we employed global sensitivity analysis technique of Marino et al. [24] to calculate the partial rank correlation coefficients of the model parameters on which the basic reproduction number depends. The result of sensitivity analysis informs which model parameters are most important in determining the output variable (in this case, ℛ0) and hence should be given much attention in terms of measurement for accuracy. The PRCC of results are shown in Figure 1.Figure 1
Plot of PRCC ofℛ0 with respect to model parameters.Local sensitivity analysis discusses the impact of small changes in parameter values on the spread of the disease. To study the sensitivity of the disease spread to parameter changes, the basic reproduction numberℛ0 is used, since it determines the persistence or possible eradication of the disease. The normalized forward sensitivity index is used to perform the analysis. It is defined as follows.Letℛ0 be a differentiable function xi. Then, the normalized forward sensitivity index of R0 relative to xi is given by(30)ηℛ0xi=∂ℛ0∂xi×xiℛ0.Local sensitivity indexes are presented in Table2. This index measures the relative percentage change in ℛ0 due to a percentage change in xi. The sensitivity indexes of the parameters determining ℛ0 are calculated and presented in Table 2. A positive sensitivity index of a parameter implies that an increase (decrease) in the value of the parameter will result in an increase (decrease) in ℛ0. On the other hand, a negative sensitivity index implies that an increase (decrease) in the value of the parameter will result in a decrease (an increase) in ℛ0.Table 2
Sensitivity indexes ofℛ0 with respect to model parameters.
C
d
P
d
P
h
β
D
D
β
d
δ
d
μ
d
μ
h
σ
d
γ
γ
1
−0.009
−2.010
−0.0143
−0.012
1.000
−1.000
1.056
3.024
−0.011
−3.967
0.943


## 3.8. Bifurcation Analysis

To study the stability of endemic equilibrium points, the center manifold theory as described in Theorem 4.1 of [25] is used as an alternative to the Lyapunov indirect method.Choosingβd=βd∗=ψ3+γ1Sh0ψ4+γ1Sh0/1−δdSd0γ1Sh0+σd+ψ3 as a bifurcation parameter leads the Jacobian model to have a zero simple eigenvalue. The right and left eigenvectors associated with this zero eigenvalue are given by w=w2η1,1,η3,η4,η5,η6T and v=0,v2,ζ3v2,0,0,0 with η3ζ3+1w2v2=1, where(31)η1=Sd0φ1γ1−φ2γ+μhφ2Sd0μd−μhμd+γ1Sh0,η3=βDDψ4+γ1Sh0,η4=φ2γSh0−φ1μd+γ1Sh0Sd0μd−μhμd+γ1Sh0,η5=ψ2Sh0γ1Sh0+ψ4+βDDψ4+γ1Sh0ψ5,η6=βHψ2Sh0γ1Sh0+ψ4+βDDψ4+γ1Sh0ψ5ψ6,ζ3=γ1Sh0+ψ3γ1Sh0+ψ4+βDD,φ1=Sd0ψ11+η3+γ1η5,φ2=Sh01+η3γ−ψ2+γSd0+θη5+θη6.The bifurcation coefficients are thus given by(32)a=2w2p∗−1,b=Sd01−δd1+η31+η3ζ3,where p∗=βd∗η11−δd1+η3/γ1η4+η51+η3ζ3.Clearly,b is always positive and the sign of a is dependent on p∗. The following result therefore easily follows [25].Theorem 6.
Model (1) undergoes a backward bifurcation at ℛ0=1 when p∗>1 and a forward bifurcation when p∗<1 in which case the endemic equilibrium will be locally asymptotically stable.We note that sinceℰ0 is GAS whenever ℛ0≤1, then backward bifurcation may only occur whenever ℛ approaches unity from above.


## 4. Numerical Simulation

We performed numerical simulation in this section to study the impact of various parameters on the spread of the rabies infection. To perform these experiments, the following initial values are used:Sd0=2.0, Ed0=1.5, Id0=1.0, Sh0=2.5, Ed0=0.5, and Ih0=0.5 and the baseline parameter values are taken from equation (29).To illustrate the local stability of the disease-free equilibrium, model (1) is solved for various initial values of the state variables and the results are plotted on the same graph. This procedure is done for the case when ℛ0<1 and when ℛ0>1. The fact that the DFE is LAS whenever ℛ0<1 as in Theorem 3 is demonstrated in Figure 2.Figure 2
Simulation of model (1) showing the local asymptotic stability of the ℰ0 (a) for ℛ0<1 and instability (b) for ℛ0>1.
(a)
(b)To illustrate the result in Theorem1, we solved model (1) for two different scenarios: (a) when initial conditions are less than Ph+Pd/minμh,μd and (b) when initial conditions are greater than Ph+Pd/minμh,μd. Figure 3 shows that the feasible region Ω is an invariant attractor of system (1).Figure 3
Simulation of the model illustrating the boundedness of solutions of model (1) (Theorem 1).We also solved the model for different values of the most influential model parameters in order to simulate the effect of those variables. The results are shown in Figures4–8.Figure 4
Solution curves of model (1) showing the effect of βd.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Solution curves of model (1) showing the effect of δd.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Solution curves of model (1) showing the effect of σd.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Solution curves of model (1) showing the effect of βDD.
(a)
(b)
(c)
(d)
(e)
(f)Figure 8
Solution curves of model (1) showing the effect of Pd.
(a)
(b)
(c)
(d)
(e)
(f)From Figure4, it is observed that increasing the dog-dog rabies transmission rate, βd leads to a decline (an increase) in the Susceptible dog population (Exposed and Infected dog populations) and a similar (but slight) effect is observed for the Susceptible human population (Exposed and Infected human populations). The smaller impact on human population could be attributable to the fact that infected dogs will often be culled, leading to reduced infection among humans.From Figure5, it is observed that increasing the rate of immunization of dogs, δd, leads to an increase (a decrease) in the Susceptible dog population (Exposed and Infected dog populations) and a similar effect is observed for the Susceptible human population (Exposed and Infected human populations). It is however observed that for higher rates of immunization, the Susceptible dog populations first peak before declining. This could be attributable to the fact that immunization may not confer permanent reduction in transmissibility of rabies from dogs. We also observe that irrespective of the rates of immunization, a steady state is reached, confirming the global stability of the disease-free equilibrium.


## 5. Findings

The following are the findings of the study:(i)
The rabies-free equilibrium point is globally asymptotically stable whenever the basic reproduction numberℛ0 is less than unity.(ii)
While the spread of disease is positively correlated with the dog-to-dog transmission rateβd, recruitment rate of dogs Pd, rate at which exposed dogs develop clinical rabiesβDD, and disease-induced death rate among σd, and natural death rate of humans μh, it is negatively correlated with the natural death rate of dogs μd, immunization of dogs δd, culling rate Cd, rate of consumption of dogs γ1, as indicated in the sensitivity indexes of ℛ0.(iii)
Efforts at controlling the spread of rabies should focus more on dog population than on humans.(iv)
Under certain conditions, the model may exhibit either forward or backward bifurcation nearℛ0=1. The possibility of forward bifurcation suggests the existence of locally asymptotically stable endemic equilibrium when ℛ0>1.


## 6. Discussion and Conclusion

In this paper, a nonlinear mathematical model has been proposed to study the dynamics of rabies infection among dogs and in a human population that is exposed to dog bites. It is observed that the model has a unique disease-free equilibrium point which is locally asymptotically stable whenever the basic reproduction number is less than unity. The disease-free equilibrium point,ℰ0, is also globally asymptotically stable whenever ℛ0≤1. It is also observed that the most influential factors on the spread of the disease are probability of infection upon contact among dogs βd, rate of immunization of dogs δd, rate of recruitment of dogs Pd, rate at which exposed dogs develop clinical rabies βDD, disease-induced death rate within the dog population σd, and natural mortality rate of dogs, μd. If the rate of immunization in the dog population is increased, it will have an increasing effect on the susceptible human population as shown in Figure 5. If the recruitment rate within the dog population, disease-induced death rate of dogs, rate of contact among dogs, and the rate at which exposed dogs develop clinical rabies each is reduced, with all other factors remaining constant, there will be a positive effect on the susceptible human population as shown in Figures 4 and 6–8. These actions will lead to a reduction in the reproduction number. However, even though it has been accepted that the eradication of infectious diseases in a community should target the reduction of ℛ0 to less than unity, this may not be the only case for this system. It has also been observed that under certain conditions, backward bifurcation may occur around ℛ0=1. From the analysis, we noted that since ε0 is GAS whenever ℛ0≤1, then backward bifurcation may only occur whenever ℛ0 approaches unity from above. This has a major implication for the spread of the disease once the basis of reducing ℛ0 below the threshold point, ℛ0=1, may not just be enough for effectively controlling the spread of the disease because whenever ℛ0 approaches unity from above, there is the likelihood of coexistence of endemic equilibrium which may result in reinfection of the rabies virus. The conclusions, however, of Asamoah et al. [5]; Hou et al. [12]; and Wang and Lou [19] are that vaccination, which may lead to immunity, as a better control method of rabies infection in dogs will still be relevant in the presence of backward bifurcation.



---

*Source: 1804270-2020-03-30.xml*