# A Mathematical Model Analysis for the Transmission Dynamics of Leptospirosis Disease in Human and Rodent Populations

Authors: Habtamu Ayalew Engida; David Mwangi Theuri; Duncan Gathungu; John Gachohi; Haileyesus Tessema Alemneh
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

This work is aimed at formulating and analyzing a compartmental mathematical model to investigate the impact of rodent-born leptospirosis on the human population by considering a load of pathogenic agents of the disease in an environment and the incidence rate of human infection due to the interaction between infected rodents and the environment. Firstly, the basic properties of the model, the equilibria points, and their stability analysis are studied. We also found the basic reproduction numberR0 of the model using the next-generation matrix approach. From the stability analysis, we obtained that the disease-free equilibrium (DFE) is globally asymptotically stable if R0<1 and unstable otherwise. The local stability of endemic equilibrium is performed using the phenomenon of the center manifold theory, and the model exhibits forward bifurcation. The most sensitive parameters on the model outcome are also identified using the normalized forward sensitivity index. Finally, numerical simulations of the model are performed to show the stability behavior of endemic equilibrium and the varying effect of the human transmission rates, human recovery rate, and the mortality rate rodents on the model dynamics. The model is simulated using the forward fourth-order Runge-Kutta method, and the results are presented graphically. From graphical stability analysis, we observed that all trajectories of the model solutions evolve towards the unique endemic equilibrium over time when R0>1. Our numerical results revealed that decreasing the transmission rates and increasing the rate of recovery and reduction of the rodent population using appropriate intervention mechanisms have a significant role in reducing the spread of disease infection in the population.


---

## Body

## 1. Introduction

Leptospirosis is an infectious bacterial disease caused by a pathogenic spirochete bacteria calledLeptospira interrogans, and it occurs throughout the world but is most common in tropical and subtropical regions, especially in developing countries including South-East Asia countries and sub-Sahara Africa [1–3]. It is one of the major direct zoonosis diseases worldwide that affects humans and animals [4, 5]. Recent studies reported that more than one million human cases with an estimated 60,000 deaths occur worldwide each year due to the Leptospira interrogans [1, 4].Human infection occurs through direct contact with infected animal reservoir urine, tissues, or other body fluid, or more commonly by contact with infected animal (reservoirs) urine-contaminated environment [1, 6, 7]. Person-to-person transmission of the disease occurs very rarely. Globally, the main animal reservoirs of human leptospirosis are rodents (rats and mice), especially in urban slum environments [8–10]. The high concentration of Leptospira bacteria shed in slum environments mainly occurs due to the high prevalence of infection in rat population [6].The incubation period of leptospirosis disease is typically 5 to 14 days [1, 11]. The symptoms and signs of the disease are nonfixed and may be confused with other diseases’ symptoms (dengue, hantavirus, malaria, melioidosis, influenza, etc.) due to the initially nonspecific presentation. Consequently, its cases are sometimes underrecognized [11]. In humans, leptospirosis illness can be demonstrated by two stages. The primary (an acute) stage is illustrated by mild illness with nonspecific signs like higher fever, headache, and conjunctival suffusion [1, 12], while syndromes in the second (severe illness or immune) stage of leptospirosis include jaundice, kidney failure, haemorrhage (especially pulmonary), meningitis, cardiac arrhythmias, respiratory insufficiency, and hemodynamic collapse [1, 11].Human vaccine against leptospirosis is not widely practiced and accessible only in a few developed countries which protect only against the serovar in the vaccine, and regular boosting is required [11]. However, prevention and control interventions can be implemented to avoid the risk of leptospirosis infections (especially for people whose working environment exposes them to the risk of infections), to minimize the acquisition of the infection, or to reduce the human infections in infected population. Preventive measures include personal protective equipment (PPE) (wearing rubber boots, waterproof overalls/dressings to cover wounds or skin, goggles, and rubber gloves) and personal hygiene. Environmental modifications such as draining wet areas or improvements in urban slum environment to reduce the concentration of pathogenic spirochete leptospiries from the environment [1, 12].Patients with early leptospirosis infection (mild illness) can be treated through antibiotic doxycycline, ampicillin or amoxicillin, azithromycin, or clarithromycin. Patients with severe illness can be treated through IV penicillin and ceftriaxone drugs [1, 11].Recently, a mathematical model has become an important tool in underlying mechanisms of disease transmission and spread and has been used to predict outcomes of diseases in communities, help to explain key factors in the disease transmission process, suggest effective control and preventive measures, and provide an estimate for the seriousness and potential scale of the epidemic [13]. In particular, deterministic mathematical models have been used to study the dynamics of infectious diseases within human and vector hosts in the population by dividing the individuals into different stages and making assumptions about the nature and time rate of transfer from one compartment to another [14].A few researchers have been investigated a two-strain deterministic mathematical model, one strain for the human population and another strain for vectors (or animal population), for the transmission process of leptospirosis disease in various forms to demonstrate the dynamical spread of the disease between the two population groups.Pimpunchat et al. [15] proposed the susceptible, infected, removed (SIR) model of leptospirosis which consists of both human and vector host populations. In their model, the human population is divided into three compartmental, susceptible, infected, and recovered humans, whereas the vector population was divided into two classes, susceptible and infected vectors. They performed numerical solutions to analyze and examine the behavior of the disease. Likewise, another compartmental mathematical model was developed by incorporating exposed cases in both human and vector populations for transmission dynamics of leptospirosis disease in [16]. They illustrated numerical simulations of the formulated model by using the well-known numerical method, the Runge-Kutta of fourth-order in the MATLAB program. Authors in [17] presented a compartmental mathematical model for the dynamic behavior of the leptospirosis disease with saturated incidence. They described the model using nonlinear ordinary differential equations by subdividing the host population of humans into three classes, susceptible individuals, infected individuals, and recovered individuals, and also by subdividing vectors (animal) population into two compartments, susceptible and infected vectors. They investigated the stability analysis of the disease-free and endemic equilibrium by using the basic reproduction number. From the results of stability analysis, they obtained that the disease-free equilibrium is stable both locally and globally whenever the basic reproduction number is less than one. Similarly, authors in [18] investigated a mathematical model of the leptospirosis disease using a compartmental approach. However, in view of above models, the contraction of Leptospira interrogans in environment for transmission had not been incorporated.Although the human infection of the disease occurs in multiple different ways, most of the infections acquired directly or indirectly from environment (soil or water) contaminated with Leptospira especially in urban slum areas [7–9].Motivated from [5, 17], this study examines a two-strain model, SEIRS model for the human population and another SIRS model for rodent (vector) population taking account the concentration leptospires population in environment and the transmission rates between susceptible individuals and the load of the pathogen in the environment and between susceptible rodents and infected rodents for the dynamics behavior of leptospirosis disease. Thus, in this work, we present two host populations (namely, humans and rodents) and the pathogenic population with its two growth controls (namely, the natural death rate of Leptospira and carrying capacity controls). The rest of this work is organized as follows. The leptospirosis model is described and formulated in Section 2. Section3 is concerned with the basic properties of the model, local and global stability of disease-free equilibrium, the phenomena of bifurcation analysis, and local stability of endemic equilibrium. Section 4 is devoted to the sensitivity analysis for the basic reproduction number R0 of the model to each of the parameter values. The numerical simulations and discussions are presented in Section 5. Our conclusions are given in Section 6.


## 2. Model Formulation

Using system of nonlinear differential equations, we build up a compartmental mathematical model for the transmission process of leptospirosis epidemic. The model includes human, vector (rodent), and bacterial populations. The total human population denoted byNht is subdivided into four compartments: Sht, Eht, Iht, and Rht. Thus,
(1)Nht=Sht+Eht+Iht+Rht.The total rodent population denoted byNvt is subgrouped into three groups: Svt, Ivt, and Rvt with
(2)Nvt=Svt+Ivt+Rvt.Moreover, the state variables of the model are described in Table1.Table 1
Description of variables of the leptospirosis disease model.
VariableDescriptionShtThe susceptible individuals; individuals who are healthy but can be infectedEhtLatently infected individuals; individuals having the disease but not showing the symptomsIhtThe infectious individuals that are showing symptoms of the disease and can transmit disease to othersRhtThose that have recovered from the disease and have got temporary immunitySvtSusceptible rodentsIvtInfected rodentsRvtRodents that are recovered from the diseaseBlLeptospirosis causing bacterial populationIn the formulation of the model, the following assumptions are considered:(i)
Each parameter of the model cannot be negative(ii)
We assume that susceptible human and rodent populations increase at constant rates are given by recruitment of the individualsΛ and Π, respectively(iii)
We assume that susceptible humans can be infected in two different ways: through either direct contact with infected rodents urine or contact with contaminated environments(iv)
Assume there is a homogeneous mixing between human and rodent populationsThe susceptible individuals can acquire infection with the incidence rateλh=λhB+λhv with λhB=β2Bl/κ+Bl and λhv=β1Iv. The susceptible rodents can acquire infection with the force of infection β3Ih. The population Bl increases its size in contaminated environment from the release of bacteria by infected humans and rodents with the rate of τ1 and τ2, respectively. Further, description of the parameters of the model is summarized in Table 2.Table 2
Description of parameters of the leptospirosis disease model.
ParameterDescriptionΛThe recruitment rate of susceptible humansΠThe recruitment rate of susceptible rodentsβ1The transmission coefficient of the disease from infected rodents to susceptible humansβ2The transmission rate of infection from environment to humanβ3The transmission rate of infection from infectious individuals to susceptible rodentsγThe disease waning immunity for humansμThe natural death rate of the human populationμvThe natural death rate of rodent populationμbThe natural death rate of leptospirosis causing bacteria populationθThe rate of the exposed human move to infected classαThe death rate due to the disease infectionδThe recovery rate from symptomatic infectiousρThe disease waning immunity for rodentsσRate of recovery from leptospirosis rodent infectionκThe concentration of pathogenic population in environmentτ1The rate of increase of bacteria byIhτ2The rate of increase of bacteria byIvFrom the flowchart diagram in Figure1, we obtained the following set of differential equations for the compartmental model of the leptospirosis disease.
(3)dShdt=Λ+γRh−λh+μSh,dEhdt=λhSh−θ+μEh,dIhdt=θEh−α+δ+μIh,dRhdt=δIh−γ+μRh,dSvdt=Π+ρRv−β3Ih+μvSv,dIvdt=β3IhSv−σ+μvIv,dRvdt=σIv−ρ+μvRv,dBldt=τ1Ih+τ2Iv−μbBl,Figure 1
Compartmental flow diagram of leptospirosis disease transmission.whereλh=β2Bl/κ+Bl+β1Iv. With the initial conditions:
(4)Sh0≥0,Eh0≥0,Ih0≥0,Rh0≥0,Sv0≥0,Iv0≥0,Rv0≥0,Bl0≥0.


## 3. Model Analysis

### 3.1. Positivity of Solutions

Theorem 1.
LetΩl0=Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Rv0,Bl0∈R+8∪0 be the initial condition for equation (3). Then, the set of solutions Sht,Eht,Iht,Rht,Svt,Ivt,Rvt,Blt of equation (3) is nonnegative for all t≥0.Proof.
Given that a set of nonnegative initial conditionsSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Rv0,Bl0, consider equation (3). Let t1=supt>0:Sht0>0,Eht0>0,Iht0>0,Rht0>0,Svt0>0,Ivt0>0,Rvt0>0,Blt0>0,∀t0in0,t. From the first equation of the system (3), it follows that
(5)Shdt+λht+μSh=Λ+γRht,which can be rewritten as
(6)ddtShteμt+∫0tλhτdτ=Λ+γRhteμt+∫0tλhτdτ.
Thus,(7)Sht1eμt1+∫0t1λhτdτ−Sh0=∫0t1Λ+γRhteμz+∫0zλhτdτdz.
Therefore,(8)Sht1=Sh0e−μt1−∫0t1λhτdτ+e−μt1−∫0t1λhτdτ×∫0t1eμz+∫0zλhτdτΛ+γRhzdz≥0.
Similarly, the fifth equation of the system (3) yields that
(9)Svt1=Sl0e−μlt1−∫0t1λlτdτ+e−μlt1−∫0t1λlτdτ×∫0t1eμlz+∫0zλlτdτΠ+ρRvzdz≥0.
Next, from the second equation of model (3), it follows that
(10)dEhdt≥−θ+μEh.
Applying the comparison theorem and variation of formula for (10) yields
(11)Eht≥Eh0e−θ+μt≥0,∀t≥0.
Similarly,(12)Iht≥Ih0e−α+δ+μt≥0,∀t≥0,Rht≥Rh0e−γ+μt≥0,∀t≥0,Ivt≥Iv0e−σ+μvt≥0,∀t≥0,Rvt≥Rv0e−ρ+μvt≥0,∀t≥0,Blt≥Bl0e−μbt≥0,∀t≥0.
Therefore, all solutions of the system (3) remain nonnegative for all nonnegative initial conditions.


### 3.2. Invariant Region

Theorem 2.
All solutions of the leptospirosis model equation (3) remain in Ωl=Sh,Eh,Ih,Rh,Sv,Iv,Rv,Bl∈R+8:0≤Nht≤Λ/μ,0≤Nvt≤Π/μv,0≤Blt≤τ∗/μblΛ/μ+Π/μv.Proof.
From equation (1), we get
(13)dNhdt=dShdt+dEhdt+dIhdt+dRhdt.
The first four equations of the system (3) and equation (13) give rise to
(14)dNhdt=Λ−μNh−αIh.
Since the parameter,α, and the state variable Ih are nonnegative in (14), we have
(15)dNhdt≤Λ−μNh.
LetNh∗t be the solution of the ODE dNh/dt=Λ−μNh with an initial condition. Then,
(16)dNh∗dt=Λ−μNh∗withNh∗0=Nh0∗∀t≥0.
Thus, equation (16) has a unique solution. By applying separation of variable on the differential equation (16), we get
(17)dNh∗Λ−μNh∗=dt.
Integrating (17) on both sides yields
(18)Nh∗t=Nh∗0e−μt+Λμ1−e−μt.
We now apply the comparison theorem [19] on the differential equation (15) and then Nht≤Nh∗t. Thus,
(19)Nht≤Nh∗0e−μt+Λμ1−e−μt.
Ast⟶∞, the population size Nh⟶Λ/μ.
Hence,(20)Nht≤Λμfor0≤Nh∗0≤Λμ,∀t≥0.
On the other hand, ifNh∗0>Λ/μ, the solution decrease to Λ/μ as t⟶∞.
Also, the last equation of (3) and equations (20) and (23) imply that
(21)dBldt≤τ∗Λμ+Πμv−μbBl,
whereτ∗=maxτ1,τ2.
Equation (21) yields
(22)Blt≤Bl0e−μblt+τ∗μblΛμ+Πμv1−e−μbt∀t≥0.
Hence,(23)Ωh=Sh,Eh,Ih,Rh,Bl∈R+5:0≤Nht≤Λμ,0≤Bl≤τ∗μblΛμ+Πμv,is positively invariant region.
Finally, the 5th equation of (3) implies that
(24)dNvdt≤Π−μNv,withNv0=Nv,0.
Then, equation (24) has a unique solution. Similarly, applying the comparison theorem and integrating the differential equation (21) yield
(25)Nvt=Nv0e−μvt+πμv1−e−μvt.
Therefore,(26)0≤Nvt≤Πμvfor0≤Nv0≤Πμv,Ωv=Sv,Iv,Rv∈R+3:0≤Nvt≤Πμv.
Therefore, the biologically feasible region for (3) is given by
(27)Ωl=Ωh×Ωv⊂R+5×R+3.
Therefore, every solution of the differential equation model with initial conditions inΩl remains in Ωl∀t≥0. Thus, the region Ωl is positively invariant with respect to the system (3).


### 3.3. Disease-Free Equilibrium

The disease-free equilibrium (DFE) of the model (3) only exists in the absence of infections. We compute DFE by setting left-hand side of the system (3) equal to zero. Since Eh∗=Ih∗=Rh∗=Iv∗=Rv∗=Bl∗=0 at DFE, we have
(28)Λ+γRh∗−β1Bl∗κ+Bl∗+β2Iv∗+μSh∗=0,Π+ρRv∗−β3Il∗+μvSv∗=0.Thus, the system (28) yields Sh∗=Λ/μ and Sv∗=Π/μv. Therefore, the DFE point of the model (3) is given by
(29)ε0l∗=Sh∗,0,0,0,Sv∗,0,0,0=Λμ,0,0,0,Πμv,0,0,0.To compute the basic reproduction number for our model (3), we follow the method presented in [20]. Considering only the infective compartments W=Eh,Ih,Iv,Bl in (3), then we have dW/dt=Ft−Vt,
(30)Ft=λhSh0β3SvIh0,Vt=ε1Eh−θEh+ε2Ihε4Iv−τ1Ih−τ2Iv+μblBl,where(31)λh=β2Blκ+Bl+β1Iv,ε1=θ+μ,ε2=α+δ+μ,The Jacobian matrices ofFt and of Vt at ε0l∗ are, respectively, given by
(32)F=00β1Λμβ2Πμκ00000β3Πμv000000,V=ε1000−θε20000ε300−τ1−τ2μb.Therefore, the next-generation matrix of the model is given by(33)3G=FV−1=00β1Λμβ2Λμκ00000β3Πμv0000001ε1000θε1ε21ε200001ε40θτ1ε1ε2μbτ1ε2μbτ2ε4μb1μb=y1y2y3y40000y5y6000000,where
(34)y1=β2Λθτ1κμμbε1ε2,y2=β2Λτ1κμμbε2,y3=Λβ2τ2+β1κμbκμμbε4,y4=β2Λκμμb,y5=β3Πθμvε1ε2,y6=β3Πμvε2.ConsideringG−λI=0, we have
(35)y1−λy2y3y40−λ00y5y6−λ0000−λ=0⟹λ2λ2−λy1−y3y5=0.Thus,(36)λ1=λ2=0,λ3=12y1−y12+4y3y5<0,λ4=12y1+y12+4y3y5>0,are the eigenvalues of the matrixG. Therefore, the basic reproduction number of the model is given by
(37)R0=ρG=12β2Λθτ1κμμbε1ε2+β2Λθτ1κμμbε1ε22+4β3ΠΛθβ2τ2+β1κμbκμμvμbε1ε2ε4.The basic reproduction numberR0 can be written in terms of the basic reproduction numbers corresponding to the vector transmission Rv, human transmission Rh, and environment-to-human transmission Rb as
(38)R0=12Rb+Rb2+4Rhv,whereRb=β2Λθτ1/κμμbε1ε2 and Rhv=Rh+Rv with Rh=β3β2ΠΛθτ2/κμμvμbε1ε2ε4 and Rv=β3β1ΠΛθ/μμvε1ε2ε4.


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

Theorem 3.
The disease-free equilibrium point,ε0l∗, is locally asymptotically stable whenever R0<1 and unstable otherwise.Proof.
The Jacobian matrices of the model (3) at DFE ε0l∗ is given by
(39)Jε0l∗=−μ00γ0−β1Λμ0−β2Λμκ0−ε1000β1Λμ0β2Λμκ0θ−ε20000000δ−ε3000000−β3Πμv0−μv0ρ000β3Πμv00−ε40000000σ−ε5000τ100τ20−μb,
whereε1=θ+μ,ε2=α+δ+μ,ε3=γ+μ,ε4=σ+μv, and ε5=ρ+μv.
Expanding the determinant of the characteristic equationJε0l∗−λI=0 by the first, fourth, fifth, and seventh columns, respectively, we obtain four of the eigenvalues of Jε0l∗: λ1=−μ,λ2=−γ+μ,λ3=−μvandλ4=−ρ+μv. The other four eigenvalues will be determined from the following polynomial equation:
(40)−ε1−λ0β1Λμβ2Λμκθ−ε2−λ000β3Πμv−ε4−λ00τ1τ2−μb−λ=0,⟹P4λ=λ4+A1λ3+A2λ2+A3λ+A4=0,
where(41)A1=ε1+ε2+ε4+μb,A2=ε1μb+ε1ε2+ε1ε4+ε2μb+ε4μb+ε2ε4,A3=ε1ε2μb+ε1ε4μb+ε1ε2ε4+ε2ε4μb−β2Λτ1θμκ−β1β3ΛΠθμμv=ε4μbε1+ε2+ε1ε2μb1−β2Λθτ1κμμbε1ε2+ε1ε2ε41−β1β3ΛΠθμμvε1ε2ε4=ε4μbε1+ε2+ε1ε2μb1−Rb+ε1ε2ε41−Rv,A4=ε1ε2ε4μb−β2Λθτ1ε4μκ+β1β3ΛΠθμbμμv+β2β3ΛΠθτ2κμμv=ε1ε2ε4μb−β2Λθτ1ε4μκ+β3ΠΛθβ2τ2+β1κμbκμμv=ε1ε2ε4μb1−β2Λθτ1ε4μκ+β3ΠΛθβ2τ2+β1κμbκμμv=ε1ε2ε4μb1−Rb+Rhv.
Since all the model parameters are positive, the above four eigenvalues of the matrixJε0l∗ are negative real numbers. To ensure that all the remaining four eigenvalues of the matrix have negative real parts, the Routh-Hurwitz stability criterion (see Theorem 5.1 in [21]) requires A1>0,A2>0,A3>0,A4>0,A1A2−A3>0, and A1A2A3−A32+A12A4>0. Among these, A1>0 and A2>0 are obvious.
Similarly,(42)A1A2−A3=A2ε4+μb+ε22ε1+ε4+μb+ε1ε1μb+ε1ε2+ε2μb+ε1ε4+ε2ε4+β2Λτ1θμκ+β1β3ΛΠθμμv>0.
The other three inequalities hold whenR0>1; the details are described as follows.
SinceRb<R0, 1−Rb>0 for R0<1. Similarly, 1−Rv=1−R01−Rb+R0+Rb+Rh>0 for R0<1. Thus,
(43)A3=ε4μbε1+ε2+ε1ε2μb1−Rb+ε1ε2ε41−Rv>0ifR0<1.
Next,(44)A4=ε1ε2ε4μb1−Rb+Rhv=ε1ε2ε4μb1−R01−Rb+R0>0ifR0<1.
Finally, we need to show thatA1A2A3−A32+A12A4>0. Let
(45)ϕ0=ε1μb,ϕ1=ε1ε2,ϕ2=ε1ε4,ϕ3=ε2μb,ϕ4=ε4μb,ϕ5=ε2ε4,ψ0=ε4μbε1+ε2,ψ1=ε1ε2μb1−Rb,ψ3=ε1ε2ε41−Rv.
Then,A2=ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5, and A3=ψ0+ψ1+ψ2.
Now,(46)A1A2A3=ε1ψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ε2ψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ε4ψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+μbψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5,A32+A12A4=ψ0+ψ1+ψ22+ε1+ε2+ε4+μb2ε1ε2ε4μb1−Rb−Rhv=ψ02+ψ12+ψ22+2ψ0ψ1+ψ0ψ2+ψ1ψ2+ε12+ε22+ε42+μb2+2ε1ε2+ε1ε4+ε1μb+ε2ε4+ε2μb+ε4μbε1ε2ε4μb1−Rb−Rhv.
After manipulating the parameters algebraically, we obtain(47)A1A2A3−A32+A12A4=Ω0+Ω1,
where(48)Ω0=ε1ψ0ϕ0+ϕ1+ϕ2+ϕ4+ψ1ϕ0+ϕ1+ϕ4+ψ2ϕ1+ϕ2+ϕ4+ε2ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ψ2ϕ0+ϕ1+ϕ2+ϕ4+ϕ5+ε4ψ0ϕ1+ϕ2+ϕ4+ψ2ϕ1+ϕ2+ϕ4+ϕ5+μbψ0ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5>0,Ω1=ε1ϕ3RbA3+ψ0+ϕ5A3Rv+ε12+ε22+A2ϕ1ϕ41+Rh+ϕ4ϕ52+ϕ02+ε2A2+ε42+μb2+A2ϕ1ϕ4Rb+Rhv>0.
Clearly,Ω0>0 and Ω1>0.
Therefore,A1A2A3−A32+A12A4>0 if R0<1. This completes the proof.


### 3.5. Global Stability of the Disease-Free Equilibrium

The method illustrated in [22, 23] is used to investigate the global asymptotic stability (GAS) of DFE point of the model. Firstly, the model (3) must be written in the form:
(49)dXdt=FX,Y,dYdt=GX,Y,GX,0=0.whereX=Sh,Rh,Sv,Rv represents the number of uninfected individuals and Y=Eh,Ih,Iv,Bl denotes the number of infected individuals. Let X∗ be the disease-free equilibrium of the system dX/dt=FX,0. Then, X∗=Λ/μ,0,Π/μv,0. The DFE point of the model ε0l∗=X∗,0=Λ/μ,0,0,0,Π/μv,0,0,0 is guaranteed to be GAS if R0<1 (which is locally asymptotically stable (LAS)), and the following two conditions H1 and H2 hold:H1: For dX/dt=FX,0, X∗ is globally asymptotically stable.H2: GX,Y=AY−G⋆X,Y,G⋆X,Y≥0,∀X,Y∈Ωl, where A=∂Gε0l∗/∂Y is an M-matrix and Ωl is the region where the model makes biological sense. If the model equation (3) satisfies the above two conditions, then the following theorem holds.Theorem 4.
The disease-free equilibrium point,ε0l∗=Λ/μ,0,0,0,Π/μv,0,0,0, is GAS for the model (3) provided that R0<1 (LAS) and conditions H1 and H2 hold.Proof.
We only need to show that the conditions H1 and H2 hold when R0<1. Considering the system (3), we have
(50)FX,0=Sh′tRh′tSv′tRv′t=Λ+γRh−μSh−γ+μRhΠ+ρRv−μvSv−ρ+μvRv,X∗=Λμ,0,Πμv,0.
Note that the second and the fourth equations of the system (50) are first-order linear ODEs, and their solutions are given by, respectively,
(51)Rht=Rh,0e−γ+μt,(52)Rvt=Rv,0e−ρ+μvt.
Comparing equation (51) and the first equation of the system (50) yields
(53)S′ht=Λ+γRv,0e−ρ+μv−μSht.
Solving the first-order linear ODE (53), its solution is found as
(54)Sht=Λμ+Sh,0e−μt−Rh,0e−γ+μt.
Similarly,(55)Sv′t=Π+ρRv,0e−ρ+μvt−μvSvt,and its solution is given by
(56)Svt=Πμv+Sv,0e−μvt−Rv,0e−ρ+μvt.
Now, suppose that the time,t⟶∞, we need to show that X⟶X∗. Obviously, Rht⟶0, Rvt⟶0 as t⟶∞. Thus, Sht⟶Λ/μ, Svt⟶Π/μv as t⟶∞, regardless of the values of Sh0, Rh0 and Sv0, Rv0, respectively. Thus, all points with respect to this condition converge at X0∗=Λ/μ,0,Π/μv,0. Hence, X∗=Λ/μ,0,Π/μv,0 is GAS.
Next, we have(57)GX,Y=G1X,YG2X,YG3X,YG4X,Y=β2Blκ+Bl+β1IvSh−θ+μEhθEh−α+δ+μIhβ3SvIh−σ+μvIvτ1Ih+τ2Iv−μbBl.
We can then obtain(58)A=∂GX∗,0∂Y=−θ+μ0β1Λμβ2Λκμθ−α+δ+μ000β3Πμv−σ+μv00τ1τ2−μb,which is clearly an M-matrix
(59)G⋆X,Y=AY−GX,Y=β1IvΛμ−Sh+β2BlκΛμ−κShκ+Bl0β3IhΠμv−Sv0.
Since all parameters are nonnegative,0≤κSh/κ+Bl≤Sh≤Λ/μ∵κ/κ+Bl≤1, 0≤Sv≤Π/μv. It follows that Λ/μ−Sh≥0, Λ/μ−κSh/κ+Bl≥0, and Π/μv−Sv≥0.
Hence,G⋆X,Y≥0∀X,Y∈Γ. Therefore, the disease-free equilibrium point, ε0l∗=Λ/μ,0,0,0,Π/μv,0,0,0, of the model (3) is globally asymptotically stable. This completes the proof.Remark 5.
The model (3) does not exhibit a backward bifurcation at R0=1 when R0<1, since DFE is the only stable (or positive) equilibrium point for R0<1 in the Theorem 4.The backward bifurcation has been observed in the leptospirosis epidemic model [16].


### 3.6. Endemic Equilibrium (EE)

The EE point of our model is computed by setting the system of differential equations in (3) to zero. Thus, it is obtained by dSh∗/dt=dEh∗/dt=dIh∗/dt=dRh∗/dt=dSv∗/dt=dIv∗/dt=dRv∗/dt=dBl∗/dt=0, which is given by El∗=Sh∗,Eh∗,Ih∗,Rh∗,Sv∗,Iv∗,Rv∗,Bl∗, where all disease states are considered to be positive, and thus, Ih∗ must be greater than zero for all the other states to be positive. Now, by setting the system of equations in (3) to zero at steady state, we have
(60)dSh∗dt=Λ+γRh∗−β2Bl∗κ+Bl∗+β1Iv∗+μSh∗=0,dEhdt=β2Bl∗κ+Bl∗+β1Iv∗Sh∗−θ+μEh∗=0,dIh∗dt=θEh∗−α+δ+μIh∗=0,dRh∗dt=δIh∗−γ+μRh∗=0,dSv∗dt=Π+ρRv∗−β3Ih∗+μvSv∗=0,dIv∗dt=β3Ih∗Sv∗−σ+μvIv∗=0,dRv∗dt=σIv∗−ρ+μvRv∗=0,dBl∗dt=τ1Ih∗+τ2Iv∗−μbBl∗=0.Due to the complexity of the model system (3), all other state variables at the steady state are expressed in terms of the steady state of Ih which is denoted as Ih∗. Accordingly, we found as
(61)El∗=Sh∗,Eh∗,Ih∗,Rh∗,Sv∗,Iv∗,Rv∗,Bl∗,=Λ+γRh∗λ∗+μ,ε2θIh∗,Ih∗,δε3Ih∗,Π+ρRv∗β3Ih∗+μv,β3Πε5Ih∗μvε4ε5+β3ε6Ih∗,β3ΠσIh∗μvε4ε5+β3ε6Ih∗,τ1Ih∗+τ2Iv∗μb,where(62)ε1=θ+μ,ε2=α+δ+μ,ε3=γ+μ,ε4=σ+μv,ε5=ρ+μv,ε6=σ+ρ+μv.By combining the equations of (61) and the first and second equations of the model system (3), we obtained a septic (or heptic) polynomial equation for Ih∗ which is described as follows:
(63)PIh∗=Ih∗AIh∗6+BIh∗5+CIh∗4+DIh∗3+EIh∗2+FIh∗+G=0,where(64)A=−Ψ3Ψ7η4+Ψ72θ+μα+δ+μγ+μμ<0always,B=−η4Ψ2Ψ7+Ψ3Ψ6+Ψ7η3θ+μα+δ+μγ+μμη32κμbRb−1+η3τ1σ+μvρ+μvμvRv−1−η3η2τ1+η1τ2−Ψ6<0ifR0<1,C=−η4Ψ1Ψ7+Ψ2Ψ6+Ψ3Ψ5+Ψ7Ω2−Ψ7Ψ5θ+μα+δ+μγ+μμ+Ψ6Ω3=−η4Ψ1Ψ7+Ψ2Ψ6+Ψ3Ψ5+Ψ6Ω3+γ+μθ+μα+δ+μμη2η3κμb2Rb−1+σ+μvρ+μvμvRv−2η3κμb+η2τ1+η1τ2<0ifR0<1,D=−η4Ψ1Ψ6+Ψ2Ψ5+Ψ3Ψ4+Ψ7Ω1−Ψ4θ+μα+δ+μγ+μμ+Ψ6Ω2+Ψ5Ω3,E=−η4Ψ1Ψ5+Ψ2Ψ4+Ψ6Ω1+Ψ5Ω2+Ψ4Ω3,F=−η4Ψ1Ψ4+Ψ5Ω1+Ψ4Ω2,G=Ψ4Ω1<0ifR0<1,with(65)η1=Πβ3ρ+μv,η2=μvσ+μvρ+μv,η3=β3μvσ+ρ+μv,η4=θγα+μ+μα+δ+μθ+γ+μ,Ψ1=η2β2η2τ1+η1τ2+β1κη1μb,Ψ2=η2τ1+η1τ2β2η3+β1η1+η3β2η2τ1+β1κη1μb,Ψ3=η3τ1β2η3+β1η1,Ψ4=κη22μb,Ψ5=2κη2η3μb+η2η2τ1+η1τ2,Ψ6=η3κη3μb+2η2τ1+η1τ2,Ψ7=η32τ1,Ω1=γ+μΨ1Λθ−Ψ4θ+μα+δ+μμ=η2κμbμvμθ+μα+δ+μγ+μσ+μvρ+μvR0−1R0−Rb+1,Ω2=γ+μΨ2Λθ−Ψ5θ+μα+δ+μμ=γ+μθ+μα+δ+μμη2η3κμb2Rb−1+σ+μvρ+μvμvRv−1η3κμb+η2τ1+η1τ2,Ω3=γ+μΨ3Λθ−Ψ6θ+μα+δ+μμ=η3θ+μα+δ+μμη32κμbRb−1+η3τ1σ+μvρ+μvμvRv−1−η2τ1+η1τ2.


### 3.7. Local Stability of the Endemic Equilibrium

Theorem 6.
The endemic equilibriumEl∗ of the model (3) is locally asymptotically stable for R0>1 (but near 1) (based on Theorem 4.1 by [24]).Proof.
Now, consider the system of equations (3). To investigate the nature of the bifurcation, we use the method introduced in Theorem 4.1 by [24]. To apply this method, the following simplification and change of variables are made first of all. Let Sh=x1, Eh=x2, Ih=x3, Rh=x4, Sv=x5, Iv=x6, Rv=x7, and Bl=x8.
By using vector notation,(66)X=x1,x2,x3,x4,x5,x6,x7,x8T.
The model system (3) can be written as
(67)dXdt=FX,
whereF=f1,f2,f3,f4,f5,f6,f7,f8T.
The following system of differential equations is established:(68)dx1dt=Λ+γx4−β2x8κ+x8+β1x6+μx1=f1,dx2dt=β2x8κ+x8+β1x6+μx1−θ+μx2=f2,dx3dt=θx2−α+δ+μx3=f3,dx4dt=δx3−γ+μx4=f4,x5dt=Π+ρx7−β3x3+μvx5=f5,dx6dt=β3x3x5−σ+μvx6=f6,dx7dt=σx6−ρ+μvx7=f7,dx8dt=τ1x3+τ2x6−μbx8=f8.
Choosingβ1∗ as bifurcation parameter and solving for β1∗ from equation (37) when R0=1, we obtained as
(69)β1∗=κμμbμvε1ε2ε4−β3ΠΛθτ1κμb+β2Πτ2β3ΠΛθκμb.
Consideringβ1=β1∗, the Jacobian of the system (68) evaluated at ε0∗l,β1∗ is given by
(70)Jε0l∗,β1∗=−μ00γ0−J10−J20−J3000J10J20θ−J40000000δ−J500000−J60−μv0ρ000−J600−J70000000σ0000τ100τ20−μb,
whereJ1=β1∗Λ/μ, J2=β2Λ/κμ, J3=θ+μ, J4=α+δ+μ, J5=γ+μ, J6=β3Π/μv, J7=σ+μv, and J8=ρ+μv.
The characteristic polynomial equation of the matrixJε0l∗,β1∗ is given by
(71)P∗λ=λλ+μλ+μvλ+γ+μλ+ρ+μvλ3+A1λ2+A2λ+A3=0,
where(72)A1=θ+μ+α+δ+μ+σ+μv+μb,A2=θ+μμb+α+δ+μ+σ+μv+μbα+δ+μ+σ+μv+α+δ+μσ+μv,A3=J7μbJ3+J4+J3J4μb1−Rb+J3J4J71−Rv,=σ+μvμbθ+2μ+α+δ+θ+μα+δ+μμbRv+θ+μα+δ+μσ+μvμbRb.
Thus,λ1=0,λ2=−μ,λ3=−μv,λ4=−γ+μ, and λ5=−ρ+μv are the five eigenvalues of Jε0l∗,β1∗. The remaining three eigenvalues are the zeros of the polynomial equation λ3+A1λ2+A2λ+A3=0, since A1>0,A2>0,A3>0,A1A2−A3>0 for R0=1 (see the proof of theorem (3.3.3)). We conclude that zero is a simple eigenvalue of Jε0l∗,β1∗ and all other eigenvalues of Jε0∗l,β1∗ have negative real numbers or real parts. This is an indication that the center manifold theorem is applicable. We first compute the left and right eigenvectors of the Jacobian matrix Jε0∗l,β1∗. let Y=y1,y2,y3,y4,y5,y6,y7,y8T be the right eigenvector of Jε0∗l,β1∗ corresponding to λ=0. Then, Y is computed as
(73)−μ00γ0−J10−J20−J3000J10J20θ−J40000000δ−J500000−J60−μv0ρ000−J600−J70000000σ0000τ100τ20−μby1y2y3y4y5y6y7y8=00000000.
Then, equation (73) can be expressed as
(74)−μy1−γy4−J1y6−J2y8=0,−J6y3−μvy5+ρy7=0,−J3y2+J1y6+J2y8=0,J6y3−J7y6=0,θy2−J4y3=0,ρy6−J8y7=0,δy3−J5y4=0,τ1y3+τ2y6−μby8=0.
The components ofY can be obtained from (74) as
(75)y1=−y3θγα+μ+μα+δ+μθ+γ+μθμγ+μ<0,y2=y3α+δ+μθ>0,y3=y3>0,y4=y3δγ+μ>0,y5=y3β3Πμvσ+ρ+μvμvσ+μvρ+μv>0,y6=y3β3Πμvσ+μv>0,y7=y3β3Πσμvσ+μvρ+μv>0,y8=y3μvτ1σ+μv+β3Πσ+μvτ2μvμbσ+μv>0.
Next, the left eigenvectorZ=z1,z2,z3,z4,z5,z6,z7,z8T of the Jacobian matrix Jε0∗l,β1∗ corresponding to λ=0 is computed as
(76)−μ00000000−J3θ0000000−J4δ−J6J60τ1γ00−J500000000−μv000−J1J1000−J4στ20000ρ0−J80−J2J200000−μbz1z2z3z4z5z6z7z8=00000000.
Thus,Z can be found as
(77)Z=0,z2,z2θ+μθ,0,0,z2μvμμbκθ+μα+δ+μ−β2Λθτ1β3Πμμbκθ,0,z2β2Λμκ,where z2>0 and it is calculated to ensure that the eigenvectors satisfy the condition Z. Y=1, for R0=1.
The local stability near the bifurcation pointβ1=β1∗ is determined by the signs of two associated constants a and b. Based on the center manifold theorem 4.1 in [24], we compute the coefficients a and b, as
(78)a=∑k,i,j=1nzkyiyj∂2fk∂xi∂xjε0l∗,β1∗,b=∑k,i=1nzkyi∂2fk∂xiβ1∗ε0l∗,β1∗.
Since the first, fourth, fifth, and seventh components ofZ are zero, we do not need the derivatives of f1,f4,f5, and f7. The derivatives of f2,f3,f6, and f8 are ∂2f2/∂x1∂x6=∂2f2/∂x6∂x1=β1, ∂2f2/∂x1∂x8=∂2f2/∂x8∂x1=β2/κ, ∂2f2/∂x6∂β1∗ε0l∗,β1∗=Λ/μ, and ∂2f2/∂x82ε0l∗,β1∗=−2β2Λ/κ2μ, and any other second order partial derivative of fk is zero, ∀k=1,2,3,4,5,6,7,8.
It follows that(79)a=∑k,i,j=1nzkyiyj∂2fk∂xi∂xjε0l∗,β1∗=z2y1y6∂2f2∂x1∂x6ε0l∗,β1∗+z2y1y8∂2f2∂x1∂x8ε0l∗,β1∗+z2y82∂2f2∂x82ε0l∗,β1∗.
Substituting the values for partial derivatives andy1,y8, we obtain
(80)a=−z2y32β3β1ΠσΦ+β2Φμvτ1σ+μv+β3Πτ2κμb+β2Λκ2μμvτ1σ+μv+β3Πτ2μvμbσ+μv2<0.
whereΦ=θγα+μ+μα+δ+μθ+γ+μ/θμvμγ+μσ+μv and
(81)b=z2y6∂2f2∂x6∂β1∗ε0l∗,β1∗=z2y3β3ΠΛμvμσ+μv>0always.
Sincea <0 and b >0 and with the help of Theorem 4.1 by Castillo-Chavez and Song in [24], it is confirmed that the system follows a forward bifurcation at R0=1 (see Figure 2), and hence, the endemic equilibrium is locally asymptotically stable for R0>1, but sufficiently close to 1
The epidemiological implication of this shows that the leptospirosis disease can be eradicated in infected population as long as whenR0<1. The condition R0<1 is necessary and sufficient for disease elimination for a model that undergoes the forward bifurcation [25].Figure 2
Forward bifurcation.


## 3.1. Positivity of Solutions

Theorem 1.
LetΩl0=Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Rv0,Bl0∈R+8∪0 be the initial condition for equation (3). Then, the set of solutions Sht,Eht,Iht,Rht,Svt,Ivt,Rvt,Blt of equation (3) is nonnegative for all t≥0.Proof.
Given that a set of nonnegative initial conditionsSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Rv0,Bl0, consider equation (3). Let t1=supt>0:Sht0>0,Eht0>0,Iht0>0,Rht0>0,Svt0>0,Ivt0>0,Rvt0>0,Blt0>0,∀t0in0,t. From the first equation of the system (3), it follows that
(5)Shdt+λht+μSh=Λ+γRht,which can be rewritten as
(6)ddtShteμt+∫0tλhτdτ=Λ+γRhteμt+∫0tλhτdτ.
Thus,(7)Sht1eμt1+∫0t1λhτdτ−Sh0=∫0t1Λ+γRhteμz+∫0zλhτdτdz.
Therefore,(8)Sht1=Sh0e−μt1−∫0t1λhτdτ+e−μt1−∫0t1λhτdτ×∫0t1eμz+∫0zλhτdτΛ+γRhzdz≥0.
Similarly, the fifth equation of the system (3) yields that
(9)Svt1=Sl0e−μlt1−∫0t1λlτdτ+e−μlt1−∫0t1λlτdτ×∫0t1eμlz+∫0zλlτdτΠ+ρRvzdz≥0.
Next, from the second equation of model (3), it follows that
(10)dEhdt≥−θ+μEh.
Applying the comparison theorem and variation of formula for (10) yields
(11)Eht≥Eh0e−θ+μt≥0,∀t≥0.
Similarly,(12)Iht≥Ih0e−α+δ+μt≥0,∀t≥0,Rht≥Rh0e−γ+μt≥0,∀t≥0,Ivt≥Iv0e−σ+μvt≥0,∀t≥0,Rvt≥Rv0e−ρ+μvt≥0,∀t≥0,Blt≥Bl0e−μbt≥0,∀t≥0.
Therefore, all solutions of the system (3) remain nonnegative for all nonnegative initial conditions.


## 3.2. Invariant Region

Theorem 2.
All solutions of the leptospirosis model equation (3) remain in Ωl=Sh,Eh,Ih,Rh,Sv,Iv,Rv,Bl∈R+8:0≤Nht≤Λ/μ,0≤Nvt≤Π/μv,0≤Blt≤τ∗/μblΛ/μ+Π/μv.Proof.
From equation (1), we get
(13)dNhdt=dShdt+dEhdt+dIhdt+dRhdt.
The first four equations of the system (3) and equation (13) give rise to
(14)dNhdt=Λ−μNh−αIh.
Since the parameter,α, and the state variable Ih are nonnegative in (14), we have
(15)dNhdt≤Λ−μNh.
LetNh∗t be the solution of the ODE dNh/dt=Λ−μNh with an initial condition. Then,
(16)dNh∗dt=Λ−μNh∗withNh∗0=Nh0∗∀t≥0.
Thus, equation (16) has a unique solution. By applying separation of variable on the differential equation (16), we get
(17)dNh∗Λ−μNh∗=dt.
Integrating (17) on both sides yields
(18)Nh∗t=Nh∗0e−μt+Λμ1−e−μt.
We now apply the comparison theorem [19] on the differential equation (15) and then Nht≤Nh∗t. Thus,
(19)Nht≤Nh∗0e−μt+Λμ1−e−μt.
Ast⟶∞, the population size Nh⟶Λ/μ.
Hence,(20)Nht≤Λμfor0≤Nh∗0≤Λμ,∀t≥0.
On the other hand, ifNh∗0>Λ/μ, the solution decrease to Λ/μ as t⟶∞.
Also, the last equation of (3) and equations (20) and (23) imply that
(21)dBldt≤τ∗Λμ+Πμv−μbBl,
whereτ∗=maxτ1,τ2.
Equation (21) yields
(22)Blt≤Bl0e−μblt+τ∗μblΛμ+Πμv1−e−μbt∀t≥0.
Hence,(23)Ωh=Sh,Eh,Ih,Rh,Bl∈R+5:0≤Nht≤Λμ,0≤Bl≤τ∗μblΛμ+Πμv,is positively invariant region.
Finally, the 5th equation of (3) implies that
(24)dNvdt≤Π−μNv,withNv0=Nv,0.
Then, equation (24) has a unique solution. Similarly, applying the comparison theorem and integrating the differential equation (21) yield
(25)Nvt=Nv0e−μvt+πμv1−e−μvt.
Therefore,(26)0≤Nvt≤Πμvfor0≤Nv0≤Πμv,Ωv=Sv,Iv,Rv∈R+3:0≤Nvt≤Πμv.
Therefore, the biologically feasible region for (3) is given by
(27)Ωl=Ωh×Ωv⊂R+5×R+3.
Therefore, every solution of the differential equation model with initial conditions inΩl remains in Ωl∀t≥0. Thus, the region Ωl is positively invariant with respect to the system (3).


## 3.3. Disease-Free Equilibrium

The disease-free equilibrium (DFE) of the model (3) only exists in the absence of infections. We compute DFE by setting left-hand side of the system (3) equal to zero. Since Eh∗=Ih∗=Rh∗=Iv∗=Rv∗=Bl∗=0 at DFE, we have
(28)Λ+γRh∗−β1Bl∗κ+Bl∗+β2Iv∗+μSh∗=0,Π+ρRv∗−β3Il∗+μvSv∗=0.Thus, the system (28) yields Sh∗=Λ/μ and Sv∗=Π/μv. Therefore, the DFE point of the model (3) is given by
(29)ε0l∗=Sh∗,0,0,0,Sv∗,0,0,0=Λμ,0,0,0,Πμv,0,0,0.To compute the basic reproduction number for our model (3), we follow the method presented in [20]. Considering only the infective compartments W=Eh,Ih,Iv,Bl in (3), then we have dW/dt=Ft−Vt,
(30)Ft=λhSh0β3SvIh0,Vt=ε1Eh−θEh+ε2Ihε4Iv−τ1Ih−τ2Iv+μblBl,where(31)λh=β2Blκ+Bl+β1Iv,ε1=θ+μ,ε2=α+δ+μ,The Jacobian matrices ofFt and of Vt at ε0l∗ are, respectively, given by
(32)F=00β1Λμβ2Πμκ00000β3Πμv000000,V=ε1000−θε20000ε300−τ1−τ2μb.Therefore, the next-generation matrix of the model is given by(33)3G=FV−1=00β1Λμβ2Λμκ00000β3Πμv0000001ε1000θε1ε21ε200001ε40θτ1ε1ε2μbτ1ε2μbτ2ε4μb1μb=y1y2y3y40000y5y6000000,where
(34)y1=β2Λθτ1κμμbε1ε2,y2=β2Λτ1κμμbε2,y3=Λβ2τ2+β1κμbκμμbε4,y4=β2Λκμμb,y5=β3Πθμvε1ε2,y6=β3Πμvε2.ConsideringG−λI=0, we have
(35)y1−λy2y3y40−λ00y5y6−λ0000−λ=0⟹λ2λ2−λy1−y3y5=0.Thus,(36)λ1=λ2=0,λ3=12y1−y12+4y3y5<0,λ4=12y1+y12+4y3y5>0,are the eigenvalues of the matrixG. Therefore, the basic reproduction number of the model is given by
(37)R0=ρG=12β2Λθτ1κμμbε1ε2+β2Λθτ1κμμbε1ε22+4β3ΠΛθβ2τ2+β1κμbκμμvμbε1ε2ε4.The basic reproduction numberR0 can be written in terms of the basic reproduction numbers corresponding to the vector transmission Rv, human transmission Rh, and environment-to-human transmission Rb as
(38)R0=12Rb+Rb2+4Rhv,whereRb=β2Λθτ1/κμμbε1ε2 and Rhv=Rh+Rv with Rh=β3β2ΠΛθτ2/κμμvμbε1ε2ε4 and Rv=β3β1ΠΛθ/μμvε1ε2ε4.


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

Theorem 3.
The disease-free equilibrium point,ε0l∗, is locally asymptotically stable whenever R0<1 and unstable otherwise.Proof.
The Jacobian matrices of the model (3) at DFE ε0l∗ is given by
(39)Jε0l∗=−μ00γ0−β1Λμ0−β2Λμκ0−ε1000β1Λμ0β2Λμκ0θ−ε20000000δ−ε3000000−β3Πμv0−μv0ρ000β3Πμv00−ε40000000σ−ε5000τ100τ20−μb,
whereε1=θ+μ,ε2=α+δ+μ,ε3=γ+μ,ε4=σ+μv, and ε5=ρ+μv.
Expanding the determinant of the characteristic equationJε0l∗−λI=0 by the first, fourth, fifth, and seventh columns, respectively, we obtain four of the eigenvalues of Jε0l∗: λ1=−μ,λ2=−γ+μ,λ3=−μvandλ4=−ρ+μv. The other four eigenvalues will be determined from the following polynomial equation:
(40)−ε1−λ0β1Λμβ2Λμκθ−ε2−λ000β3Πμv−ε4−λ00τ1τ2−μb−λ=0,⟹P4λ=λ4+A1λ3+A2λ2+A3λ+A4=0,
where(41)A1=ε1+ε2+ε4+μb,A2=ε1μb+ε1ε2+ε1ε4+ε2μb+ε4μb+ε2ε4,A3=ε1ε2μb+ε1ε4μb+ε1ε2ε4+ε2ε4μb−β2Λτ1θμκ−β1β3ΛΠθμμv=ε4μbε1+ε2+ε1ε2μb1−β2Λθτ1κμμbε1ε2+ε1ε2ε41−β1β3ΛΠθμμvε1ε2ε4=ε4μbε1+ε2+ε1ε2μb1−Rb+ε1ε2ε41−Rv,A4=ε1ε2ε4μb−β2Λθτ1ε4μκ+β1β3ΛΠθμbμμv+β2β3ΛΠθτ2κμμv=ε1ε2ε4μb−β2Λθτ1ε4μκ+β3ΠΛθβ2τ2+β1κμbκμμv=ε1ε2ε4μb1−β2Λθτ1ε4μκ+β3ΠΛθβ2τ2+β1κμbκμμv=ε1ε2ε4μb1−Rb+Rhv.
Since all the model parameters are positive, the above four eigenvalues of the matrixJε0l∗ are negative real numbers. To ensure that all the remaining four eigenvalues of the matrix have negative real parts, the Routh-Hurwitz stability criterion (see Theorem 5.1 in [21]) requires A1>0,A2>0,A3>0,A4>0,A1A2−A3>0, and A1A2A3−A32+A12A4>0. Among these, A1>0 and A2>0 are obvious.
Similarly,(42)A1A2−A3=A2ε4+μb+ε22ε1+ε4+μb+ε1ε1μb+ε1ε2+ε2μb+ε1ε4+ε2ε4+β2Λτ1θμκ+β1β3ΛΠθμμv>0.
The other three inequalities hold whenR0>1; the details are described as follows.
SinceRb<R0, 1−Rb>0 for R0<1. Similarly, 1−Rv=1−R01−Rb+R0+Rb+Rh>0 for R0<1. Thus,
(43)A3=ε4μbε1+ε2+ε1ε2μb1−Rb+ε1ε2ε41−Rv>0ifR0<1.
Next,(44)A4=ε1ε2ε4μb1−Rb+Rhv=ε1ε2ε4μb1−R01−Rb+R0>0ifR0<1.
Finally, we need to show thatA1A2A3−A32+A12A4>0. Let
(45)ϕ0=ε1μb,ϕ1=ε1ε2,ϕ2=ε1ε4,ϕ3=ε2μb,ϕ4=ε4μb,ϕ5=ε2ε4,ψ0=ε4μbε1+ε2,ψ1=ε1ε2μb1−Rb,ψ3=ε1ε2ε41−Rv.
Then,A2=ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5, and A3=ψ0+ψ1+ψ2.
Now,(46)A1A2A3=ε1ψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ε2ψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ε4ψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+μbψ0ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ψ2ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5,A32+A12A4=ψ0+ψ1+ψ22+ε1+ε2+ε4+μb2ε1ε2ε4μb1−Rb−Rhv=ψ02+ψ12+ψ22+2ψ0ψ1+ψ0ψ2+ψ1ψ2+ε12+ε22+ε42+μb2+2ε1ε2+ε1ε4+ε1μb+ε2ε4+ε2μb+ε4μbε1ε2ε4μb1−Rb−Rhv.
After manipulating the parameters algebraically, we obtain(47)A1A2A3−A32+A12A4=Ω0+Ω1,
where(48)Ω0=ε1ψ0ϕ0+ϕ1+ϕ2+ϕ4+ψ1ϕ0+ϕ1+ϕ4+ψ2ϕ1+ϕ2+ϕ4+ε2ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ψ2ϕ0+ϕ1+ϕ2+ϕ4+ϕ5+ε4ψ0ϕ1+ϕ2+ϕ4+ψ2ϕ1+ϕ2+ϕ4+ϕ5+μbψ0ϕ2+ϕ3+ϕ4+ϕ5+ψ1ϕ0+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5>0,Ω1=ε1ϕ3RbA3+ψ0+ϕ5A3Rv+ε12+ε22+A2ϕ1ϕ41+Rh+ϕ4ϕ52+ϕ02+ε2A2+ε42+μb2+A2ϕ1ϕ4Rb+Rhv>0.
Clearly,Ω0>0 and Ω1>0.
Therefore,A1A2A3−A32+A12A4>0 if R0<1. This completes the proof.


## 3.5. Global Stability of the Disease-Free Equilibrium

The method illustrated in [22, 23] is used to investigate the global asymptotic stability (GAS) of DFE point of the model. Firstly, the model (3) must be written in the form:
(49)dXdt=FX,Y,dYdt=GX,Y,GX,0=0.whereX=Sh,Rh,Sv,Rv represents the number of uninfected individuals and Y=Eh,Ih,Iv,Bl denotes the number of infected individuals. Let X∗ be the disease-free equilibrium of the system dX/dt=FX,0. Then, X∗=Λ/μ,0,Π/μv,0. The DFE point of the model ε0l∗=X∗,0=Λ/μ,0,0,0,Π/μv,0,0,0 is guaranteed to be GAS if R0<1 (which is locally asymptotically stable (LAS)), and the following two conditions H1 and H2 hold:H1: For dX/dt=FX,0, X∗ is globally asymptotically stable.H2: GX,Y=AY−G⋆X,Y,G⋆X,Y≥0,∀X,Y∈Ωl, where A=∂Gε0l∗/∂Y is an M-matrix and Ωl is the region where the model makes biological sense. If the model equation (3) satisfies the above two conditions, then the following theorem holds.Theorem 4.
The disease-free equilibrium point,ε0l∗=Λ/μ,0,0,0,Π/μv,0,0,0, is GAS for the model (3) provided that R0<1 (LAS) and conditions H1 and H2 hold.Proof.
We only need to show that the conditions H1 and H2 hold when R0<1. Considering the system (3), we have
(50)FX,0=Sh′tRh′tSv′tRv′t=Λ+γRh−μSh−γ+μRhΠ+ρRv−μvSv−ρ+μvRv,X∗=Λμ,0,Πμv,0.
Note that the second and the fourth equations of the system (50) are first-order linear ODEs, and their solutions are given by, respectively,
(51)Rht=Rh,0e−γ+μt,(52)Rvt=Rv,0e−ρ+μvt.
Comparing equation (51) and the first equation of the system (50) yields
(53)S′ht=Λ+γRv,0e−ρ+μv−μSht.
Solving the first-order linear ODE (53), its solution is found as
(54)Sht=Λμ+Sh,0e−μt−Rh,0e−γ+μt.
Similarly,(55)Sv′t=Π+ρRv,0e−ρ+μvt−μvSvt,and its solution is given by
(56)Svt=Πμv+Sv,0e−μvt−Rv,0e−ρ+μvt.
Now, suppose that the time,t⟶∞, we need to show that X⟶X∗. Obviously, Rht⟶0, Rvt⟶0 as t⟶∞. Thus, Sht⟶Λ/μ, Svt⟶Π/μv as t⟶∞, regardless of the values of Sh0, Rh0 and Sv0, Rv0, respectively. Thus, all points with respect to this condition converge at X0∗=Λ/μ,0,Π/μv,0. Hence, X∗=Λ/μ,0,Π/μv,0 is GAS.
Next, we have(57)GX,Y=G1X,YG2X,YG3X,YG4X,Y=β2Blκ+Bl+β1IvSh−θ+μEhθEh−α+δ+μIhβ3SvIh−σ+μvIvτ1Ih+τ2Iv−μbBl.
We can then obtain(58)A=∂GX∗,0∂Y=−θ+μ0β1Λμβ2Λκμθ−α+δ+μ000β3Πμv−σ+μv00τ1τ2−μb,which is clearly an M-matrix
(59)G⋆X,Y=AY−GX,Y=β1IvΛμ−Sh+β2BlκΛμ−κShκ+Bl0β3IhΠμv−Sv0.
Since all parameters are nonnegative,0≤κSh/κ+Bl≤Sh≤Λ/μ∵κ/κ+Bl≤1, 0≤Sv≤Π/μv. It follows that Λ/μ−Sh≥0, Λ/μ−κSh/κ+Bl≥0, and Π/μv−Sv≥0.
Hence,G⋆X,Y≥0∀X,Y∈Γ. Therefore, the disease-free equilibrium point, ε0l∗=Λ/μ,0,0,0,Π/μv,0,0,0, of the model (3) is globally asymptotically stable. This completes the proof.Remark 5.
The model (3) does not exhibit a backward bifurcation at R0=1 when R0<1, since DFE is the only stable (or positive) equilibrium point for R0<1 in the Theorem 4.The backward bifurcation has been observed in the leptospirosis epidemic model [16].


## 3.6. Endemic Equilibrium (EE)

The EE point of our model is computed by setting the system of differential equations in (3) to zero. Thus, it is obtained by dSh∗/dt=dEh∗/dt=dIh∗/dt=dRh∗/dt=dSv∗/dt=dIv∗/dt=dRv∗/dt=dBl∗/dt=0, which is given by El∗=Sh∗,Eh∗,Ih∗,Rh∗,Sv∗,Iv∗,Rv∗,Bl∗, where all disease states are considered to be positive, and thus, Ih∗ must be greater than zero for all the other states to be positive. Now, by setting the system of equations in (3) to zero at steady state, we have
(60)dSh∗dt=Λ+γRh∗−β2Bl∗κ+Bl∗+β1Iv∗+μSh∗=0,dEhdt=β2Bl∗κ+Bl∗+β1Iv∗Sh∗−θ+μEh∗=0,dIh∗dt=θEh∗−α+δ+μIh∗=0,dRh∗dt=δIh∗−γ+μRh∗=0,dSv∗dt=Π+ρRv∗−β3Ih∗+μvSv∗=0,dIv∗dt=β3Ih∗Sv∗−σ+μvIv∗=0,dRv∗dt=σIv∗−ρ+μvRv∗=0,dBl∗dt=τ1Ih∗+τ2Iv∗−μbBl∗=0.Due to the complexity of the model system (3), all other state variables at the steady state are expressed in terms of the steady state of Ih which is denoted as Ih∗. Accordingly, we found as
(61)El∗=Sh∗,Eh∗,Ih∗,Rh∗,Sv∗,Iv∗,Rv∗,Bl∗,=Λ+γRh∗λ∗+μ,ε2θIh∗,Ih∗,δε3Ih∗,Π+ρRv∗β3Ih∗+μv,β3Πε5Ih∗μvε4ε5+β3ε6Ih∗,β3ΠσIh∗μvε4ε5+β3ε6Ih∗,τ1Ih∗+τ2Iv∗μb,where(62)ε1=θ+μ,ε2=α+δ+μ,ε3=γ+μ,ε4=σ+μv,ε5=ρ+μv,ε6=σ+ρ+μv.By combining the equations of (61) and the first and second equations of the model system (3), we obtained a septic (or heptic) polynomial equation for Ih∗ which is described as follows:
(63)PIh∗=Ih∗AIh∗6+BIh∗5+CIh∗4+DIh∗3+EIh∗2+FIh∗+G=0,where(64)A=−Ψ3Ψ7η4+Ψ72θ+μα+δ+μγ+μμ<0always,B=−η4Ψ2Ψ7+Ψ3Ψ6+Ψ7η3θ+μα+δ+μγ+μμη32κμbRb−1+η3τ1σ+μvρ+μvμvRv−1−η3η2τ1+η1τ2−Ψ6<0ifR0<1,C=−η4Ψ1Ψ7+Ψ2Ψ6+Ψ3Ψ5+Ψ7Ω2−Ψ7Ψ5θ+μα+δ+μγ+μμ+Ψ6Ω3=−η4Ψ1Ψ7+Ψ2Ψ6+Ψ3Ψ5+Ψ6Ω3+γ+μθ+μα+δ+μμη2η3κμb2Rb−1+σ+μvρ+μvμvRv−2η3κμb+η2τ1+η1τ2<0ifR0<1,D=−η4Ψ1Ψ6+Ψ2Ψ5+Ψ3Ψ4+Ψ7Ω1−Ψ4θ+μα+δ+μγ+μμ+Ψ6Ω2+Ψ5Ω3,E=−η4Ψ1Ψ5+Ψ2Ψ4+Ψ6Ω1+Ψ5Ω2+Ψ4Ω3,F=−η4Ψ1Ψ4+Ψ5Ω1+Ψ4Ω2,G=Ψ4Ω1<0ifR0<1,with(65)η1=Πβ3ρ+μv,η2=μvσ+μvρ+μv,η3=β3μvσ+ρ+μv,η4=θγα+μ+μα+δ+μθ+γ+μ,Ψ1=η2β2η2τ1+η1τ2+β1κη1μb,Ψ2=η2τ1+η1τ2β2η3+β1η1+η3β2η2τ1+β1κη1μb,Ψ3=η3τ1β2η3+β1η1,Ψ4=κη22μb,Ψ5=2κη2η3μb+η2η2τ1+η1τ2,Ψ6=η3κη3μb+2η2τ1+η1τ2,Ψ7=η32τ1,Ω1=γ+μΨ1Λθ−Ψ4θ+μα+δ+μμ=η2κμbμvμθ+μα+δ+μγ+μσ+μvρ+μvR0−1R0−Rb+1,Ω2=γ+μΨ2Λθ−Ψ5θ+μα+δ+μμ=γ+μθ+μα+δ+μμη2η3κμb2Rb−1+σ+μvρ+μvμvRv−1η3κμb+η2τ1+η1τ2,Ω3=γ+μΨ3Λθ−Ψ6θ+μα+δ+μμ=η3θ+μα+δ+μμη32κμbRb−1+η3τ1σ+μvρ+μvμvRv−1−η2τ1+η1τ2.


## 3.7. Local Stability of the Endemic Equilibrium

Theorem 6.
The endemic equilibriumEl∗ of the model (3) is locally asymptotically stable for R0>1 (but near 1) (based on Theorem 4.1 by [24]).Proof.
Now, consider the system of equations (3). To investigate the nature of the bifurcation, we use the method introduced in Theorem 4.1 by [24]. To apply this method, the following simplification and change of variables are made first of all. Let Sh=x1, Eh=x2, Ih=x3, Rh=x4, Sv=x5, Iv=x6, Rv=x7, and Bl=x8.
By using vector notation,(66)X=x1,x2,x3,x4,x5,x6,x7,x8T.
The model system (3) can be written as
(67)dXdt=FX,
whereF=f1,f2,f3,f4,f5,f6,f7,f8T.
The following system of differential equations is established:(68)dx1dt=Λ+γx4−β2x8κ+x8+β1x6+μx1=f1,dx2dt=β2x8κ+x8+β1x6+μx1−θ+μx2=f2,dx3dt=θx2−α+δ+μx3=f3,dx4dt=δx3−γ+μx4=f4,x5dt=Π+ρx7−β3x3+μvx5=f5,dx6dt=β3x3x5−σ+μvx6=f6,dx7dt=σx6−ρ+μvx7=f7,dx8dt=τ1x3+τ2x6−μbx8=f8.
Choosingβ1∗ as bifurcation parameter and solving for β1∗ from equation (37) when R0=1, we obtained as
(69)β1∗=κμμbμvε1ε2ε4−β3ΠΛθτ1κμb+β2Πτ2β3ΠΛθκμb.
Consideringβ1=β1∗, the Jacobian of the system (68) evaluated at ε0∗l,β1∗ is given by
(70)Jε0l∗,β1∗=−μ00γ0−J10−J20−J3000J10J20θ−J40000000δ−J500000−J60−μv0ρ000−J600−J70000000σ0000τ100τ20−μb,
whereJ1=β1∗Λ/μ, J2=β2Λ/κμ, J3=θ+μ, J4=α+δ+μ, J5=γ+μ, J6=β3Π/μv, J7=σ+μv, and J8=ρ+μv.
The characteristic polynomial equation of the matrixJε0l∗,β1∗ is given by
(71)P∗λ=λλ+μλ+μvλ+γ+μλ+ρ+μvλ3+A1λ2+A2λ+A3=0,
where(72)A1=θ+μ+α+δ+μ+σ+μv+μb,A2=θ+μμb+α+δ+μ+σ+μv+μbα+δ+μ+σ+μv+α+δ+μσ+μv,A3=J7μbJ3+J4+J3J4μb1−Rb+J3J4J71−Rv,=σ+μvμbθ+2μ+α+δ+θ+μα+δ+μμbRv+θ+μα+δ+μσ+μvμbRb.
Thus,λ1=0,λ2=−μ,λ3=−μv,λ4=−γ+μ, and λ5=−ρ+μv are the five eigenvalues of Jε0l∗,β1∗. The remaining three eigenvalues are the zeros of the polynomial equation λ3+A1λ2+A2λ+A3=0, since A1>0,A2>0,A3>0,A1A2−A3>0 for R0=1 (see the proof of theorem (3.3.3)). We conclude that zero is a simple eigenvalue of Jε0l∗,β1∗ and all other eigenvalues of Jε0∗l,β1∗ have negative real numbers or real parts. This is an indication that the center manifold theorem is applicable. We first compute the left and right eigenvectors of the Jacobian matrix Jε0∗l,β1∗. let Y=y1,y2,y3,y4,y5,y6,y7,y8T be the right eigenvector of Jε0∗l,β1∗ corresponding to λ=0. Then, Y is computed as
(73)−μ00γ0−J10−J20−J3000J10J20θ−J40000000δ−J500000−J60−μv0ρ000−J600−J70000000σ0000τ100τ20−μby1y2y3y4y5y6y7y8=00000000.
Then, equation (73) can be expressed as
(74)−μy1−γy4−J1y6−J2y8=0,−J6y3−μvy5+ρy7=0,−J3y2+J1y6+J2y8=0,J6y3−J7y6=0,θy2−J4y3=0,ρy6−J8y7=0,δy3−J5y4=0,τ1y3+τ2y6−μby8=0.
The components ofY can be obtained from (74) as
(75)y1=−y3θγα+μ+μα+δ+μθ+γ+μθμγ+μ<0,y2=y3α+δ+μθ>0,y3=y3>0,y4=y3δγ+μ>0,y5=y3β3Πμvσ+ρ+μvμvσ+μvρ+μv>0,y6=y3β3Πμvσ+μv>0,y7=y3β3Πσμvσ+μvρ+μv>0,y8=y3μvτ1σ+μv+β3Πσ+μvτ2μvμbσ+μv>0.
Next, the left eigenvectorZ=z1,z2,z3,z4,z5,z6,z7,z8T of the Jacobian matrix Jε0∗l,β1∗ corresponding to λ=0 is computed as
(76)−μ00000000−J3θ0000000−J4δ−J6J60τ1γ00−J500000000−μv000−J1J1000−J4στ20000ρ0−J80−J2J200000−μbz1z2z3z4z5z6z7z8=00000000.
Thus,Z can be found as
(77)Z=0,z2,z2θ+μθ,0,0,z2μvμμbκθ+μα+δ+μ−β2Λθτ1β3Πμμbκθ,0,z2β2Λμκ,where z2>0 and it is calculated to ensure that the eigenvectors satisfy the condition Z. Y=1, for R0=1.
The local stability near the bifurcation pointβ1=β1∗ is determined by the signs of two associated constants a and b. Based on the center manifold theorem 4.1 in [24], we compute the coefficients a and b, as
(78)a=∑k,i,j=1nzkyiyj∂2fk∂xi∂xjε0l∗,β1∗,b=∑k,i=1nzkyi∂2fk∂xiβ1∗ε0l∗,β1∗.
Since the first, fourth, fifth, and seventh components ofZ are zero, we do not need the derivatives of f1,f4,f5, and f7. The derivatives of f2,f3,f6, and f8 are ∂2f2/∂x1∂x6=∂2f2/∂x6∂x1=β1, ∂2f2/∂x1∂x8=∂2f2/∂x8∂x1=β2/κ, ∂2f2/∂x6∂β1∗ε0l∗,β1∗=Λ/μ, and ∂2f2/∂x82ε0l∗,β1∗=−2β2Λ/κ2μ, and any other second order partial derivative of fk is zero, ∀k=1,2,3,4,5,6,7,8.
It follows that(79)a=∑k,i,j=1nzkyiyj∂2fk∂xi∂xjε0l∗,β1∗=z2y1y6∂2f2∂x1∂x6ε0l∗,β1∗+z2y1y8∂2f2∂x1∂x8ε0l∗,β1∗+z2y82∂2f2∂x82ε0l∗,β1∗.
Substituting the values for partial derivatives andy1,y8, we obtain
(80)a=−z2y32β3β1ΠσΦ+β2Φμvτ1σ+μv+β3Πτ2κμb+β2Λκ2μμvτ1σ+μv+β3Πτ2μvμbσ+μv2<0.
whereΦ=θγα+μ+μα+δ+μθ+γ+μ/θμvμγ+μσ+μv and
(81)b=z2y6∂2f2∂x6∂β1∗ε0l∗,β1∗=z2y3β3ΠΛμvμσ+μv>0always.
Sincea <0 and b >0 and with the help of Theorem 4.1 by Castillo-Chavez and Song in [24], it is confirmed that the system follows a forward bifurcation at R0=1 (see Figure 2), and hence, the endemic equilibrium is locally asymptotically stable for R0>1, but sufficiently close to 1
The epidemiological implication of this shows that the leptospirosis disease can be eradicated in infected population as long as whenR0<1. The condition R0<1 is necessary and sufficient for disease elimination for a model that undergoes the forward bifurcation [25].Figure 2
Forward bifurcation.


## 4. Sensitivity Analysis

In this section, we investigated the sensitivity of the parameters for the basic reproduction of the model using the idea presented in [26–28]. It is important to carry out the sensitivity of the basic reproduction number R0 for its parameters. This will give parameters with a high impact on the leptospirosis model (3) and therefore allow to target on control measures to reduce the transmission of the disease. To measure the sensitivity index of R0 to a given parameter β, we use the following relation:
(82)ΩβR0=∂R0∂ββR0.An analytical expression for the sensitivity index of each parameter involved inR0 is derived as follows:
(83)Ωβ1R0=RvR02R0−Rb>0∵R0>Rb,Ωβ2R0=Rb2R0−Rb+Rb2+2Rv2R02R0−Rb>0,ΩΛR0=Rb2R0−Rb+Rb2+2Rhv2R02R0−Rb>0,ΩΠR0=RhvR02R0−Rb>0,ΩμR0=−Rb2R0−Rb+Rb2+2Rhvε1ε2+με1+ε22R02R0−Rb<0,Ωβ3R0=RhvR02R0−Rb>0,ΩθR0=μRb2R0−Rb+Rb2+2Rhv2R02R0−Rb>0,ΩμvR0=−ε4+μvRhvR02R0−Rb<0,ΩαR0=−αε2Rb2R0−Rb+Rb2+2Rhv2R02R0−Rb<0,Ωτ1R0=Rb2R0−Rb+Rb22R02R0−Rb>0,ΩδR0=−δε2Rb2R0−Rb+Rb2+2Rhv2R02R0−Rb<0,ΩσR0=−σRhvε4R02R0−Rb<0,Ωτ2R0=RhR02R0−Rb>0,ΩκR0=−RbR0+RhR02R0−Rb<0,ΩμbR0=−RbR0+RhR02R0−Rb<0.The values of the sensitivity indices ofR0 to each of its parameter values are computed and listed in Table 3 in decreasing order of sensitivity indices from the most to the least sensitive parameter. We used an algorithm Markov Chain Monte Carlo (MCMC) to fit parameter values of the model. The parameter values were used from previously published papers, and others were estimated as shown in Table 3.Table 3
Sensitivity indices ofR0 to each of the parameter values.
ParameterValueSensitivity indexΛ0.270.5421μv0.0029−0.0306μ0.0009−0.0058Π20.4579β30.00070.4579σ0.064−0.438β10.000330.4167δ0.072−0.3457α0.04−0.1921β20.08150.5κ10000−0.1254τ10.060.0842τ20.20.0412μb0.05−0.1254θ0.0920.0004879The negative and positive signs of sensitivity index values in Table3 represent the negative and positive (respectively) influence of parameters on R0. Thus, from the table, it can be observed that R0 increases when the values of Λ,Π,β1,β2,β3,θ,τ1,τ2are increase and vice versa. On the other hand, increasing the values of μ,μv,μb,σ,δ,α,κ decreases R0 and vice versa. From the results of sensitivity analysis, the parameters Λ,Π,β1,β3,μ,μv,δ are the most influencing on R0. Further, deceasing Λ by 10% decreases R0 by 5.54%. Increasing Π or β3 by 10% increases R0 by 4.47%. In contrast, increasing δ by 10% decreases R0 by 3.2%. Similarly, decreasing (or increasing) μ and μv by a given percentage increases (or decreases, respectively) R0 by half of the percentage. Epidemiological implications of this is that preventive and control efforts should be targeted on parameters Λ,Π,β1,β3,μv,δ to reduce the burden of the disease in infected population since it is not reasonable biologically or ethically to use the mortality rate of humans as a control measure. In particular, reducing the total number of rodent population can minimize the transmission of disease, since μv is the reduction rate of rodents and it has a highly negative influence on R0. Further, the effect of the most sensitive parameters on disease dynamics behavior is illustrated graphically in the next section.


## 5. Numerical Simulations and Discussions

In this section, we perform some numerical results for the leptospirosis model (3) to demonstrate the results of the qualitative analysis of the model that has been already performed in the previous sections. To do this, we used the fourth-order Runge-Kutta method in the MATLAB program. The parameter values used in the simulations are given in Table 4. These parameter values have taken from previously published papers, and others are estimated as shown in Table 4.Table 4
Description of parameter values used in the simulations.
ParameterDescriptionValueUnitSourceΛThe recruitment rate of susceptible humansμ×N0Humans day-1AssumedΠThe recruitment rate of susceptible rodents2Rodents day-1[5]β1The transmission coefficient betweenSh and Iv0.00033(Rodents day)-1Assumedβ2The transmission coefficient betweenSh and Bb0.0815Day-1Assumedβ3The transmission coefficient betweenSv and Ih0.0007(Humans day)-1AssumedμThe natural death rate of susceptible humans0.0009Day-1[29]; [5]αThe disease-related death rate of humans0.04Day-1AssumedθThe rate at which latently infected individual transfer to infectious0.092Day-1[30]γLeptospirosis disease waning immunity fromRh0.089Day-1[31]δThe rate at which infectious individual transfer to human recovery class0.072Day-1AssumedσThe rate of recovery from infectious rodent0.064Day-1AssumedρLeptospirosis disease waning immunity fromRv0.083Day-1AssumedμvThe natural death rate of the rodents0.0029Day-1[29]μbThe natural death rate of the bacteria0.05Day-1[6]τ1The rate of increase of bacterial population due toIh0.06No.ofpathogenhumansdayAssumedτ2The rate of increase of bacterial population due toIv0.2No.ofpathogenrodentsdayAssumedκThe concentration of pathogenic population10000No.ofpathogenAssumedThe initial conditions are chosen asSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Rv0,Bm0=270,20,10,0,510,10,0,100. Using equation (63), parameter values of the Table 4 and the initial conditions, we found that R0≈2.8552>1, a unique positive endemic equilibrium which given by El∗=Sh∗,Eh∗,Ih∗,Rh∗,Sv∗,Iv∗,Rv∗,Bl∗≈47.668,6.5229,5.3154,4.2571,628.6425,34.9632,26.0494,146.2314 with unique positive value Ih∗≈5.3154; the other six roots of the polynomial equation (63) are obtained as I0=0,I1≈−9658.1994,I2≈−6509.1639,I3≈−54.7671,I4≈−47.8513,I5≈−47.2545.From graphical illustration of stability analysis of EE, we observed that all trajectories of the model solutions eventually move towards the steady-stateEl∗ as shown the in Figures 3(a)–3(i) since El∗ the unique positive EE is locally asymptotically stable for R0>1. This shows that the disease persists in population because R0=2.8552>1. On the other hand, if changing the values of β1,β2,β3,μv, and δ into 0.0001, 0.0515, 0.0003, 0.0035, and 0.092, respectively, then it yields R0≈0.932<1. In this case, the model has the DFE which is globally asymptotically stable, and hence, El∗ is unstable. Thus, all trajectories of the model solutions tend to the DFE over time as shown in Figures 4(a)–4(f). Biologically, the implication is that the disease will die out from the community over time, while the susceptible humans and susceptible rodents eventually approach their maximum value for this particular study. Therefore, we should reduce the value of R0 as much as possible to eradicate the disease rapidly.Figure 3
Plots showing the stability behavior ofEE, E1∗=47.668,6.5229,5.3154,4.2571,628.6425,34.9632,26.0494,146.2314 of the model (3) as t varies.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 4
Plots showing the stability behavior ofDFE,ε0l∗=300,0,0,0,571.4,0,0,0 of the model (3) as time changes.
(a)(b)(c)(d)(e)(f)


### 5.1. Effect of the Transmission Rates of Humans and Rodents (β1 and β3) on Dynamics Behavior of the Disease

The transmission ratesβ1 and β3 play a significant role on expansion of the disease transmission by increasing the incidence rates of humans and rodents in the infected population. In this section, we illustrate the sensitivity of the model (3) in human, rodent, and bacterial populations using different values of the parameters. In Figures 5(a) and 5(b), it can be seen that decreasing the human transmission rate β1 increases the sizes of susceptible humans and susceptible rodents (eventually) and vice versa. On the other hand, reducing this parameter reduces the sizes of infected individuals, infected rodents, recovered rodents, and bacterial population from the beginning and eventually as shown Figures 5(c)–5(f), respectively. This was performed by taking the values of β1 as β1=0.00017,β1=0.00023,β1=0.00033,β1=0.0005, and β1=0.0008. In the same manner, we observed the effects for the rodent transmission rate β3 by setting its values as β3=0.0005,β3=0.0006,β3=0.0007,β3=0.00085, and β3=0.00099 as shown in Figure 6. A decrease (or increase) in the value of human transmission rate β3 will cause decrease (or increase, respectively) in the size of the infected individuals, infected rodents, and recovered rodents as indicated in Figures 6(b), 6(d), and 6(e), respectively. Similarly, in Figure 6(f), it can be seen that decreasing β3 decreases the growth bacterial population in an environment with time as they are directly proportional. However, decreasing (or increasing) the human transmission rate β3 increases (or decreases, respectively) the size of the number of susceptible humans and susceptible rodents as shown in Figures 6(a) and 6(a), respectively. In view of the simulations of the model in Figures 5 and 6, the smaller parameter values (as possible) will result in the smaller R0 value and will decrease the incidence rates of human and rodent infections. This shows that reducing the expansion of human infection depends on the reduction of the transmission rates β1 and β3 on susceptible individuals and susceptible rodents. Therefore, by minimizing the values of β1 and β3 significantly, the number of infected humans, infected rodents, and the growth of leptospires in the environment can be reduced and could finally vanish. This can be achieved by treating infected individuals (using treatment control), avoiding contact with the contaminated environment (soil or water), or by reducing infected rodents as well as the total number of rodent population.Figure 5
Plots showing the sensitivity of model (3) using different values of β1.
(a)(b)(c)(d)(e)(f)Figure 6
Plots showing the sensitivity of model (3) using different values of β3.
(a)(b)(c)(d)(e)(f)


### 5.2. Effect of the Human Recovery Rate (δ) on Dynamics Behavior of the Disease

The infectious individuals can be reduced by increasing the value of the human recovery rateδ using control interventions (treatment of disease infection). In this section, we simulate the model to observe the effect of δ in human, rodent, and bacterial populations by varying the value δ. From Figure 7, we observed that δ is directly proportional with susceptible and recovered individuals as revealed in Figures 7(a) and 7(c). Moreover, it can be seen in Figures 7(b), 7(e), and 7(f) that the population of infected humans and infected rodents and population of bacterial increase (or decrease) from the beginning and eventually as δ decreases (or increases, respectively), whereas the size of susceptible rodents increases (or decreases) from the beginning and eventually as δ increases (or decreases, respectively) with time as shown in Figure 7(d). From an epidemiological perspective, human leptospirosis infection can be minimized by reducing the number of infected individuals significantly using treatment control effort on infected individuals.Figure 7
Plots showing the sensitivity of model (3) using different values of δ.
(a)(b)(c)(d)(e)(f)


### 5.3. Effect of the Mortality Rate of Rodents (μv)

From Figures5(b), 6(c), and 7(d), increasing the number of susceptible rodents depends on the reduction of the parameters β1,β3, and δ on susceptible humans, susceptible rodents, and infected individuals, respectively. However, reducing these parameters may cause an increase in the total number of rodent population slowly (as can be seen in the figures). Thus, the control mechanism for growth of the rodent population (for instance, rodenticide) should be implemented to reduce the total number of rodent population. The effect of mortality rate of rodents μv in human, rodent, and bacterial population is illustrated in this section by using the different values of μv as μv=0.001,μv=0.0029,μv=0.0059,μv=0.009, and μv=0.014. It can be observed, from Figures 8(b)–8(f), that the number of infected individuals, the total number of rodents (SV+Iv+RV), and the size of bacterial population decrease as μv increases. Conversely, decreasing the mortality rate of rodents decreases the number of susceptible humans as shown in Figure 8(a).Figure 8
Plots showing the sensitivity of model (3) using different values of μv.
(a)(b)(c)(d)(e)(f)Therefore, controlling the major direct zoonotic infections in a human population depends on the reduction of animal reservoirs, as in the case of leptospirosis infection. Reduction of the rodent population may also reduce the number of human infections as well as the concentration ofLeptospira interrogans in the environment. Thus, the rodent born leptospirosis infections of humans can be controlled by implementing rodents’ reduction mechanism. Further, since R0≈2.8552>1, all trajectories of the model solutions eventually move towards the unique positive El∗. Thus, it should be intensified to use optimal efforts on interventions to eliminate the disease in the infected population as much as possible.


## 5.1. Effect of the Transmission Rates of Humans and Rodents (β1 and β3) on Dynamics Behavior of the Disease

The transmission ratesβ1 and β3 play a significant role on expansion of the disease transmission by increasing the incidence rates of humans and rodents in the infected population. In this section, we illustrate the sensitivity of the model (3) in human, rodent, and bacterial populations using different values of the parameters. In Figures 5(a) and 5(b), it can be seen that decreasing the human transmission rate β1 increases the sizes of susceptible humans and susceptible rodents (eventually) and vice versa. On the other hand, reducing this parameter reduces the sizes of infected individuals, infected rodents, recovered rodents, and bacterial population from the beginning and eventually as shown Figures 5(c)–5(f), respectively. This was performed by taking the values of β1 as β1=0.00017,β1=0.00023,β1=0.00033,β1=0.0005, and β1=0.0008. In the same manner, we observed the effects for the rodent transmission rate β3 by setting its values as β3=0.0005,β3=0.0006,β3=0.0007,β3=0.00085, and β3=0.00099 as shown in Figure 6. A decrease (or increase) in the value of human transmission rate β3 will cause decrease (or increase, respectively) in the size of the infected individuals, infected rodents, and recovered rodents as indicated in Figures 6(b), 6(d), and 6(e), respectively. Similarly, in Figure 6(f), it can be seen that decreasing β3 decreases the growth bacterial population in an environment with time as they are directly proportional. However, decreasing (or increasing) the human transmission rate β3 increases (or decreases, respectively) the size of the number of susceptible humans and susceptible rodents as shown in Figures 6(a) and 6(a), respectively. In view of the simulations of the model in Figures 5 and 6, the smaller parameter values (as possible) will result in the smaller R0 value and will decrease the incidence rates of human and rodent infections. This shows that reducing the expansion of human infection depends on the reduction of the transmission rates β1 and β3 on susceptible individuals and susceptible rodents. Therefore, by minimizing the values of β1 and β3 significantly, the number of infected humans, infected rodents, and the growth of leptospires in the environment can be reduced and could finally vanish. This can be achieved by treating infected individuals (using treatment control), avoiding contact with the contaminated environment (soil or water), or by reducing infected rodents as well as the total number of rodent population.Figure 5
Plots showing the sensitivity of model (3) using different values of β1.
(a)(b)(c)(d)(e)(f)Figure 6
Plots showing the sensitivity of model (3) using different values of β3.
(a)(b)(c)(d)(e)(f)


## 5.2. Effect of the Human Recovery Rate (δ) on Dynamics Behavior of the Disease

The infectious individuals can be reduced by increasing the value of the human recovery rateδ using control interventions (treatment of disease infection). In this section, we simulate the model to observe the effect of δ in human, rodent, and bacterial populations by varying the value δ. From Figure 7, we observed that δ is directly proportional with susceptible and recovered individuals as revealed in Figures 7(a) and 7(c). Moreover, it can be seen in Figures 7(b), 7(e), and 7(f) that the population of infected humans and infected rodents and population of bacterial increase (or decrease) from the beginning and eventually as δ decreases (or increases, respectively), whereas the size of susceptible rodents increases (or decreases) from the beginning and eventually as δ increases (or decreases, respectively) with time as shown in Figure 7(d). From an epidemiological perspective, human leptospirosis infection can be minimized by reducing the number of infected individuals significantly using treatment control effort on infected individuals.Figure 7
Plots showing the sensitivity of model (3) using different values of δ.
(a)(b)(c)(d)(e)(f)


## 5.3. Effect of the Mortality Rate of Rodents (μv)

From Figures5(b), 6(c), and 7(d), increasing the number of susceptible rodents depends on the reduction of the parameters β1,β3, and δ on susceptible humans, susceptible rodents, and infected individuals, respectively. However, reducing these parameters may cause an increase in the total number of rodent population slowly (as can be seen in the figures). Thus, the control mechanism for growth of the rodent population (for instance, rodenticide) should be implemented to reduce the total number of rodent population. The effect of mortality rate of rodents μv in human, rodent, and bacterial population is illustrated in this section by using the different values of μv as μv=0.001,μv=0.0029,μv=0.0059,μv=0.009, and μv=0.014. It can be observed, from Figures 8(b)–8(f), that the number of infected individuals, the total number of rodents (SV+Iv+RV), and the size of bacterial population decrease as μv increases. Conversely, decreasing the mortality rate of rodents decreases the number of susceptible humans as shown in Figure 8(a).Figure 8
Plots showing the sensitivity of model (3) using different values of μv.
(a)(b)(c)(d)(e)(f)Therefore, controlling the major direct zoonotic infections in a human population depends on the reduction of animal reservoirs, as in the case of leptospirosis infection. Reduction of the rodent population may also reduce the number of human infections as well as the concentration ofLeptospira interrogans in the environment. Thus, the rodent born leptospirosis infections of humans can be controlled by implementing rodents’ reduction mechanism. Further, since R0≈2.8552>1, all trajectories of the model solutions eventually move towards the unique positive El∗. Thus, it should be intensified to use optimal efforts on interventions to eliminate the disease in the infected population as much as possible.


## 6. Conclusion

In this work, we presented and analyzed a deterministic mathematical model for the transmission process of leptospirosis disease, which includes human, vector (rodent), and bacterial populations. Firstly, the qualitative properties of the model are studied. Based on the next-generation matrix approach, the basic reproduction numberR0 of the formulated model is computed to perform sensitivity analysis of the model parameters. From the stability analysis of the model, we found that the disease-free equilibrium (DFE) is globally asymptotically stable if R0<1 and unstable otherwise. The local stability of an endemic equilibrium is performed based on the general center manifold theory, and EE is locally asymptotically stable when R0<1. The model exhibits forward bifurcation. This shows that the leptospirosis disease can be eradicated in infected population as long as when R0<1. The most sensitive parameters for R0 are Λ,Π,μ,β1,β3,μv, and δ, among other parameters. This illustrates that minimizing the number of infected individuals depends on the reduction of β1 on susceptible individuals and the increase of recovery rate δ. This shows that minimizing the human infection depends on the reduction of β3 on susceptible vectors and the increase of mortality rate of rodents μv. Epidemiological implications of this illustrate that preventive and control efforts should be targeted on the parameters β1, β3, and μv to reduce the burden of the disease in human populations since it is not reasonable ethically or biologically to use the death rate of humans as a control measure. Finally, the model is simulated to show the changing effects of the most influencing parameters on the disease expansion as well as the stability behaviors of the steady states. From graphical stability behavior of EE, all trajectories of the model solutions eventually move to the unique EE for R0≈2.8552>1. In this case, the disease persists in the population. On the other hand, all trajectories of the model solutions gradually move towards the DFE for R0<1. The biological implication of this shows that the disease will die from a community as long as when the basic reproduction number of the model is less than unity. Therefore, the value of R0 should be reduced as much as possible to eradicate the disease rapidly in an infected population. Moreover, our numerical results of the model demonstrate that reducing the transmission of human leptospirosis infection can be achieved by implementing the presentation and control interventions. Thus, our study recommends the following intervention efforts to minimize the risk factors and control the disease infection that should be implemented in human, rodent, and bacterial populations in an infected community:
(i)
Reduce the rodent population using rodenticide or resource reduction. Implementation of this intervention also reduces the shedding rate of leptospire into the environment(ii)
Prevention interventions on susceptible humans. Avoid contact with infected rodents and rodents urine-contaminated environment (floodwater, lakes, rivers, or soil) using an appropriate personal protective equipment (PPE) such as wearing rubber boots, waterproof dressings to cover wounds or skin, goggles, and rubber gloves as well as maintaining good personal hygiene. These practices will minimize the acquisition of the infection in the community(iii)
Use the treatment control on infected individuals. Although the leptospirosis infection shows nonfixed symptoms and signs initially, its clinical manifestation can be detected by microbiological laboratory tests. This control mechanism can minimize the expansion of the disease by treating infected individuals through antibiotics that are recommended for treatment of the disease in the infected population following clinical rules(iv)
Increase mortality rate ofLeptospira interrogans in environment (especially, urban slum environment) using appropriate environmental modificationsThese clinical interpretations (interventions) can be implemented by heightening the awareness of public health-care sectors on the spread and impacts of the disease on human health and by providing information for infected populations to use an effort to minimize the risk factors of the disease through intensifying successive education (or training) on the disease. However, to illustrate the effectiveness and cost-effectiveness of these intervention efforts an optimal control intervention should be studied.



---

*Source: 1806585-2022-09-17.xml*