# Dynamics ofMycobacterium and bovine tuberculosis in a Human-Buffalo Population

Authors: A. S. Hassan; S. M. Garba; A. B. Gumel; J. M.-S. Lubuma
Journal: Computational and Mathematical Methods in Medicine
 (2014)

---

## Abstract

A new model for the transmission dynamics ofMycobacterium tuberculosis and bovine tuberculosis in a community, consisting of humans and African buffalos, is presented. The buffalo-only component of the model exhibits the phenomenon of backward bifurcation, which arises due to the reinfection of exposed and recovered buffalos, when the associated reproduction number is less than unity. This model has a unique endemic equilibrium, which is globally asymptotically stable for a special case, when the reproduction number exceeds unity. Uncertainty and sensitivity analyses, using data relevant to the dynamics of the two diseases in the Kruger National Park, show that the distribution of the associated reproduction number is less than unity (hence, the diseases would not persist in the community). Crucial parameters that influence the dynamics of the two diseases are also identified. Both the buffalo-only and the buffalo-human model exhibit the same qualitative dynamics with respect to the local and global asymptotic stability of their respective disease-free equilibrium, as well as with respect to the backward bifurcation phenomenon. Numerical simulations of the buffalo-human model show that the cumulative number ofMycobacterium tuberculosis cases in humans (buffalos) decreases with increasing number of bovine tuberculosis infections in humans (buffalo).


---

## Body

## 1. Introduction

Mycobacterium tuberculosis (MTB) and bovine tuberculosis (BTB) are chronic bacterial diseases, classified amongst the closely related species that form theM. tuberculosis complex (MTBC) [1]. The human MTB is caused bytubercle bacillus (M. tuberculosis), while BTB is caused bybovine bacillus (M. bovis) [2]. MTB and BTB affect a wide range of hosts, including domestic livestock (such as cattle, goats, sheep, deer, and bison), wildlife (such as badgers, deer, bison, and African buffalo) which can either be reservoir or spillover, and humans [3].MTB remains a major global health problem affecting millions of people each year [4]. It is ranked second to human immunodeficiency virus (HIV) among the leading causes of death worldwide [4]. For instance, in the year 2012, there were 8.6 million new MTB cases and 1.3 million MTB deaths globally [4]. Similarly, BTB remains a serious problem for animal and human health in many developing countries [5]. Its widespread distribution has drastic negative socioeconomic impact, affecting public health, international trade, tourism, animal mortality, and milk production [6]. For example, in Argentina, the annual loss due to BTB is estimated to be US$ 63 million [7]. A benefit/cost analyses of BTB eradication in the United States showed an actual cost of US$ 538 million between 1917 and 1992 (current programs cost approximately US$ 3.5–4.0 millionper year [5]).The African buffalo transmits BTB to humans,via aerosol or oral (as a result of consuming raw unpasteurized milk) [1]. Furthermore, BTB can be transmitted from human to human by direct contact [1]. As in cattle, the main sources of BTB transmission in buffalo are direct contact, aerosol, oral, through a bite, or contamination of a skin wound [3] (other means of transmission, such as vertical and pseudovertical [8], also occur). Similarly, MTB can be transmitted from human to human, or from human to buffalo,via coughing or sneezing [1]. In humans, MTB is regarded to be an airborne disease [9]. It typically affects the lungs (pulmonary TB), but can affect other parts of the body also (extrapulmonary TB) [3]. Common signs and symptoms of MTB include coughing, chest pain, fever, weakness, and weight loss. The incubation period of MTB is approximately 2 to 12 weeks. African buffalos infected with BTB show clinical signs only when the disease has reached an advanced stage (the clinical signs of BTB in buffalo at such stage include coughing, debilitation, poor body condition or emaciation, and lagging when chased by helicopter [3, 8]). The incubation period for BTB is approximately 9 months to a year, and infections can remain dormant for years (and reactivate during periods of stress or in old age) [6].BTB is typically controlled using isolation or quarantine of infected herds, test-and-slaughter policy, and pasteurization of milk [10]. In South Africa’s Kruger National Park (KNP), other control measures, such as culling, vaccination, and their combination, are also used [10] (a demographic map of KNP and a herd of African Buffalos [11] are shown in Figure 1). Similarly, MTB in humans is controlledvia standard six-month course of four antimicrobial drugs [12–15]. The World Health Organization embarked on numerous global initiatives, such as “Stop TB Partnership,” “International Standards of Tuberculosis Care and Patient’s Care,” and the “Global Plan to Stop TB,” with the hope of minimizing the burden of TB worldwide [12].African buffalos and demographic map of Kruger National Park [11].
(a)
(b)Several mathematical models have been developed and used to gain insight into the transmission dynamics of BTB or MTB in populations (see, for instance, [8, 10, 12, 13, 16–19] and some of the references therein). However, none of these studies incorporate humans in the transmission dynamics of BTB. The purpose of the current study is to design, and analyse, a new realistic model (which extends some of the aforementioned studies in the literature) for BTB-MTB transmission dynamics. The objective is to gain insight into the qualitative dynamics of the two diseases in a buffalo-human population.The paper is organized as follows. The new model for BTB and MTB transmission dynamics in a community consisting of human and buffalo is presented in Section2. The buffalo-only model is rigorously analysed in Section 3, and the full buffalo-human model is analysed in Section 4. Numerical simulations are also reported.


## 2. Model Formulation

The model to be designed is based on the transmission dynamics of MTB and BTB in a population consisting of humans and African buffalos. The total human population at timet, denoted by NH(t), is subdivided into seven mutually exclusive compartments of susceptible humans (SH(t)), exposed humans (who have been infected with MTB but have not yet shown clinical symptoms of the disease) (EH1(t)), exposed humans with BTB (EH2(t)), infected humans with clinical symptoms of MTB (IH1(t)), infected humans with clinical symptoms of BTB (IH2(t)), and humans who recovered from MTB (RH1(t)) or BTB (RH2(t)), so that
(1)NH(t)=SH(t)+EH1(t)+EH2(t)+IH1(t)+IH2(t)+RH1(t)+RH2(t).
Similarly, the total buffalo population (in the herd) at time t, denoted by NB(t), is split into susceptible (SB(t)), early-exposed with BTB (EB1(t)), early-exposed with MTB (EM1(t)), advanced-exposed with BTB (EB2(t)), advanced-exposed with MTB (EM2(t)), infected with clinical symptoms of BTB (IBB(t)), infected with clinical symptoms of MTB (IMB(t)), and recovered from BTB (RBB(t)) or MTB (RMB(t)), so that
(2)NB(t)=SB(t)+EB1(t)+EM1(t)+EB2(t)+EM2(t)+IBB(t)+IMB(t)+RBB(t)+RMB(t).
The susceptible human population (SH(t)) is increased by the recruitment of people (either by birth or immigration) into the human-buffalo poulation (at a rate ΠH). The population is decreased by infection with MTB (at a rate λH) or BTB (at a rate λB1), where
(3)λH=βHNH(ηH1EH1+IH1),λB1=λHB+θMMλB,
with
(4)λHB=βHNH(ηH2EH2+IH2),λB=βBNB(ηB1EB1+ηB2EB2+IBB).
In (3) and (4), βH and βB represent the effective contact rates (i.e., contacts capable of leading to MTB or BTB infection), respectively. Furthermore, 0≤ηH1<1 and 0≤ηH2<1 are modification parameters accounting for the assumed reduction in infectiousness of exposed humans, in comparison to infected humans with clinical symptoms of MTB or BTB, respectively. Similarly, 0≤ηB1<1 and 0≤ηB2<1 are modification parameters accounting for the assumed reduction in infectiousness of exposed buffalos, in comparison to infected buffalos with clinical symptoms of BTB. The modification parameter 0≤θMM<1 accounts for the assumed reduced likelihood of susceptible humans acquiring BTB infection, in comparison to susceptible buffalos acquiring BTB infection. Natural death is assumed to occur in all human compartments at a rate μH. Thus, the rate of change of the susceptible human population is given by
(5)dSHdt=ΠH-(λH+λB1+μH)SH.
The population of exposed humans with MTB (EH1(t)) is generated by the infection of susceptible humans with MTB (at the rate λH) and is decreased by the development of clinical symptoms of MTB (at a rate σ1), exogenous reinfection (at a rate θH1λH, where 0≤θH1<1 accounts for the assumption that reinfection of exposed humans with MTB occurs at a rate lower than primary infection of susceptible humans with MTB) and natural death, so that
(6)dEH1dt=λHSH-(σ1+θH1λH+μH)EH1.
Similarly, the population of exposed humans with BTB (EH2(t)) is increased by the infection of susceptible humans with BTB (at the rate λB1) and is reduced by the development of clinical symptoms of BTB (at a rate σ2), exogenous reinfection (at a rate θH2λB1, with 0≤θH2<1 similarly defined as θH1) and natural death. Thus,
(7)dEH2dt=λB1SH-(σ2+θH2λB+μH)EH2.
The population of humans with clinical symptoms of MTB (IH1(t)) increases following the development of clinical symptoms of MTB by exposed humans (at the rate σ1) and exogenous reinfection of exposed and recovered humans (at the rates θH1λH and θRHλH, resp., with 0≤θRH<1). This population is decreased by recovery (at a rate γ1), natural death, and MTB-induced death (at a rate δH1), so that
(8)dIH1dt=σ1EH1+(θH1EH1+θRHRH1)λH-(γ1+μH+δH1)IH1.
The population of infected humans with clinical symptoms of BTB (IH2(t)) is generated by the development of clinical symptoms of BTB by exposed humans (at the rate σ2) and reinfection of exposed and recovered humans (at the rates θH2λB1 and θRBλB1, resp., with 0≤θRB<1). This population is decreased by recovery (at a rate γ2), natural death, and BTB-induced death (at a rate δH2). This gives
(9)dIH2dt=σ2EH2+(θH2EH2+θRBRH2)λB1-(γ2+μH+δH2)IH2.
The population of humans who recovered from MTB (RH1(t)) is generated by the recovery of humans with clinical symptoms of MTB (at the rate γ1). It is decreased by exogenous reinfection (at the rate θRHλH) and natural death. Hence,
(10)dRH1dt=γ1IH1-(θRHλH+μH)RH1.
It should be mentioned that, since MTB-infected humans do not completely eliminate the bacteria from their body (usually the bacteria hide in the bone marrow), “recovery” in this case implies (or represents) a long period of latency (which could last for a lifetime) [19, 30].Similarly, the population of humans who recovered from BTB(RH2(t)) is generated by the recovery of humans with clinical symptoms of BTB (at the rate γ2) and is decreased by reinfection (at the rate θRBλB1) and natural death, so that
(11)dRH2dt=γ2IH2-(θRBλB1+μH)RH2.
The population of susceptible buffalos (SB(t)) is generated by the recruitment of buffalos (either by birth or restocking from other herds) at a rate ΠB. It is assumed that all recruited buffalos are susceptible. The population of susceptible buffalos is decreased by acquisition of BTB infection (following effective contact with a human or buffalo infected with BTB), at the rate λB (where, λB=θBBλHB+λBB, with the modification parameters 0≤θBB<1 accounting for the expected reduced likelihood of humans transmitting of BTB to buffalo, in relation to BTB transmission from a human to another human) or MTB (following effective contact with a human infected with MTB) at a reduced rate θHHλH (where 0≤θHH<1 is a modification parameter accounting for the assumed reduction in the transmissibility of MTB from humans to buffalos, in comparison to MTB transmission from humans to humans), and by natural death (at a rate μB, buffalos in each epidemiological compartment suffer natural death at this rate). Thus,
(12)dSBdt=ΠB-(λB+θHHλH+μB)SB.
An important feature of BTB transmission within the buffalo population is that an infected buffalo could be in early or advanced stage of infection. This is owing to the fact that the clinical symptoms of BTB usually take months to manifest in buffalos [6]. Thus, BTB infections can remain dormant for years and reactivate during periods of stress or in old age [6]. These (early and advanced-exposed stage) features are incorporated in the model being develop. The population of buffalos early-exposed to BTB (EB1(t)) is increased by the infection of susceptible buffalos with BTB (at the rate λB). This population is decreased by exogenous reinfection with BTB (at a rate θEBλB, with 0≤θEB<1), progression to the advanced-exposed class (at a rate κ1), and natural death. This gives
(13)dEB1dt=λBSB-(θEBλB+κ1+μB)EB1.The population of buffalos early-exposed to MTB is increased by the infection of susceptible buffalos with MTB (at the rateθHHλH, where 0≤θHH<1 is as defined above). The population is decreased by exogenous reinfection (at a rate θEBλH), progression to the advanced-exposed MTB class (at a rate κ2), and natural death. This gives
(14)dEM1dt=θHHλHSB-(θEBλH+κ2+μB)EM1.The population of buffalos at advanced-exposed BTB class(EB2(t)) is increased by the progression of buffalos in the early-exposed BTB class (at the rate κ1). It is decreased by exogenous reinfection (at a rate θEBλB), development of clinical symptoms of BTB (at a rate σB2), and natural death, so that
(15)dEB2dt=κ1EB1-(θEBλB+σB2+μB)EB2.
Similarly, the population of buffalos at advanced-exposed MTB class (EM2(t)) is generated by the progression of buffalos in the early-exposed MTB class (at the rate κ2). It is decreased by exogenous reinfection (at a rate θEBλH), development of clinical symptoms of MTB (at a rate σM2), and natural death. Hence,
(16)dEM2dt=κ2EM1-(θEBλH+σM2+μB)EM2.
The population of buffalos with clinical symptoms of BTB (IBB(t)) is increased by the development of clinical symptoms of exposed buffalos with BTB (at the rate σB2) and by the exogenous reinfection of exposed and recovered buffalos (at the rates θEBλB and θRBλB, resp.). It is decreased by recovery (at a rate γB1), natural death, and BTB-induced mortality (at a rate δB). Thus,
(17)dIBBdt=σB2EB2+(EB1+EB2)θEBλB+θRBλBRBB-(γB1+μB+δB)IBB.
The population of buffalos with clinical symptoms of MTB (IMB(t)) is increased by the development of clinical symptoms of exposed buffalos with MTB (at the rate σM2) and by the exogenous reinfection of exposed and recovered buffalos (at the rates θEBλH and θRBλH, resp.). It is decreased by recovery (at a rate γM1), natural death, and MTB-induced mortality (at a rate δM). Thus,
(18)dIMBdt=σM2EM2+(EM1+EM2)θEBλH+θRBλHRMB-(γM1+μB+δM)IMB.
The population of buffalos who recovered from BTB (RBB(t)) is increased following the recovery of buffalos with clinical symptoms of BTB (at the rate γB1). It is decreased by reinfection (at the rate θRBλB) and natural death, so that
(19)dRBBdt=γB1IBB-(θRBλB+μB)RBB.
Finally, the population of buffalos who recovered from MTB (RMB(t)) is generated by the recovery of buffalos with MTB (at the rate γM1) and is decreased following reinfection (at the rate θRBλH) and natural death. This gives
(20)dRMBdt=γM1IMB-(θRBλH+μB)RMB.
It is assumed that recovered buffalos and humans acquire permanent natural immunity against BTB or MTB infection so that recovered buffalos and humans do not return to their respective susceptible class (albeit buffalos and humans in recovered classes can acquire reinfection).Thus, based on the above assumptions and formulations, the model for the BTB-MTB transmission dynamics in a human-buffalo population is given by the following deterministic system of nonlinear differential equations (a flow diagram of the model is depicted in Figure2, and the associated variables and parameters are described in Tables 1 and 2, resp.):
(21)HumanComponentdSHdt=ΠH-(λH+λB1+μH)SH,dEH1dt=λHSH-(θH1λH+σ1+μH)EH1,dEH2dt=λB1SH-(θH2λB1+σ2+μH)EH2,dIH1dt=σ1EH1+(θH1EH1+θRHRH1)λHiiiiiiiiii-(γ1+μH+δH1)IH1,dIH2dt=σ2EH2+(θH2EH2+θRBRH2)λB1iiiiiiiiii-(γ2+μH+δH2)IH2,dRH1dt=γ1IH1-(θRHλH+μH)RH1,dRH2dt=γ2IH2-(θRBλB1+μH)RH2.BuffaloComponentdSBdt=ΠB-(θHHλH+λB+μB)SB,dEB1dt=λBSB-(θEBλB+κ1+μB)EB1,dEM1dt=θHHλHSB-(θEBλH+κ2+μB)EM1,dEB2dt=κ1EB1-(θEBλB+σB2+μB)EB2,dEM2dt=κ2EM1-(θEBλH+σM2+μB)EM2,dIBBdt=σB2EB2+(EB1+EB2)θEBλB+θRBλBRBBiiiiiiiiii-(γB1+μB+δB)IBB,dIMBdt=σM2EM2+(EM1+EM2)θEBλH+θRBλHRMBiiiiiiiiiii-(γM1+μB+δM)IMB,dRBBdt=γB1IBB-(θRBλB+μB)RBB,dRMBdt=γM1IMB-(θRBλH+μB)RMB.Table 1
Description of the variables of the BTB-MTB model (21).
Variable
Interpretation
S
H
Population of susceptible humans
E
H
1
Population of humans exposed to MTB
E
H
2
Population of humans exposed to BTB
I
H
1
Population of infected humans with clinical symptoms of MTB
I
H
2
Population of infected humans with clinical symptoms of BTB
R
H
1
Population of humans who recovered from MTB
R
H
2
Population of humans who recovered from BTB
S
B
Population of susceptible buffalos
E
B
1
Population of buffalos early-exposed to BTB
E
M
1
Population of buffalos early-exposed to MTB
E
B
2
Population of buffalos at advanced-exposed BTB stage
E
M
2
Population of buffalos at advanced-exposed MTB stage
I
B
B
Population of buffalos with clinical symptoms of BTB
I
M
B
Population of buffalos with clinical symptoms of MTB
R
B
B
Population of buffalos who recovered from BTB
R
M
B
Population of buffalos who recovered from MTBTable 2
Description of parameters of the BTB-MTB model (21).
Parameter
Interpretation
Π
H
Recruitment rate of humans
Π
B
Recruitment rate of buffalos
μ
H
Natural death rate of humans
μ
B
Natural death rate of buffalos
β
H
Transmission rate of MTB
β
B
Transmission rate of BTB
η
H
1
Modification parameter for the reduction in infectiousness of exposed humans in comparison to humans with clinical symptoms of MTB
η
H
2
Modification parameter for the reduction in infectiousness of exposed humans in comparison to humans with clinical symptoms of BTB
η
B
1
,
η
B
2
Modification parameters for the reduction in infectiousness of exposed buffalos in comparison to buffalos with clinical symptoms of BTB
θ
H
H
,
θ
B
B
Modification parameters for the reduction in transmissibility of MTB to buffalos in comparison to humans
θ
M
M
Modification parameters for the reduction in transmissibility of BTB to humans in comparison to buffalos
γ
i
(
i
=
1
,
2
)
Recovery rate of humans
γ
B
1
,
γ
M
1
Recovery rate of buffalos
σ
i
(
i
=
1
,
2
)
Progression rate fromEHi to IHi class
κ
1
Progression rate fromEB1 to EB2 class
κ
2
Progression rate fromEM1 to EM2 class
σ
B
2
Progression rate fromEB2 to IBB class
σ
M
2
Progression rate fromEM2 to IMB class
θ
H
i
(
i
=
1
,
2
)
Exogenous reinfection rate for humans in theEHi class
θ
R
B
,
θ
R
H
Exogenous reinfection rate for recovered humans
θ
E
B
Exogenous reinfection rate for buffalos in the exposed and recovered classes, respectively
δ
H
1
,
δ
H
2
Disease-induced death rate for humans
δ
B
,
δ
M
Disease-induced death rate for buffalosFigure 2
Schematic diagram of the BTB-MTB model (21).The model (21) is, to the authors’ knowledge, the first to incorporate humans and MTB dynamics in the transmission dynamics of BTB in a human-buffalo community. Furthermore, it extends numerous models for BTB transmission in the literature, such as those in [8, 12, 13, 16–19], by,inter alia,(i)
including the dynamics of early and advanced-exposed buffalos (exposed buffalo classes were not considered in the models in [8, 12, 16–18]),(ii)
allowing for BTB and MTB transmission by exposed buffalos and humans (this was not considered in [8, 12, 16–19]),(iii)
including the dynamics of humans (this was not considered in [8, 13, 18, 19]),(iv)
allowing for the reinfection of exposed and recovered buffalos and humans (this was not considered in [8, 12, 13, 16, 18]),(v)
allowing for the transmission of both BTB and MTB in both the buffalo and human populations (this was not considered in the models in [8, 12, 13, 16–18]).The model (21) will now be rigorously analysed to gain insight into its dynamical features. Before doing so, it is instructive, however, to consider the dynamics within the buffalo population only as below.


## 3. Analysis of Buffalo-Only Model

Consider the model (21) in the absence of humans (buffalo-only model), obtained by setting the human components to zero (i.e., setting SH=EH1=EH2=IH1=IH2=RH1=RH2=λH=θHH=0 in (21)), given by
(22)dSBdt=ΠB-(λB+μB)SB,dEB1dt=λBSB-(θEBλB+κ1+μB)EB1,dEM1dt=-(κ2+μB)EM1,dEB2dt=κ1EB1-(θEBλB+σB2+μB)EB2,dEM2dt=κ2EM1-(σM2+μB)EM2,dIBBdt=σB2EB2+(EB1+EB2)θEBλB+θRBλBRBBIIIIIIII-(γB1+μB+δB)IBB,dIMBdt=σM2EM2-(γM1+μB+δM)IMB,dRBBdt=γB1IBB-(θRBλB+μB)RBB,dRMBdt=γM1IMB-μBRMB,
where, now,
(23)λB=βBNB(ηB1EB1+ηB2EB2+IBB).
The buffalo-only model (22) is fitted using data obtained from South Africa’s Kruger National Park [29], as shown in Figure 3 (from which it is evident that the model mimics the data reasonably well).Figure 3
Data fit of the simulation of the buffalo-only model (22), using data obtained from South Africa’s Kruger National Park (Table 4) [29]. Parameter values used are as given in Table 3.It is worth stating that since there are no humans in the dynamics of the buffalo-only model (22), MTB is not transmitted to susceptible buffalos. Furthermore, it is clear from the third equation in (22) that
(24)EM1(t)⟶0ast⟶∞.
Substituting (24) in the fifth equation in (22) shows that
(25)EM2(t)⟶0ast⟶∞.
Similarly, by substituting (EM1,EM2)=(0,0) into the equations for IMB and RMB in (22), it follows that
(26)(IMB(t),RMB(t))⟶(0,0)ast⟶∞.
Thus, the buffalo-only model reduces to the following (limited) model at steady-state:
(27)dSBdt=ΠB-(λB+μB)SB,dEB1dt=λBSB-(θEBλB+κ1+μB)EB1,dEB2dt=κ1EB1-(θEBλB+σB2+μB)EB2,dIBBdt=σB2EB2+(EB1+EB2)θEBλB+θRBλBRBBIIIIIIIII-(γB1+μB+δB)IBB,dRBBdt=γB1IBB-(θRBλB+μB)RBB.Lemma 1.
The following biologically feasible region of the buffalo-only model (27)
(28)Γ={ΠBμB(SB,EB1,EB2,IBB,RBB)∈R+5:SB+EB1II+EB2+IBB+RBB≤ΠBμB}
is positively invariant and attracting.Proof.
Adding the equations in the buffalo-only model system (27) gives
(29)dNB(t)dt=ΠB-μBNB(t)-δBIBB(t),
so that
(30)dNB(t)dt≤ΠB-μBNB(t).
It follows from (30) and the Gronwall inequality that
(31)NB(t)≤NB(0)e-μB(t)+ΠBμB[1-e-μB(t)].
In particular, NB(t)≤ΠB/μB if NB(0)≤ΠB/μB. Thus, Γ is positively invariant. Hence, it is sufficient to consider the dynamics of the buffalo-only model (27) in Γ (where the model can be considered to be epidemiologically and mathematically well-posed [31]).Theorem 2.
Let the initial dataSB(0)>0, EB1(0)>0, EB2(0)>0, IBB(0)>0, RBB(0)>0. Then, the solutions SB(t), EB1(t), EB2(t), IBB(t), and RBB(t) of the buffalo-only model (27) are positive for all t≥0.Proof.
It is clear from the first equation of the buffalo-only model (27) that
(32)dSBdt≥-(λB+μB)SB,
so that
(33)SB(t)≥SB(0)exp⁡[-∫0t‍(λB+μB)du]>0,∀t>0.
Using similar approach, it can be shown that EB1(t)>0, EB2(t)>0, IBB(t)>0, and RBB(t)>0, for all t>0.


### 3.1. Asymptotic Stability of Disease-Free Equilibrium (DFE)

#### 3.1.1.  Local Asymptotic Stability

The DFE of the buffalo-only model (27) is given by
(34)E0=(SB*,EB1*,EB2*,IBB*,RBB*)=(ΠBμB,0,0,0,0).
The linear stability of E0 can be established using the next generation operator method on the system (22) [32, 33]. The matrices F (for the new infection terms) and V (of the transition terms) associated with the system (27) are given, respectively, by
(35)F=[βBηB1βBηB2βB000000],V=[K100-κ1K300-σB2K5],
where K1=κ1+μB, K3=σB2+μB, and K5=γB1+μB+δB. It follows that thebasic reproduction number of the buffalo-only model (27), denoted by R0, is given by
(36)R0=βB[ηB1K3K5+κ1(ηB2K5+σB2)]K1K3K5.
Hence, using Theorem 2 of [33], the following result is established.Lemma 3.
The DFE,E0, of the buffalo-only model (27), is locally asymptotically stable (LAS) if R0<1 and unstable if R0>1.The threshold quantity,R0, represents the average number of secondary cases of BTB in the buffalo population that one BTB-infected buffalo can generate if introduced into a completely susceptible buffalo population [31, 34, 35].


#### 3.1.2. Interpretation ofR0

The threshold quantity,R0, can be interpreted as follows. It is worth recalling, first of all, that susceptible buffalos can acquire BTB infection following effective contact with either early-exposed buffalo with BTB (EB1(t)), advanced-exposed buffalo with BTB (EB2(t)), or infected buffalo with clinical symptoms of BTB (IBB(t)). It follows that the number of BTB infections generated by an early-exposed buffalo (near the DFE) is given by the product of the infection rate of an early-exposed buffalo (βBηB1/NB*) and the average duration of stay in the early-exposed class (1/K1). Thus, the average number of BTB infections generated by early-exposed buffalos is given by
(37)βBηB1SB*K1NB*.Similarly, the number of BTB infections generated by an advanced-exposed buffalo (near the DFE) is given by the product of the infection rate of advanced-exposed buffalos (βBηB2/NB*), the probability that early-exposed buffalo survived the early-exposed class and move to the advanced-exposed class (κ1/K1), and the average duration of stay in the advanced-exposed class (1/K3). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by
(38)βBηB2κ1SB*K1K3NB*.Furthermore, the number of BTB infections generated by an infected buffalo with clinical symptoms of BTB (near the DFE) is given by the product of the infection rate of buffalos with clinical symptoms of BTB (βB/NB*), the probability that an advanced-exposed buffalo survived the advanced-exposed class and move to the symptomatic class IBB (κ1σB2/K1K3), and the average duration of stay in the symptomatic class IBB (1/K5). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by
(39)βBκ1σB2SB*K1K3K5NB*.The sum of the terms in (37), (38), and (39) gives R0. That is, the average number of new infections generated by infected buffalos (early-exposed, advanced-exposed, or symptomatic) is given by (noting that SB*=ΠB/μB and NB*=ΠB/μB)
(40)R0=βB[ηB1K3K5+κ1(ηB2K5+σB2)]K1K3K5.The epidemiological implication of Lemma3 is that BTB can be effectively controlled in (or eliminated from) the buffalo population (herd) if the initial sizes of the state variables of the buffalo-only model (27) are in the basin of attraction of the DFE (E0). It is worth mentioning, however, that TB models with exogenous reinfection are often shown to exhibit the phenomenon of backward bifurcation (where the stable DFE coexists with a stable endemic equilibrium when R0<1 [12, 17, 25, 36]). The epidemiological implication of this phenomenon is that the classical requirement of R0<1 is, although necessary, no longer sufficient for diseases elimination [12, 36]. Thus, the presence of backward bifurcation in the transmission dynamics of a disease makes its effective control in a population more difficult. Hence, it is instructive to explore the possibility of such phenomenon in the buffalo-only model (22). This is investigated below.Theorem 4.
The buffalo-only model (22) undergoes backward bifurcation at R0=1 whenever the bifurcation coefficient, a, given by (A.9) (in Appendix A), is positive.The proof of Theorem4, based on using centre manifold theory [17, 33], is given in Appendix A. It should be noted that, in the absence of reinfection of exposed and recovered buffalos (i.e., the case of the model (22) with θEB=θRB=0), the backward bifurcation coefficient, a, given by (A.9) in Appendix A, reduces to (it should be noted, from Appendix A, that A1>0, A2>0,βB*>0, and, from Theorem 2, that all parameters of the buffalo-only model (22) are nonnegative):
(41)-2βB*μB2ΠB{1+μB(K5+σB1)+σB2γB1A1μBσB2iiiiiiiiiiiiiiiii+γM1σM2K2+K2μBK6+σB2A2K2μBσM2},<0,
where K2=κ2+μB, K4=σM2+μB and K6=γM1+μB+δM. Since the bifurcation coefficient, a, is automatically negative, it follows from the analyses in Appendix A, and Theorem 4.1 of [17], that the buffalo-only model (22) does not undergo backward bifurcation in the absence of reinfection (this result is consistent with that in [12, 17, 36], on the transmission dynamicsMycobacterium tuberculosis in human populations). This result is summarized below.Lemma 5.
The buffalo-only model (22) does not undergo backward bifurcation at R0=1 in the absence of the reinfection of exposed and recovered buffalos (θEB=θRB=0).Hence, this study shows that the reinfection of exposed and recovered buffalos causes the phenomenon of backward bifurcation in the transmission dynamics of BTB and MTB in a buffalo-only population. To further confirm the absence of backward bifurcation in this case, a global asymptotic stability result is established for the DFE(E0) of the buffalo-only model (27) in the absence of reinfection (i.e., θEB=θRB=0) below.


### 3.2. Global Asymptotic Stability of the DFE

Consider the buffalo-only model (27) in the absence of the reinfection of exposed (θEB=0) and recovered (θRB=0) buffalos.Theorem 6.
The DFE,E0, of the buffalo-only model (27) with θEB=θRB=0 is globally asymptotically stable (GAS) in Γ if R0≤1.Proof.
Consider the buffalo-only model (27) in the absence of reinfection (θEB=θRB=0). Furthermore, let R0≤1. Consider the linear Lyapunov function F=a0EB1+a1EB2+a2IBB, where
(42)a0=R0,a1=βB(ηB2K5+σB2)K3K5,a2=βBK5,
with Lyapunov derivative given by (where a dot represents differentiation with respect to time t)
(43)F˙=a0E˙B1+a1E˙B2+a2I˙BB,=a0[βBNB(ηB1EB1+ηB2EB2+IBB)SB-K1EB1]+a1(κ1EB1-K3EB2)iii+a2(σB2EB2-K5IBB),=(a0βBηB1SBNB-a0K1+a1κ1)EB1+(a0βBηB2SBNB-a1K3+a2σB2)EB2+(a0βBSBNB-a2K5)IBB,≤βB(ηB1EB1+ηB2EB2+IBB)(R0-1)iiiiiiiiiiiiiiiiiisinceSB(t)≤NB(t)∀tinΓ,≤0ifR0≤1.
Since all the parameters and variables of model (27) are nonnegative (Theorem 2), it follows that F˙≤0 for R0≤1 with F˙=0 if and only if EB1=EB2=IBB=0. Thus, it follows, by LaSalle’s Invariance Principle [37], that
(44)(EB1(t),EB2(t),IBB(t))⟶(0,0,0)ast⟶∞.
Since lim⁡t→∞⁡sup⁡IBB(t)=0 (from (44)), it follows that, for sufficiently small ϖ*>0, there exists a constant M>0, such that lim⁡t→∞⁡sup⁡IBB(t)≤ϖ* for all t>M. Hence, it follows from the fifth equation of the buffalo-only model (27) that, for t>M,R˙BB≤γB1ϖ*-μBRBB. Thus, by comparison theorem [38], RBB∞=lim⁡t→∞⁡sup⁡RBB≤γB1ϖ*/μB, so that, by letting, ϖ*→0,
(45)RBB∞=lim⁡t→∞sup⁡RBB(t)≤0.
Similarly, it can be shown that
(46)RBB∞=lim⁡t→∞inf⁡RBB(t)≥0.
Thus, it follows from (45) and (46), that RBB∞≥0≥RBB∞. Hence,
(47)lim⁡t→∞RBB(t)=0.
Furthermore, substituting (44) in the first equation of (27) shows that
(48)SB(t)⟶ΠBμBast⟶∞.
Thus, by combining equations (44), (47), and (48), it follows that every solution of the equations of the buffalo-only model (27), with θEB=θRB=0 and initial conditions in Γ, approaches E0 as t→∞ (whenever R0≤1).Theorem6 shows that, in the absence of the reinfection of exposed and recovered buffalos (i.e., θEB=θRB=0), BTB can be eliminated from the buffalo-only population if the reproduction number of the model (R0) can be brought to (and maintained at) a value less than unity. Figure 4(a) depicts the solution profiles of the buffalo-only model (27), generated using various initial conditions, showing convergence to the DFE E0 when R0<1 (in line with Theorem 6).Simulations of the buffalo-only model (27), showing the total number of infected buffalos with clinical symptoms of BTB (IBB(t)(t)) at time t as a function of time. Parameter values used are as given in Table 3 with (a) βB=0.00733 (so that, R0=0.7036<1) and (b) βB= 0.0733, δB=0 (so that, R~0=8.6050>1).
(a)
(b)


### 3.3. Existence of Endemic Equilibria: Special Case

In this section, the existence of nontrivial (endemic) equilibria (where the components of the infected variables of the model are nonzero) of the buffalo-only model (27) is explored for the special case without reinfection (i.e., θEB=θRB=0). Solving the equations of the buffalo-only model (27) at steady-state gives the following general form of the endemic equilibrium (denoted by E1):
(49)E1=(SB**,EB1**,EB2**,IBB**,RBB**),
where
(50)SB**=ΠBλB**+μB,EB1**=λB**ΠBK1(λB**+μB),EB2**=κ1λB**ΠBK1K3(λB**+μB),IBB**=σB2κ1λB**ΠBK1K3K5(λB**+μB),RBB**=γB1σB2κ1λB**ΠBK1K3K5μB(λB**+μB),
with the force of infection at steady-state (λB**) given by
(51)λB**=βBNB**(ηB1EB1**+ηB2EB2**+IBB**).Using (50) in the expression for λB** in (51) shows that the nonzero equilibrium of the model (22) satisfies the linear equation:
(52)b1λB**+b2=0,
where b1=K5μB(K3+κ1)+σB2κ1(μB+γB2) and b2=K1K3K5μB(1-R0). Clearly, the coefficient b1 is always positive, and b2 is positive (negative) if R0 is less than (greater than) unity, respectively. Thus, the linear system (52) has a unique positive solution, given by λB**=-b2/b1, whenever R0>1. Further, the force of infection for buffalos (λB**) is negative whenever R0<1 (which is biologically meaningless). Hence, the buffalo-only model (27) has no positive equilibrium in this case. These results are summarized below.Theorem 7.
The buffalo-only model (27), with θEB=θRB=0, has a unique endemic equilibrium, E1, whenever R0>1, and no endemic equilibrium otherwise.


#### 3.3.1. Global Asymptotic Stability of Endemic Equilibrium

The global asymptotic stability of the unique endemic equilibrium (E1) of the buffalo-only model is explored for the special case without reinfection (θEB=θRB=0) and BTB-induced death in buffalos (δB=0). It is convenient to define
(53)Γ1={(SB,EB1,EB2,IBB,RBB)∈Γ:EB1=EB2iiiii=IBB=RBB=0},
the stable manifold of the DFE (E0) of the buffalo-only model (27).Theorem 8.
The unique endemic equilibrium (E1) of the buffalo-only model (27), with θEB=θRB=δB=0, is GAS in Γ∖Γ1 if R~0=R0|δB=0>1.Proof.
Consider the buffalo-only model (27) with θEB=θRB=δB=0. For this case, it follows from Theorem 7 that the buffalo-only model (27) has a unique endemic equilibrium whenever R~0>1. Furthermore, setting δB=0 in model (27) shows that NB(t)→ΠB/μB as t→∞. Consider the following nonlinear Lyapunov function (of Goh-Volterra type) for the subsystem of model (27) involving the state variables SB, EB1, EB2, and IBB (noting that NB(t) is now replaced by its limiting value ΠB/μB):
(54)F=SB-SB**-SB**ln⁡(SBSB**)+EB1-EB1**-EB1**ln⁡(EB1EB1**)+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[EB2-EB2**-EB2**ln⁡(EB2EB2**)]+β~BSB**IBB**σB2EB2**[IBB-IBB**-IBB**ln⁡(IBBIBB**)],
where β~B=μBβB/ΠB. The Lyapunov derivative of F is given by
(55)F˙=S˙B-SB**SBS˙B+E˙B1-EB1**EB1E˙B1+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)(E˙B2-EB2**EB2E˙B2)+β~BSB**IBB**σB2EB2**(I˙BB-IBB**IBBI˙BB),=ΠB-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB-SB**SB[ΠB-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB]+β~B(ηB1EB1+ηB2EB2+IBB)SB-K1EB1-EB1**EB1[β~B(ηB1EB1+ηB2EB2+IBB)SB-K1EB1]+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[κ1EB1-K3EB2-EB2**EB2(κ1EB1-K3EB2)]+β~BSB**IBB**σB2EB2**[σB2EB2-K5IBB-IBB**IBBDDDDDDDD×(σB2EB2-K5IBB)σB2EB2-K5IBB-IBB**IBB].
Using the following steady-state relations (obtained from (27)),
(56)ΠB=β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**,κ1EB1**=K3EB2**,β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**=K1EB1**,σB2EB2**=K5IBB**,γB1IBB**=μBRBB**,
the Lyapunov derivative can be simplified to
(57)F˙=β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**-μBSB-SB**SB[β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**iiiiiiiiiiii-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB]-K1EB1-EB1**EB1[β~B(ηB1EB1+ηB2EB2+IBB)SBddddddddddddiiiii-K1EB1[β~B(ηB1EB1+ηB2EB2+IBB)SB]+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[κ1EB1-K3EB2-EB2**EB2(κ1EB1-K3EB2)]+β~BSB**IBB**σB2EB2**×[σB2EB2-K5IBB-IBB**IBB(σB2EB2-K5IBB)].
Thus,(58)F˙=μBSB**(2-SB**SB-SBSB**)+β~BηB1EB1**SB**(2-SB**SB-SBSB**)+β~BηB2EB2**SB**(3-SB**SB-EB1EB2**EB1**EB2-SBEB1**EB2SB**EB1EB2**)+β~BIBB**SB**×(4-SB**SB-EB1EB2**EB1**EB2-EB2IBB**EB2**IBB-SBEB1**IBBSB**EB1IBB**).
Finally, since the arithmetic mean exceeds the geometric mean, it follows then that
(59)μBSB**(2-SB**SB-SBSB**)≤0,β~BηB1EB1**SB**(2-SB**SB-SBSB**)≤0,β~BηB2EB2**SB**(3-SB**SB-EB1EB2**EB1**EB2-SBEB1**EB2SB**EB1EB2**)≤0,β~BIBB**SB**(4-SB**SB-EB1EB2**EB1**EB2-EB2IBB**EB2**IBB-SBEB1**IBBSB**EB1IBB**)≤0.
Furthermore, since all the model parameters are nonnegative, it follows that F˙≤0 for R~0>1. Thus, F˙ is a Lyapunov function for the subsystem of model (27) on Γ∖Γ1. Therefore, it follows, by LaSalle’s Invariance Principle [37], that
(60)lim⁡t→∞SB(t)=SB**,lim⁡t→∞EB1(t)=EB1**,lim⁡t→∞EB2(t)=EB2**,lim⁡t→∞IBB(t)=IBB**.
Since IBB(t)→IBB** as t→∞, it follows from the equation for dRBB/dt in (27) that RBB(t)→γB1IBB**/μB=RBB** as t→∞. The proof is concluded using similar arguments as in the proof of Theorem 6.The epidemiological implication of Theorem8 is that BTB will be endemic in the buffalo population if R~0>1 (and θEB=θRB=δB=0). Figure 4(b) depicts the solutions of model (27) for the case when R~0>1 and θEB=θRB=δB=0, showing convergence of the initial solutions to the unique endemic equilibrium (in line with Theorem 8).


### 3.4. Sensitivity and Uncertainty Analyses

In this section, sensitivity and uncertainty analyses will be carried out, using Latin hypercube sampling (LHS) and partial correlation coefficient (PRCC) [39–41], to assess the effect of uncertainty in the estimate of the parameter values used to simulate the buffalo-only model (on the simulation results obtained) and to determine the key parameters that drive the dynamics of the disease in the buffalo-human population. The ranges and baseline values of the parameters of the buffalo-only model, given in Table 3 with θRB=θEB=0 (i.e., in the absence of backward bifurcation), will be used in these analyses. Each parameter of the buffalo-only model (22) is assumed to obey a uniform distribution [42]. Following [42], a total of 1000 LHS runs (N=1000) are carried out. Furthermore, the following initial conditions (which are consistent with the dynamics of African buffalo in the Kruger National Park [29]): (SB(0), EB1(0), EM1(0), EB2(0), EM2(0), IBB(0), IMB(0), RBB(0), RMB(0))=(28000,100,100,20,20,10,10,100,100)) are used in the simulations.Table 3
Ranges and baseline values for parameters of the BTB-MTB model (21).
Parameter
Range (day−1)
Baseline value (day−1)
Reference
Π
H
[26, 80]
53
[20]
Π
B
[2, 4]
3
[20–22]
μ
H
[0.0000274, 0.0000549]
0.000047
[14, 23–25]
μ
B
(0.00009477, 0.00011583)
0.0001053
[22, 26]
β
H
[0.00011, 0.000959]
0.000535
[14, 17]
β
B
(0.006597, 0.008063)
0.00733
[10]
η
H
1
[0, 1)
0.5
Fitted
η
H
2
[0, 1)
0.5
Fitted
η
B
1
(0.4455, 0.5045)
0.45
Fitted
η
B
2
(0.495, 0.605)
0.55
Fitted
θ
B
B
[0, 1)
0.5
Fitted
θ
M
M
[0, 1]
0.5
Assumed
θ
H
H
[0, 1)
0.5
Assumed
γ
i
(
i
=
1
,
2
)
(0.0000823, 0.000823)
0.000453
[14, 23]
γ
B
1
(0.00774, 0.00946)
0.0086
[10]
γ
M
1
(0.13374, 0.160086)
0.1486
[12]
σ
i
(
i
=
1
,
2
)
(0.0000822, 0.00247)
0.001
[23, 27]
κ
1
(0.45, 0.55)
0.5
[27]
κ
2
(0.45, 0.55)
0.5
[27]
σ
B
2
(0.25, 0.35)
0.3
[27]
σ
M
2
(0.36, 0.44)
0.4
[27]
θ
H
i
(
i
=
1
,
2
)
[0, 0.1]
0.00271
[23]
θ
R
B
,
θ
R
H
[0.002439, 0.002981]
0.00271
[23]
θ
E
B
[0.002439, 0.002981]
0.00271
[23]
δ
H
1
,
δ
H
2
[0.000115, 0.000822]
0.0002
[14, 17, 23, 25]
δ
B
[0.0018, 0.0022]
0.002
[28]
δ
M
[0.0018, 0.0022]
0.002
[14, 23]Table 4
Number of symptomatic buffalos with BTB at Kruger National Park [29].
Year
Number of symptomatic buffalos [29]
2001
35
2002
135
2003
185
2004
238
2005
230Figure5 depicts a box plot of R0, as a function of the number of LHS runs carried out (N=1000), from which it is evident that the distribution of R0 lies in the range R0∈[0.25,0.75] (each box plot displays the upper and lower quartile ranges of R0. a horizontal line within the box is the median value, and values of R0 beyond the whiskers are outliers [43]). Thus, since the distribution of the reproduction number of the buffalo-only model is less than unity, it follows (from Theorem 6) that the BTB-MTB outbreaks (in the buffalo-human population) will die out with time (in other words, the disease will be effectively controlled). The PRCC values of the parameters of the buffalo-only model (22), using R0 as the response function, are depicted in Figure 6. It follows from Figure 6 that the top three parameters that most influences the value of R0 (hence the disease dynamics) are the BTB transmission rate (βB), the recovery rate of buffalos (γB1), and the BTB-induced mortality in buffalos (δB).Figure 5
Box plot ofR0 as a function of the number of LHS runs carried out for the buffalo-only model (22), using parameter values and ranges given in Table 3 with θEB=θRB=0.Figure 6
PRCC values of the parameters of the buffalo-only model (22), using R0 as the output function. Parameter values used are as given in Table 3.Similarly, Figure7 depicts the box plot of the buffalo-only model (22) using total number of symptomatic buffalos (IBB+IMB) as the response function. This figure shows a distribution of the number of symptomatic buffalos lying in the range [20–130]. Hence, this study shows that, using the parameter values and ranges relevant to BTB-MTB dynamics at the Kruger National Park, a BTB outbreak could cause no more than 130 confirmed cases (of BTB and MTB) in the park. The associated PRCC values (with the total number of symptomatic buffalos as the output) are depicted in Figure 8, from which it is evident that, in this scenario, the top three parameters (that most influences the output) are the buffalo recruitment rate (ΠB), the natural (μB), and the disease-induced (δB) death rate of buffalos. Hence, this study shows variability in the top-ranked PRCC values on the chosen response/output function.Figure 7
Box plot of the total number of symptomatic buffalos (IBB+IMB) as a function of the number of LHS runs for the buffalo-only model (22), using parameter values and ranges given in Table 3 with θEB=θRB=0.Figure 8
PRCC values of the parameters of the buffalo-only model (22), using total number of symptomatic buffalos (IBB+IMB) as the output function. Parameter values used are as given in Table 3.Having fully studied the dynamics of the bovine-only model (22), the full BTB-MTB model (22) will now be analysed.


## 3.1. Asymptotic Stability of Disease-Free Equilibrium (DFE)

### 3.1.1.  Local Asymptotic Stability

The DFE of the buffalo-only model (27) is given by
(34)E0=(SB*,EB1*,EB2*,IBB*,RBB*)=(ΠBμB,0,0,0,0).
The linear stability of E0 can be established using the next generation operator method on the system (22) [32, 33]. The matrices F (for the new infection terms) and V (of the transition terms) associated with the system (27) are given, respectively, by
(35)F=[βBηB1βBηB2βB000000],V=[K100-κ1K300-σB2K5],
where K1=κ1+μB, K3=σB2+μB, and K5=γB1+μB+δB. It follows that thebasic reproduction number of the buffalo-only model (27), denoted by R0, is given by
(36)R0=βB[ηB1K3K5+κ1(ηB2K5+σB2)]K1K3K5.
Hence, using Theorem 2 of [33], the following result is established.Lemma 3.
The DFE,E0, of the buffalo-only model (27), is locally asymptotically stable (LAS) if R0<1 and unstable if R0>1.The threshold quantity,R0, represents the average number of secondary cases of BTB in the buffalo population that one BTB-infected buffalo can generate if introduced into a completely susceptible buffalo population [31, 34, 35].


### 3.1.2. Interpretation ofR0

The threshold quantity,R0, can be interpreted as follows. It is worth recalling, first of all, that susceptible buffalos can acquire BTB infection following effective contact with either early-exposed buffalo with BTB (EB1(t)), advanced-exposed buffalo with BTB (EB2(t)), or infected buffalo with clinical symptoms of BTB (IBB(t)). It follows that the number of BTB infections generated by an early-exposed buffalo (near the DFE) is given by the product of the infection rate of an early-exposed buffalo (βBηB1/NB*) and the average duration of stay in the early-exposed class (1/K1). Thus, the average number of BTB infections generated by early-exposed buffalos is given by
(37)βBηB1SB*K1NB*.Similarly, the number of BTB infections generated by an advanced-exposed buffalo (near the DFE) is given by the product of the infection rate of advanced-exposed buffalos (βBηB2/NB*), the probability that early-exposed buffalo survived the early-exposed class and move to the advanced-exposed class (κ1/K1), and the average duration of stay in the advanced-exposed class (1/K3). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by
(38)βBηB2κ1SB*K1K3NB*.Furthermore, the number of BTB infections generated by an infected buffalo with clinical symptoms of BTB (near the DFE) is given by the product of the infection rate of buffalos with clinical symptoms of BTB (βB/NB*), the probability that an advanced-exposed buffalo survived the advanced-exposed class and move to the symptomatic class IBB (κ1σB2/K1K3), and the average duration of stay in the symptomatic class IBB (1/K5). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by
(39)βBκ1σB2SB*K1K3K5NB*.The sum of the terms in (37), (38), and (39) gives R0. That is, the average number of new infections generated by infected buffalos (early-exposed, advanced-exposed, or symptomatic) is given by (noting that SB*=ΠB/μB and NB*=ΠB/μB)
(40)R0=βB[ηB1K3K5+κ1(ηB2K5+σB2)]K1K3K5.The epidemiological implication of Lemma3 is that BTB can be effectively controlled in (or eliminated from) the buffalo population (herd) if the initial sizes of the state variables of the buffalo-only model (27) are in the basin of attraction of the DFE (E0). It is worth mentioning, however, that TB models with exogenous reinfection are often shown to exhibit the phenomenon of backward bifurcation (where the stable DFE coexists with a stable endemic equilibrium when R0<1 [12, 17, 25, 36]). The epidemiological implication of this phenomenon is that the classical requirement of R0<1 is, although necessary, no longer sufficient for diseases elimination [12, 36]. Thus, the presence of backward bifurcation in the transmission dynamics of a disease makes its effective control in a population more difficult. Hence, it is instructive to explore the possibility of such phenomenon in the buffalo-only model (22). This is investigated below.Theorem 4.
The buffalo-only model (22) undergoes backward bifurcation at R0=1 whenever the bifurcation coefficient, a, given by (A.9) (in Appendix A), is positive.The proof of Theorem4, based on using centre manifold theory [17, 33], is given in Appendix A. It should be noted that, in the absence of reinfection of exposed and recovered buffalos (i.e., the case of the model (22) with θEB=θRB=0), the backward bifurcation coefficient, a, given by (A.9) in Appendix A, reduces to (it should be noted, from Appendix A, that A1>0, A2>0,βB*>0, and, from Theorem 2, that all parameters of the buffalo-only model (22) are nonnegative):
(41)-2βB*μB2ΠB{1+μB(K5+σB1)+σB2γB1A1μBσB2iiiiiiiiiiiiiiiii+γM1σM2K2+K2μBK6+σB2A2K2μBσM2},<0,
where K2=κ2+μB, K4=σM2+μB and K6=γM1+μB+δM. Since the bifurcation coefficient, a, is automatically negative, it follows from the analyses in Appendix A, and Theorem 4.1 of [17], that the buffalo-only model (22) does not undergo backward bifurcation in the absence of reinfection (this result is consistent with that in [12, 17, 36], on the transmission dynamicsMycobacterium tuberculosis in human populations). This result is summarized below.Lemma 5.
The buffalo-only model (22) does not undergo backward bifurcation at R0=1 in the absence of the reinfection of exposed and recovered buffalos (θEB=θRB=0).Hence, this study shows that the reinfection of exposed and recovered buffalos causes the phenomenon of backward bifurcation in the transmission dynamics of BTB and MTB in a buffalo-only population. To further confirm the absence of backward bifurcation in this case, a global asymptotic stability result is established for the DFE(E0) of the buffalo-only model (27) in the absence of reinfection (i.e., θEB=θRB=0) below.


## 3.1.1.  Local Asymptotic Stability

The DFE of the buffalo-only model (27) is given by
(34)E0=(SB*,EB1*,EB2*,IBB*,RBB*)=(ΠBμB,0,0,0,0).
The linear stability of E0 can be established using the next generation operator method on the system (22) [32, 33]. The matrices F (for the new infection terms) and V (of the transition terms) associated with the system (27) are given, respectively, by
(35)F=[βBηB1βBηB2βB000000],V=[K100-κ1K300-σB2K5],
where K1=κ1+μB, K3=σB2+μB, and K5=γB1+μB+δB. It follows that thebasic reproduction number of the buffalo-only model (27), denoted by R0, is given by
(36)R0=βB[ηB1K3K5+κ1(ηB2K5+σB2)]K1K3K5.
Hence, using Theorem 2 of [33], the following result is established.Lemma 3.
The DFE,E0, of the buffalo-only model (27), is locally asymptotically stable (LAS) if R0<1 and unstable if R0>1.The threshold quantity,R0, represents the average number of secondary cases of BTB in the buffalo population that one BTB-infected buffalo can generate if introduced into a completely susceptible buffalo population [31, 34, 35].


## 3.1.2. Interpretation ofR0

The threshold quantity,R0, can be interpreted as follows. It is worth recalling, first of all, that susceptible buffalos can acquire BTB infection following effective contact with either early-exposed buffalo with BTB (EB1(t)), advanced-exposed buffalo with BTB (EB2(t)), or infected buffalo with clinical symptoms of BTB (IBB(t)). It follows that the number of BTB infections generated by an early-exposed buffalo (near the DFE) is given by the product of the infection rate of an early-exposed buffalo (βBηB1/NB*) and the average duration of stay in the early-exposed class (1/K1). Thus, the average number of BTB infections generated by early-exposed buffalos is given by
(37)βBηB1SB*K1NB*.Similarly, the number of BTB infections generated by an advanced-exposed buffalo (near the DFE) is given by the product of the infection rate of advanced-exposed buffalos (βBηB2/NB*), the probability that early-exposed buffalo survived the early-exposed class and move to the advanced-exposed class (κ1/K1), and the average duration of stay in the advanced-exposed class (1/K3). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by
(38)βBηB2κ1SB*K1K3NB*.Furthermore, the number of BTB infections generated by an infected buffalo with clinical symptoms of BTB (near the DFE) is given by the product of the infection rate of buffalos with clinical symptoms of BTB (βB/NB*), the probability that an advanced-exposed buffalo survived the advanced-exposed class and move to the symptomatic class IBB (κ1σB2/K1K3), and the average duration of stay in the symptomatic class IBB (1/K5). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by
(39)βBκ1σB2SB*K1K3K5NB*.The sum of the terms in (37), (38), and (39) gives R0. That is, the average number of new infections generated by infected buffalos (early-exposed, advanced-exposed, or symptomatic) is given by (noting that SB*=ΠB/μB and NB*=ΠB/μB)
(40)R0=βB[ηB1K3K5+κ1(ηB2K5+σB2)]K1K3K5.The epidemiological implication of Lemma3 is that BTB can be effectively controlled in (or eliminated from) the buffalo population (herd) if the initial sizes of the state variables of the buffalo-only model (27) are in the basin of attraction of the DFE (E0). It is worth mentioning, however, that TB models with exogenous reinfection are often shown to exhibit the phenomenon of backward bifurcation (where the stable DFE coexists with a stable endemic equilibrium when R0<1 [12, 17, 25, 36]). The epidemiological implication of this phenomenon is that the classical requirement of R0<1 is, although necessary, no longer sufficient for diseases elimination [12, 36]. Thus, the presence of backward bifurcation in the transmission dynamics of a disease makes its effective control in a population more difficult. Hence, it is instructive to explore the possibility of such phenomenon in the buffalo-only model (22). This is investigated below.Theorem 4.
The buffalo-only model (22) undergoes backward bifurcation at R0=1 whenever the bifurcation coefficient, a, given by (A.9) (in Appendix A), is positive.The proof of Theorem4, based on using centre manifold theory [17, 33], is given in Appendix A. It should be noted that, in the absence of reinfection of exposed and recovered buffalos (i.e., the case of the model (22) with θEB=θRB=0), the backward bifurcation coefficient, a, given by (A.9) in Appendix A, reduces to (it should be noted, from Appendix A, that A1>0, A2>0,βB*>0, and, from Theorem 2, that all parameters of the buffalo-only model (22) are nonnegative):
(41)-2βB*μB2ΠB{1+μB(K5+σB1)+σB2γB1A1μBσB2iiiiiiiiiiiiiiiii+γM1σM2K2+K2μBK6+σB2A2K2μBσM2},<0,
where K2=κ2+μB, K4=σM2+μB and K6=γM1+μB+δM. Since the bifurcation coefficient, a, is automatically negative, it follows from the analyses in Appendix A, and Theorem 4.1 of [17], that the buffalo-only model (22) does not undergo backward bifurcation in the absence of reinfection (this result is consistent with that in [12, 17, 36], on the transmission dynamicsMycobacterium tuberculosis in human populations). This result is summarized below.Lemma 5.
The buffalo-only model (22) does not undergo backward bifurcation at R0=1 in the absence of the reinfection of exposed and recovered buffalos (θEB=θRB=0).Hence, this study shows that the reinfection of exposed and recovered buffalos causes the phenomenon of backward bifurcation in the transmission dynamics of BTB and MTB in a buffalo-only population. To further confirm the absence of backward bifurcation in this case, a global asymptotic stability result is established for the DFE(E0) of the buffalo-only model (27) in the absence of reinfection (i.e., θEB=θRB=0) below.


## 3.2. Global Asymptotic Stability of the DFE

Consider the buffalo-only model (27) in the absence of the reinfection of exposed (θEB=0) and recovered (θRB=0) buffalos.Theorem 6.
The DFE,E0, of the buffalo-only model (27) with θEB=θRB=0 is globally asymptotically stable (GAS) in Γ if R0≤1.Proof.
Consider the buffalo-only model (27) in the absence of reinfection (θEB=θRB=0). Furthermore, let R0≤1. Consider the linear Lyapunov function F=a0EB1+a1EB2+a2IBB, where
(42)a0=R0,a1=βB(ηB2K5+σB2)K3K5,a2=βBK5,
with Lyapunov derivative given by (where a dot represents differentiation with respect to time t)
(43)F˙=a0E˙B1+a1E˙B2+a2I˙BB,=a0[βBNB(ηB1EB1+ηB2EB2+IBB)SB-K1EB1]+a1(κ1EB1-K3EB2)iii+a2(σB2EB2-K5IBB),=(a0βBηB1SBNB-a0K1+a1κ1)EB1+(a0βBηB2SBNB-a1K3+a2σB2)EB2+(a0βBSBNB-a2K5)IBB,≤βB(ηB1EB1+ηB2EB2+IBB)(R0-1)iiiiiiiiiiiiiiiiiisinceSB(t)≤NB(t)∀tinΓ,≤0ifR0≤1.
Since all the parameters and variables of model (27) are nonnegative (Theorem 2), it follows that F˙≤0 for R0≤1 with F˙=0 if and only if EB1=EB2=IBB=0. Thus, it follows, by LaSalle’s Invariance Principle [37], that
(44)(EB1(t),EB2(t),IBB(t))⟶(0,0,0)ast⟶∞.
Since lim⁡t→∞⁡sup⁡IBB(t)=0 (from (44)), it follows that, for sufficiently small ϖ*>0, there exists a constant M>0, such that lim⁡t→∞⁡sup⁡IBB(t)≤ϖ* for all t>M. Hence, it follows from the fifth equation of the buffalo-only model (27) that, for t>M,R˙BB≤γB1ϖ*-μBRBB. Thus, by comparison theorem [38], RBB∞=lim⁡t→∞⁡sup⁡RBB≤γB1ϖ*/μB, so that, by letting, ϖ*→0,
(45)RBB∞=lim⁡t→∞sup⁡RBB(t)≤0.
Similarly, it can be shown that
(46)RBB∞=lim⁡t→∞inf⁡RBB(t)≥0.
Thus, it follows from (45) and (46), that RBB∞≥0≥RBB∞. Hence,
(47)lim⁡t→∞RBB(t)=0.
Furthermore, substituting (44) in the first equation of (27) shows that
(48)SB(t)⟶ΠBμBast⟶∞.
Thus, by combining equations (44), (47), and (48), it follows that every solution of the equations of the buffalo-only model (27), with θEB=θRB=0 and initial conditions in Γ, approaches E0 as t→∞ (whenever R0≤1).Theorem6 shows that, in the absence of the reinfection of exposed and recovered buffalos (i.e., θEB=θRB=0), BTB can be eliminated from the buffalo-only population if the reproduction number of the model (R0) can be brought to (and maintained at) a value less than unity. Figure 4(a) depicts the solution profiles of the buffalo-only model (27), generated using various initial conditions, showing convergence to the DFE E0 when R0<1 (in line with Theorem 6).Simulations of the buffalo-only model (27), showing the total number of infected buffalos with clinical symptoms of BTB (IBB(t)(t)) at time t as a function of time. Parameter values used are as given in Table 3 with (a) βB=0.00733 (so that, R0=0.7036<1) and (b) βB= 0.0733, δB=0 (so that, R~0=8.6050>1).
(a)
(b)


## 3.3. Existence of Endemic Equilibria: Special Case

In this section, the existence of nontrivial (endemic) equilibria (where the components of the infected variables of the model are nonzero) of the buffalo-only model (27) is explored for the special case without reinfection (i.e., θEB=θRB=0). Solving the equations of the buffalo-only model (27) at steady-state gives the following general form of the endemic equilibrium (denoted by E1):
(49)E1=(SB**,EB1**,EB2**,IBB**,RBB**),
where
(50)SB**=ΠBλB**+μB,EB1**=λB**ΠBK1(λB**+μB),EB2**=κ1λB**ΠBK1K3(λB**+μB),IBB**=σB2κ1λB**ΠBK1K3K5(λB**+μB),RBB**=γB1σB2κ1λB**ΠBK1K3K5μB(λB**+μB),
with the force of infection at steady-state (λB**) given by
(51)λB**=βBNB**(ηB1EB1**+ηB2EB2**+IBB**).Using (50) in the expression for λB** in (51) shows that the nonzero equilibrium of the model (22) satisfies the linear equation:
(52)b1λB**+b2=0,
where b1=K5μB(K3+κ1)+σB2κ1(μB+γB2) and b2=K1K3K5μB(1-R0). Clearly, the coefficient b1 is always positive, and b2 is positive (negative) if R0 is less than (greater than) unity, respectively. Thus, the linear system (52) has a unique positive solution, given by λB**=-b2/b1, whenever R0>1. Further, the force of infection for buffalos (λB**) is negative whenever R0<1 (which is biologically meaningless). Hence, the buffalo-only model (27) has no positive equilibrium in this case. These results are summarized below.Theorem 7.
The buffalo-only model (27), with θEB=θRB=0, has a unique endemic equilibrium, E1, whenever R0>1, and no endemic equilibrium otherwise.


### 3.3.1. Global Asymptotic Stability of Endemic Equilibrium

The global asymptotic stability of the unique endemic equilibrium (E1) of the buffalo-only model is explored for the special case without reinfection (θEB=θRB=0) and BTB-induced death in buffalos (δB=0). It is convenient to define
(53)Γ1={(SB,EB1,EB2,IBB,RBB)∈Γ:EB1=EB2iiiii=IBB=RBB=0},
the stable manifold of the DFE (E0) of the buffalo-only model (27).Theorem 8.
The unique endemic equilibrium (E1) of the buffalo-only model (27), with θEB=θRB=δB=0, is GAS in Γ∖Γ1 if R~0=R0|δB=0>1.Proof.
Consider the buffalo-only model (27) with θEB=θRB=δB=0. For this case, it follows from Theorem 7 that the buffalo-only model (27) has a unique endemic equilibrium whenever R~0>1. Furthermore, setting δB=0 in model (27) shows that NB(t)→ΠB/μB as t→∞. Consider the following nonlinear Lyapunov function (of Goh-Volterra type) for the subsystem of model (27) involving the state variables SB, EB1, EB2, and IBB (noting that NB(t) is now replaced by its limiting value ΠB/μB):
(54)F=SB-SB**-SB**ln⁡(SBSB**)+EB1-EB1**-EB1**ln⁡(EB1EB1**)+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[EB2-EB2**-EB2**ln⁡(EB2EB2**)]+β~BSB**IBB**σB2EB2**[IBB-IBB**-IBB**ln⁡(IBBIBB**)],
where β~B=μBβB/ΠB. The Lyapunov derivative of F is given by
(55)F˙=S˙B-SB**SBS˙B+E˙B1-EB1**EB1E˙B1+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)(E˙B2-EB2**EB2E˙B2)+β~BSB**IBB**σB2EB2**(I˙BB-IBB**IBBI˙BB),=ΠB-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB-SB**SB[ΠB-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB]+β~B(ηB1EB1+ηB2EB2+IBB)SB-K1EB1-EB1**EB1[β~B(ηB1EB1+ηB2EB2+IBB)SB-K1EB1]+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[κ1EB1-K3EB2-EB2**EB2(κ1EB1-K3EB2)]+β~BSB**IBB**σB2EB2**[σB2EB2-K5IBB-IBB**IBBDDDDDDDD×(σB2EB2-K5IBB)σB2EB2-K5IBB-IBB**IBB].
Using the following steady-state relations (obtained from (27)),
(56)ΠB=β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**,κ1EB1**=K3EB2**,β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**=K1EB1**,σB2EB2**=K5IBB**,γB1IBB**=μBRBB**,
the Lyapunov derivative can be simplified to
(57)F˙=β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**-μBSB-SB**SB[β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**iiiiiiiiiiii-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB]-K1EB1-EB1**EB1[β~B(ηB1EB1+ηB2EB2+IBB)SBddddddddddddiiiii-K1EB1[β~B(ηB1EB1+ηB2EB2+IBB)SB]+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[κ1EB1-K3EB2-EB2**EB2(κ1EB1-K3EB2)]+β~BSB**IBB**σB2EB2**×[σB2EB2-K5IBB-IBB**IBB(σB2EB2-K5IBB)].
Thus,(58)F˙=μBSB**(2-SB**SB-SBSB**)+β~BηB1EB1**SB**(2-SB**SB-SBSB**)+β~BηB2EB2**SB**(3-SB**SB-EB1EB2**EB1**EB2-SBEB1**EB2SB**EB1EB2**)+β~BIBB**SB**×(4-SB**SB-EB1EB2**EB1**EB2-EB2IBB**EB2**IBB-SBEB1**IBBSB**EB1IBB**).
Finally, since the arithmetic mean exceeds the geometric mean, it follows then that
(59)μBSB**(2-SB**SB-SBSB**)≤0,β~BηB1EB1**SB**(2-SB**SB-SBSB**)≤0,β~BηB2EB2**SB**(3-SB**SB-EB1EB2**EB1**EB2-SBEB1**EB2SB**EB1EB2**)≤0,β~BIBB**SB**(4-SB**SB-EB1EB2**EB1**EB2-EB2IBB**EB2**IBB-SBEB1**IBBSB**EB1IBB**)≤0.
Furthermore, since all the model parameters are nonnegative, it follows that F˙≤0 for R~0>1. Thus, F˙ is a Lyapunov function for the subsystem of model (27) on Γ∖Γ1. Therefore, it follows, by LaSalle’s Invariance Principle [37], that
(60)lim⁡t→∞SB(t)=SB**,lim⁡t→∞EB1(t)=EB1**,lim⁡t→∞EB2(t)=EB2**,lim⁡t→∞IBB(t)=IBB**.
Since IBB(t)→IBB** as t→∞, it follows from the equation for dRBB/dt in (27) that RBB(t)→γB1IBB**/μB=RBB** as t→∞. The proof is concluded using similar arguments as in the proof of Theorem 6.The epidemiological implication of Theorem8 is that BTB will be endemic in the buffalo population if R~0>1 (and θEB=θRB=δB=0). Figure 4(b) depicts the solutions of model (27) for the case when R~0>1 and θEB=θRB=δB=0, showing convergence of the initial solutions to the unique endemic equilibrium (in line with Theorem 8).


## 3.3.1. Global Asymptotic Stability of Endemic Equilibrium

The global asymptotic stability of the unique endemic equilibrium (E1) of the buffalo-only model is explored for the special case without reinfection (θEB=θRB=0) and BTB-induced death in buffalos (δB=0). It is convenient to define
(53)Γ1={(SB,EB1,EB2,IBB,RBB)∈Γ:EB1=EB2iiiii=IBB=RBB=0},
the stable manifold of the DFE (E0) of the buffalo-only model (27).Theorem 8.
The unique endemic equilibrium (E1) of the buffalo-only model (27), with θEB=θRB=δB=0, is GAS in Γ∖Γ1 if R~0=R0|δB=0>1.Proof.
Consider the buffalo-only model (27) with θEB=θRB=δB=0. For this case, it follows from Theorem 7 that the buffalo-only model (27) has a unique endemic equilibrium whenever R~0>1. Furthermore, setting δB=0 in model (27) shows that NB(t)→ΠB/μB as t→∞. Consider the following nonlinear Lyapunov function (of Goh-Volterra type) for the subsystem of model (27) involving the state variables SB, EB1, EB2, and IBB (noting that NB(t) is now replaced by its limiting value ΠB/μB):
(54)F=SB-SB**-SB**ln⁡(SBSB**)+EB1-EB1**-EB1**ln⁡(EB1EB1**)+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[EB2-EB2**-EB2**ln⁡(EB2EB2**)]+β~BSB**IBB**σB2EB2**[IBB-IBB**-IBB**ln⁡(IBBIBB**)],
where β~B=μBβB/ΠB. The Lyapunov derivative of F is given by
(55)F˙=S˙B-SB**SBS˙B+E˙B1-EB1**EB1E˙B1+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)(E˙B2-EB2**EB2E˙B2)+β~BSB**IBB**σB2EB2**(I˙BB-IBB**IBBI˙BB),=ΠB-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB-SB**SB[ΠB-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB]+β~B(ηB1EB1+ηB2EB2+IBB)SB-K1EB1-EB1**EB1[β~B(ηB1EB1+ηB2EB2+IBB)SB-K1EB1]+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[κ1EB1-K3EB2-EB2**EB2(κ1EB1-K3EB2)]+β~BSB**IBB**σB2EB2**[σB2EB2-K5IBB-IBB**IBBDDDDDDDD×(σB2EB2-K5IBB)σB2EB2-K5IBB-IBB**IBB].
Using the following steady-state relations (obtained from (27)),
(56)ΠB=β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**,κ1EB1**=K3EB2**,β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**=K1EB1**,σB2EB2**=K5IBB**,γB1IBB**=μBRBB**,
the Lyapunov derivative can be simplified to
(57)F˙=β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**-μBSB-SB**SB[β~B(ηB1EB1**+ηB2EB2**+IBB**)SB**+μBSB**iiiiiiiiiiii-β~B(ηB1EB1+ηB2EB2+IBB)SB-μBSB]-K1EB1-EB1**EB1[β~B(ηB1EB1+ηB2EB2+IBB)SBddddddddddddiiiii-K1EB1[β~B(ηB1EB1+ηB2EB2+IBB)SB]+(β~BηB2SB**EB2**+β~BSB**IBB**κ1EB1**)×[κ1EB1-K3EB2-EB2**EB2(κ1EB1-K3EB2)]+β~BSB**IBB**σB2EB2**×[σB2EB2-K5IBB-IBB**IBB(σB2EB2-K5IBB)].
Thus,(58)F˙=μBSB**(2-SB**SB-SBSB**)+β~BηB1EB1**SB**(2-SB**SB-SBSB**)+β~BηB2EB2**SB**(3-SB**SB-EB1EB2**EB1**EB2-SBEB1**EB2SB**EB1EB2**)+β~BIBB**SB**×(4-SB**SB-EB1EB2**EB1**EB2-EB2IBB**EB2**IBB-SBEB1**IBBSB**EB1IBB**).
Finally, since the arithmetic mean exceeds the geometric mean, it follows then that
(59)μBSB**(2-SB**SB-SBSB**)≤0,β~BηB1EB1**SB**(2-SB**SB-SBSB**)≤0,β~BηB2EB2**SB**(3-SB**SB-EB1EB2**EB1**EB2-SBEB1**EB2SB**EB1EB2**)≤0,β~BIBB**SB**(4-SB**SB-EB1EB2**EB1**EB2-EB2IBB**EB2**IBB-SBEB1**IBBSB**EB1IBB**)≤0.
Furthermore, since all the model parameters are nonnegative, it follows that F˙≤0 for R~0>1. Thus, F˙ is a Lyapunov function for the subsystem of model (27) on Γ∖Γ1. Therefore, it follows, by LaSalle’s Invariance Principle [37], that
(60)lim⁡t→∞SB(t)=SB**,lim⁡t→∞EB1(t)=EB1**,lim⁡t→∞EB2(t)=EB2**,lim⁡t→∞IBB(t)=IBB**.
Since IBB(t)→IBB** as t→∞, it follows from the equation for dRBB/dt in (27) that RBB(t)→γB1IBB**/μB=RBB** as t→∞. The proof is concluded using similar arguments as in the proof of Theorem 6.The epidemiological implication of Theorem8 is that BTB will be endemic in the buffalo population if R~0>1 (and θEB=θRB=δB=0). Figure 4(b) depicts the solutions of model (27) for the case when R~0>1 and θEB=θRB=δB=0, showing convergence of the initial solutions to the unique endemic equilibrium (in line with Theorem 8).


## 3.4. Sensitivity and Uncertainty Analyses

In this section, sensitivity and uncertainty analyses will be carried out, using Latin hypercube sampling (LHS) and partial correlation coefficient (PRCC) [39–41], to assess the effect of uncertainty in the estimate of the parameter values used to simulate the buffalo-only model (on the simulation results obtained) and to determine the key parameters that drive the dynamics of the disease in the buffalo-human population. The ranges and baseline values of the parameters of the buffalo-only model, given in Table 3 with θRB=θEB=0 (i.e., in the absence of backward bifurcation), will be used in these analyses. Each parameter of the buffalo-only model (22) is assumed to obey a uniform distribution [42]. Following [42], a total of 1000 LHS runs (N=1000) are carried out. Furthermore, the following initial conditions (which are consistent with the dynamics of African buffalo in the Kruger National Park [29]): (SB(0), EB1(0), EM1(0), EB2(0), EM2(0), IBB(0), IMB(0), RBB(0), RMB(0))=(28000,100,100,20,20,10,10,100,100)) are used in the simulations.Table 3
Ranges and baseline values for parameters of the BTB-MTB model (21).
Parameter
Range (day−1)
Baseline value (day−1)
Reference
Π
H
[26, 80]
53
[20]
Π
B
[2, 4]
3
[20–22]
μ
H
[0.0000274, 0.0000549]
0.000047
[14, 23–25]
μ
B
(0.00009477, 0.00011583)
0.0001053
[22, 26]
β
H
[0.00011, 0.000959]
0.000535
[14, 17]
β
B
(0.006597, 0.008063)
0.00733
[10]
η
H
1
[0, 1)
0.5
Fitted
η
H
2
[0, 1)
0.5
Fitted
η
B
1
(0.4455, 0.5045)
0.45
Fitted
η
B
2
(0.495, 0.605)
0.55
Fitted
θ
B
B
[0, 1)
0.5
Fitted
θ
M
M
[0, 1]
0.5
Assumed
θ
H
H
[0, 1)
0.5
Assumed
γ
i
(
i
=
1
,
2
)
(0.0000823, 0.000823)
0.000453
[14, 23]
γ
B
1
(0.00774, 0.00946)
0.0086
[10]
γ
M
1
(0.13374, 0.160086)
0.1486
[12]
σ
i
(
i
=
1
,
2
)
(0.0000822, 0.00247)
0.001
[23, 27]
κ
1
(0.45, 0.55)
0.5
[27]
κ
2
(0.45, 0.55)
0.5
[27]
σ
B
2
(0.25, 0.35)
0.3
[27]
σ
M
2
(0.36, 0.44)
0.4
[27]
θ
H
i
(
i
=
1
,
2
)
[0, 0.1]
0.00271
[23]
θ
R
B
,
θ
R
H
[0.002439, 0.002981]
0.00271
[23]
θ
E
B
[0.002439, 0.002981]
0.00271
[23]
δ
H
1
,
δ
H
2
[0.000115, 0.000822]
0.0002
[14, 17, 23, 25]
δ
B
[0.0018, 0.0022]
0.002
[28]
δ
M
[0.0018, 0.0022]
0.002
[14, 23]Table 4
Number of symptomatic buffalos with BTB at Kruger National Park [29].
Year
Number of symptomatic buffalos [29]
2001
35
2002
135
2003
185
2004
238
2005
230Figure5 depicts a box plot of R0, as a function of the number of LHS runs carried out (N=1000), from which it is evident that the distribution of R0 lies in the range R0∈[0.25,0.75] (each box plot displays the upper and lower quartile ranges of R0. a horizontal line within the box is the median value, and values of R0 beyond the whiskers are outliers [43]). Thus, since the distribution of the reproduction number of the buffalo-only model is less than unity, it follows (from Theorem 6) that the BTB-MTB outbreaks (in the buffalo-human population) will die out with time (in other words, the disease will be effectively controlled). The PRCC values of the parameters of the buffalo-only model (22), using R0 as the response function, are depicted in Figure 6. It follows from Figure 6 that the top three parameters that most influences the value of R0 (hence the disease dynamics) are the BTB transmission rate (βB), the recovery rate of buffalos (γB1), and the BTB-induced mortality in buffalos (δB).Figure 5
Box plot ofR0 as a function of the number of LHS runs carried out for the buffalo-only model (22), using parameter values and ranges given in Table 3 with θEB=θRB=0.Figure 6
PRCC values of the parameters of the buffalo-only model (22), using R0 as the output function. Parameter values used are as given in Table 3.Similarly, Figure7 depicts the box plot of the buffalo-only model (22) using total number of symptomatic buffalos (IBB+IMB) as the response function. This figure shows a distribution of the number of symptomatic buffalos lying in the range [20–130]. Hence, this study shows that, using the parameter values and ranges relevant to BTB-MTB dynamics at the Kruger National Park, a BTB outbreak could cause no more than 130 confirmed cases (of BTB and MTB) in the park. The associated PRCC values (with the total number of symptomatic buffalos as the output) are depicted in Figure 8, from which it is evident that, in this scenario, the top three parameters (that most influences the output) are the buffalo recruitment rate (ΠB), the natural (μB), and the disease-induced (δB) death rate of buffalos. Hence, this study shows variability in the top-ranked PRCC values on the chosen response/output function.Figure 7
Box plot of the total number of symptomatic buffalos (IBB+IMB) as a function of the number of LHS runs for the buffalo-only model (22), using parameter values and ranges given in Table 3 with θEB=θRB=0.Figure 8
PRCC values of the parameters of the buffalo-only model (22), using total number of symptomatic buffalos (IBB+IMB) as the output function. Parameter values used are as given in Table 3.Having fully studied the dynamics of the bovine-only model (22), the full BTB-MTB model (22) will now be analysed.


## 4. Analysis of the BTB-MTB Model

It can be shown, using the approach in Section3, that the following biologically feasible region
(61)Ω={ΠBμB(SH,EH1,EH2,IH1,IH2,RH1,RH2,SB,iiiiiiEB1,EM1,EB2,EM2,IBB,IMB,RBB,RMB)iiiiiiii∈R+16:NH≤ΠHμH,NB≤ΠBμB}
is positively invariant and attracting for the BTB-MTB model (21).


### 4.1. Local Stability of DFE

The analyses in this section will be carried out for the special case of the BTB-MTB model (21) with θMM=θBB=0. The justification for this assumption is based on the fact that contact between humans and buffalos in the Kruger National Park are tightly controlled (hence, it is reasonable to assume that buffalo-to-human or human-to-buffalo transmission of BTB is negligible). The DFE of the BTB-MTB model (21) is given by
(62)E0f=(SH*,EH1*,EH2*,IH1*,IH2*,RH1*,RH2*,SB*,iiiiiEB1*,EM1*,EB2*,EM2*,IBB*,IMB*,RBB*,RMB*)=(ΠHμH,0,0,0,0,0,0,ΠBμB,0,0,0,0,0,0,0,0).
The associated next generation matrices of the BTB-MTB model (21), denoted by Ff and Vf, are given, respectively, by(63)Ff=[βHηH10βH00000000βHηH20βH000000000000000000000000000000βBηB10βBηB20βB000000000000000000000000000000000000000000000000000],Vf=[Q10000000000Q200000000-σ10Q300000000-σ20Q40000000000K10000000000K200000000-κ10K3000000000-κ2K400000000-σB20K500000000-σM20K6].It follows then that thereproduction number of the BTB-MTB model (21), denoted by Rf, is given by
(64)Rf=ρ(FfVf-1)=max⁡{RHM,RHB,R0},
where RHM and RHB are the associatedreproduction numbers for humans infected with MTB and with BTB, respectively, given by
(65)RHM=βH(ηH1Q3+σ1)Q1Q3,RHB=βH(ηH2Q4+σ2)Q2Q4,
where Q1=σ1+μH, Q2=σ2+μH, Q3=γ1+μH+δH1, and Q4=γ2+μH+δH2, and R0 is as defined in Section 3. Thus, using the approach in Section 3.1, the following result can be established for the BTB-MTB model (21).Lemma 9.
The DFE,E0f, of model (21), with θMM=θBB=0, is LAS in Ω if Rf<1, and unstable if Rf>1.It can be shown, as in Section3, that the BTB-MTB model (21) also undergoes backward bifurcation. Unlike in the buffalo-only model (22), however, this phenomenon persists even if the bovine-associated reinfection terms (θRB and θEB) are set to zero. This is due to the reinfection of exposed and recovered humans (i.e., θH1≠0 and θH2≠0). To illustrate this fact, it is shown that the DFE (E0f) of the BTB-MTB model (21) is GAS in Ω in the absence of reinfection of exposed and recovered buffalos and humans, whenever the associated reproduction number (Rf) is less than unity.


### 4.2. Global Asymptotic Stability of DFE

Theorem 10.
The DFE,E0f, of the BTB-MTB model (21) with θH1=θH2=θRB=θRH=θBB=θMM=θEB=θRB=0, is GAS in Ω if Rf<1.Proof.
The proof of Theorem10, based on using comparison theorem [44], is given in Appendix B. Hence, the analyses in this section show that the buffalo-only model and the full BTB-MTB model (21) have essentially the same qualitative dynamics with respect to the local and global asymptotic stability of the associated disease-free equilibrium (in the absence of reinfection) as well as the backward bifurcation property established in the transmission dynamics of BTB and BTB-MTB in a buffalo-human population. In both cases, the backward bifurcation phenomenon is shown to arise due to the reinfection of the exposed and recovered host(s) (buffalos for the buffalo-only model (22), and buffalos and humans for the BTB-MTB model).


### 4.3. Numerical Simulations

The BTB-MTB model (21) is simulated, using the baseline values tabulated in Table 3 (unless otherwise stated), to assess the effect of the dynamics of BTB (MTB) on the spread of MTB (BTB) in the human (buffalo) population.


#### 4.3.1. Effect of BTB on MTB

The effect of BTB (in the human-buffalo population within the Kruger National Park) on the spread of MTB in the human population (within the park) is assessed by simulating the BTB-MTB model (21) using parameter values in Table 3, subject to the following four effectiveness levels of BTB transmission likelihood from buffalos to humans (i.e., choosing four different values of the parameter θMM, for the reduced likelihood of humans acquiring BTB infection from buffalos):(i)
no transmission of BTB from buffalos to humans:θMM=0,(ii)
low rate of transmission of BTB from buffalos to humans:θMM=0.25,(iii)
moderate rate of transmission of BTB from buffalos to humans:θMM=0.50,(iv)
high rate of transmission of BTB from buffalos to humans:θMM=0.75.The simulation results obtained, depicted in Figure9(a), show that the cumulative number of new MTB cases in humans decreases with increasing rate of BTB transmission to humans by buffalos (θMM).Cumulative number of new cases of (a) MTB infection in humans and (b) BTB infection in buffalos. Parameter values used are as given in Table3, with various values of θMM (a) or θHH (b).
(a)
(b)


#### 4.3.2. Effect of MTB on BTB

Similar plot is generated to assess the effect of MTB (in the human-buffalo population) on the spread of BTB in the buffalo population. Here, too, four transmission levels of the associated parameter(θHH) are considered, namely, none (θHH=0), low (θHH=0.25), moderate (θHH=0.50), and high (θHH=0.75). The results obtained, depicted in Figure 9(b), show that the cumulative number of new BTB infections in buffalos decreases with increasing rate of MTB transmission to buffalos by humans.


## 4.1. Local Stability of DFE

The analyses in this section will be carried out for the special case of the BTB-MTB model (21) with θMM=θBB=0. The justification for this assumption is based on the fact that contact between humans and buffalos in the Kruger National Park are tightly controlled (hence, it is reasonable to assume that buffalo-to-human or human-to-buffalo transmission of BTB is negligible). The DFE of the BTB-MTB model (21) is given by
(62)E0f=(SH*,EH1*,EH2*,IH1*,IH2*,RH1*,RH2*,SB*,iiiiiEB1*,EM1*,EB2*,EM2*,IBB*,IMB*,RBB*,RMB*)=(ΠHμH,0,0,0,0,0,0,ΠBμB,0,0,0,0,0,0,0,0).
The associated next generation matrices of the BTB-MTB model (21), denoted by Ff and Vf, are given, respectively, by(63)Ff=[βHηH10βH00000000βHηH20βH000000000000000000000000000000βBηB10βBηB20βB000000000000000000000000000000000000000000000000000],Vf=[Q10000000000Q200000000-σ10Q300000000-σ20Q40000000000K10000000000K200000000-κ10K3000000000-κ2K400000000-σB20K500000000-σM20K6].It follows then that thereproduction number of the BTB-MTB model (21), denoted by Rf, is given by
(64)Rf=ρ(FfVf-1)=max⁡{RHM,RHB,R0},
where RHM and RHB are the associatedreproduction numbers for humans infected with MTB and with BTB, respectively, given by
(65)RHM=βH(ηH1Q3+σ1)Q1Q3,RHB=βH(ηH2Q4+σ2)Q2Q4,
where Q1=σ1+μH, Q2=σ2+μH, Q3=γ1+μH+δH1, and Q4=γ2+μH+δH2, and R0 is as defined in Section 3. Thus, using the approach in Section 3.1, the following result can be established for the BTB-MTB model (21).Lemma 9.
The DFE,E0f, of model (21), with θMM=θBB=0, is LAS in Ω if Rf<1, and unstable if Rf>1.It can be shown, as in Section3, that the BTB-MTB model (21) also undergoes backward bifurcation. Unlike in the buffalo-only model (22), however, this phenomenon persists even if the bovine-associated reinfection terms (θRB and θEB) are set to zero. This is due to the reinfection of exposed and recovered humans (i.e., θH1≠0 and θH2≠0). To illustrate this fact, it is shown that the DFE (E0f) of the BTB-MTB model (21) is GAS in Ω in the absence of reinfection of exposed and recovered buffalos and humans, whenever the associated reproduction number (Rf) is less than unity.


## 4.2. Global Asymptotic Stability of DFE

Theorem 10.
The DFE,E0f, of the BTB-MTB model (21) with θH1=θH2=θRB=θRH=θBB=θMM=θEB=θRB=0, is GAS in Ω if Rf<1.Proof.
The proof of Theorem10, based on using comparison theorem [44], is given in Appendix B. Hence, the analyses in this section show that the buffalo-only model and the full BTB-MTB model (21) have essentially the same qualitative dynamics with respect to the local and global asymptotic stability of the associated disease-free equilibrium (in the absence of reinfection) as well as the backward bifurcation property established in the transmission dynamics of BTB and BTB-MTB in a buffalo-human population. In both cases, the backward bifurcation phenomenon is shown to arise due to the reinfection of the exposed and recovered host(s) (buffalos for the buffalo-only model (22), and buffalos and humans for the BTB-MTB model).


## 4.3. Numerical Simulations

The BTB-MTB model (21) is simulated, using the baseline values tabulated in Table 3 (unless otherwise stated), to assess the effect of the dynamics of BTB (MTB) on the spread of MTB (BTB) in the human (buffalo) population.


### 4.3.1. Effect of BTB on MTB

The effect of BTB (in the human-buffalo population within the Kruger National Park) on the spread of MTB in the human population (within the park) is assessed by simulating the BTB-MTB model (21) using parameter values in Table 3, subject to the following four effectiveness levels of BTB transmission likelihood from buffalos to humans (i.e., choosing four different values of the parameter θMM, for the reduced likelihood of humans acquiring BTB infection from buffalos):(i)
no transmission of BTB from buffalos to humans:θMM=0,(ii)
low rate of transmission of BTB from buffalos to humans:θMM=0.25,(iii)
moderate rate of transmission of BTB from buffalos to humans:θMM=0.50,(iv)
high rate of transmission of BTB from buffalos to humans:θMM=0.75.The simulation results obtained, depicted in Figure9(a), show that the cumulative number of new MTB cases in humans decreases with increasing rate of BTB transmission to humans by buffalos (θMM).Cumulative number of new cases of (a) MTB infection in humans and (b) BTB infection in buffalos. Parameter values used are as given in Table3, with various values of θMM (a) or θHH (b).
(a)
(b)


### 4.3.2. Effect of MTB on BTB

Similar plot is generated to assess the effect of MTB (in the human-buffalo population) on the spread of BTB in the buffalo population. Here, too, four transmission levels of the associated parameter(θHH) are considered, namely, none (θHH=0), low (θHH=0.25), moderate (θHH=0.50), and high (θHH=0.75). The results obtained, depicted in Figure 9(b), show that the cumulative number of new BTB infections in buffalos decreases with increasing rate of MTB transmission to buffalos by humans.


## 4.3.1. Effect of BTB on MTB

The effect of BTB (in the human-buffalo population within the Kruger National Park) on the spread of MTB in the human population (within the park) is assessed by simulating the BTB-MTB model (21) using parameter values in Table 3, subject to the following four effectiveness levels of BTB transmission likelihood from buffalos to humans (i.e., choosing four different values of the parameter θMM, for the reduced likelihood of humans acquiring BTB infection from buffalos):(i)
no transmission of BTB from buffalos to humans:θMM=0,(ii)
low rate of transmission of BTB from buffalos to humans:θMM=0.25,(iii)
moderate rate of transmission of BTB from buffalos to humans:θMM=0.50,(iv)
high rate of transmission of BTB from buffalos to humans:θMM=0.75.The simulation results obtained, depicted in Figure9(a), show that the cumulative number of new MTB cases in humans decreases with increasing rate of BTB transmission to humans by buffalos (θMM).Cumulative number of new cases of (a) MTB infection in humans and (b) BTB infection in buffalos. Parameter values used are as given in Table3, with various values of θMM (a) or θHH (b).
(a)
(b)


## 4.3.2. Effect of MTB on BTB

Similar plot is generated to assess the effect of MTB (in the human-buffalo population) on the spread of BTB in the buffalo population. Here, too, four transmission levels of the associated parameter(θHH) are considered, namely, none (θHH=0), low (θHH=0.25), moderate (θHH=0.50), and high (θHH=0.75). The results obtained, depicted in Figure 9(b), show that the cumulative number of new BTB infections in buffalos decreases with increasing rate of MTB transmission to buffalos by humans.


## 5. Conclusions

A new deterministic model for the transmission dynamics of BTB and MTB in a community of humans and buffalos is designed and rigorously analyzed. Some of the main findings of the study are as follows.(i)
The buffalo-only model undergoes the phenomenon of backward bifurcation. This phenomenon is caused by the exogenous reinfection of exposed and infected buffalos. In the absence of reinfection, the disease-free equilibrium of the buffalo-only model is shown to be globally asymptotically stable whenever the associated reproduction number of the model is less than unity.(ii)
In the absence of the reinfection of exposed and recovered buffalos (θEB=θRB=0), the buffalo-only model is shown to have unique endemic equilibrium whenever its reproduction number exceeds unity (R0>1). This equilibrium is shown to be globally asymptotically stable for the special case where the disease-induced mortality in buffalos is negligible (δB=0).(iii)
Detailed uncertainty analyses of the buffalo-only model, using a reasonable set of parameter values and ranges (Table3) relevant to BTB dynamics in the Kruger National Park, shows that the distribution of the associated reproduction number of the buffalo-only model is less than unity (hence, BTB outbreaks will not persist in the Park). Furthermore, such outbreak would cause no more than 120 confirmed (symptomatic) cases of BTB or MTB within the Park. Sensitivity analysis, for the case when the reproduction number (R0) is chosen as the response/output function, reveals that the three main parameters that govern the disease dynamics are the BTB transmission rate, recovery rate of buffalos, and BTB-induced mortality rate. Similarly, three parameters (recruitment rate of buffalos, natural, and BTB-induced death rates in buffalos) are identified as the main influential parameters for the case where the number of symptomatic buffalos (with BTB or MTB) is the chosen output function.(iv)
The BTB-MTB model also undergoes backward bifurcation. Unlike in the buffalo-only model, this phenomenon persists even if the bovine-associated reinfection terms (θRB and θEB) are set to zero. This is due to the reinfection of exposed and recovered humans (θH1≠0, θH2≠0, θRH≠0, θRB≠0). It is shown that this model does not undergo backward bifurcation in the absence of reinfection of exposed and recovered host(s) (buffalos and humans). For this case, it is shown that the DFE of the BTB-MTB model (21) is globally asymptotically stable, whenever the associated reproduction number is less than unity.(v)
The buffalo-only model and the full BTB-MTB model exhibit the same qualitative dynamics with respect to the local and global asymptotic stability of the associated disease-free equilibrium (in the absence of reinfection of the associated host(s)). In both models, the backward bifurcation phenomenon is shown to arises due to the reinfection of exposed and recovered host(s).(vi)
Numerical simulations of the BTB-MTB model show that an increase in the cumulative number of BTB infection leads to a marked reduction in the cumulative number of new MTB cases in humans. Similarly, an increase in the cumulative number of MTB infection led to a significant decrease in the cumulative number of new BTB cases in buffalos.



---

*Source: 912306-2014-09-02.xml*