Dynamics of Mycobacterium and bovine tuberculosis in a human-buffalo population.

A new model for the transmission dynamics of Mycobacterium 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 of Mycobacterium tuberculosis cases in humans (buffalos) decreases with increasing number of bovine tuberculosis infections in humans (buffalo).

1. Introduction

Mycobacterium tuberculosis (MTB) and bovine tuberculosis (BTB) are chronic bacterial diseases, classified amongst the closely related species that form the M. tuberculosis complex (MTBC) [1]. The human MTB is caused by tubercle bacillus (M. tuberculosis), while BTB is caused by bovine 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 million per 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 controlled via 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].

Figure 1

African buffalos and demographic map of Kruger National Park [11].

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 Section 2. 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 time t, denoted by N (t), is subdivided into seven mutually exclusive compartments of susceptible humans (S (t)), exposed humans (who have been infected with MTB but have not yet shown clinical symptoms of the disease) (E (t)), exposed humans with BTB (E (t)), infected humans with clinical symptoms of MTB (I (t)), infected humans with clinical symptoms of BTB (I (t)), and humans who recovered from MTB (R (t)) or BTB (R (t)), so that Similarly, the total buffalo population (in the herd) at time t, denoted by N (t), is split into susceptible (S (t)), early-exposed with BTB (E (t)), early-exposed with MTB (E (t)), advanced-exposed with BTB (E (t)), advanced-exposed with MTB (E (t)), infected with clinical symptoms of BTB (I (t)), infected with clinical symptoms of MTB (I (t)), and recovered from BTB (R (t)) or MTB (R (t)), so that The susceptible human population (S (t)) is increased by the recruitment of people (either by birth or immigration) into the human-buffalo poulation (at a rate Π). The population is decreased by infection with MTB (at a rate λ ) or BTB (at a rate λ ), where with In (3) and (4), β and β represent the effective contact rates (i.e., contacts capable of leading to MTB or BTB infection), respectively. Furthermore, 0 ≤ η < 1 and 0 ≤ η < 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 ≤ η < 1 and 0 ≤ η < 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 ≤ θ < 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 μ . Thus, the rate of change of the susceptible human population is given by The population of exposed humans with MTB (E (t)) is generated by the infection of susceptible humans with MTB (at the rate λ ) and is decreased by the development of clinical symptoms of MTB (at a rate σ 1), exogenous reinfection (at a rate θ λ , where 0 ≤ θ < 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 Similarly, the population of exposed humans with BTB (E (t)) is increased by the infection of susceptible humans with BTB (at the rate λ ) and is reduced by the development of clinical symptoms of BTB (at a rate σ 2), exogenous reinfection (at a rate θ λ , with 0 ≤ θ < 1 similarly defined as θ ) and natural death. Thus, The population of humans with clinical symptoms of MTB (I (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 θ λ and θ λ , resp., with 0 ≤ θ < 1). This population is decreased by recovery (at a rate γ 1), natural death, and MTB-induced death (at a rate δ ), so that The population of infected humans with clinical symptoms of BTB (I (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 θ λ and θ λ , resp., with 0 ≤ θ < 1). This population is decreased by recovery (at a rate γ 2), natural death, and BTB-induced death (at a rate δ ). This gives The population of humans who recovered from MTB (R (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 θ λ ) and natural death. Hence, 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 (R (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 θ λ ) and natural death, so that The population of susceptible buffalos (S (t)) is generated by the recruitment of buffalos (either by birth or restocking from other herds) at a rate Π. 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 λ (where, λ = θ λ + λ , with the modification parameters 0 ≤ θ < 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 θ λ (where 0 ≤ θ < 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 μ , buffalos in each epidemiological compartment suffer natural death at this rate). Thus, 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 (E (t)) is increased by the infection of susceptible buffalos with BTB (at the rate λ ). This population is decreased by exogenous reinfection with BTB (at a rate θ λ , with 0 ≤ θ < 1), progression to the advanced-exposed class (at a rate κ 1), and natural death. This gives

The population of buffalos early-exposed to MTB is increased by the infection of susceptible buffalos with MTB (at the rate θ λ , where 0 ≤ θ < 1 is as defined above). The population is decreased by exogenous reinfection (at a rate θ λ ), progression to the advanced-exposed MTB class (at a rate κ 2), and natural death. This gives

The population of buffalos at advanced-exposed BTB class (E (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 θ λ ), development of clinical symptoms of BTB (at a rate σ ), and natural death, so that Similarly, the population of buffalos at advanced-exposed MTB class (E (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 θ λ ), development of clinical symptoms of MTB (at a rate σ ), and natural death. Hence, The population of buffalos with clinical symptoms of BTB (I (t)) is increased by the development of clinical symptoms of exposed buffalos with BTB (at the rate σ ) and by the exogenous reinfection of exposed and recovered buffalos (at the rates θ λ and θ λ , resp.). It is decreased by recovery (at a rate γ ), natural death, and BTB-induced mortality (at a rate δ ). Thus, The population of buffalos with clinical symptoms of MTB (I (t)) is increased by the development of clinical symptoms of exposed buffalos with MTB (at the rate σ ) and by the exogenous reinfection of exposed and recovered buffalos (at the rates θ λ and θ λ , resp.). It is decreased by recovery (at a rate γ ), natural death, and MTB-induced mortality (at a rate δ ). Thus, The population of buffalos who recovered from BTB (R (t)) is increased following the recovery of buffalos with clinical symptoms of BTB (at the rate γ ). It is decreased by reinfection (at the rate θ λ ) and natural death, so that Finally, the population of buffalos who recovered from MTB (R (t)) is generated by the recovery of buffalos with MTB (at the rate γ ) and is decreased following reinfection (at the rate θ λ ) and natural death. This gives 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 Figure 2, and the associated variables and parameters are described in Tables 1 and 2, resp.):

Figure 2

Schematic diagram of the BTB-MTB model (21).

Table 1

Description of the variables of the BTB-MTB model (21).

Variable	Interpretation
S H	Population of susceptible humans
E H1	Population of humans exposed to MTB
E H2	Population of humans exposed to BTB
I H1	Population of infected humans with clinical symptoms of MTB
I H2	Population of infected humans with clinical symptoms of BTB
R H1	Population of humans who recovered from MTB
R H2	Population of humans who recovered from BTB
S B	Population of susceptible buffalos
E B1	Population of buffalos early-exposed to BTB
E M1	Population of buffalos early-exposed to MTB
E B2	Population of buffalos at advanced-exposed BTB stage
E M2	Population of buffalos at advanced-exposed MTB stage
I BB	Population of buffalos with clinical symptoms of BTB
I MB	Population of buffalos with clinical symptoms of MTB
R BB	Population of buffalos who recovered from BTB
R MB	Population of buffalos who recovered from MTB

Table 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
η H1	Modification parameter for the reduction in infectiousness of exposed humans in comparison to humans with clinical symptoms of MTB
η H2	Modification parameter for the reduction in infectiousness of exposed humans in comparison to humans with clinical symptoms of BTB
η B1, η B2	Modification parameters for the reduction in infectiousness of exposed buffalos in comparison to buffalos with clinical symptoms of BTB
θ HH, θ BB	Modification parameters for the reduction in transmissibility of MTB to buffalos in comparison to humans
θ MM	Modification parameters for the reduction in transmissibility of BTB to humans in comparison to buffalos
γ i  (i = 1, 2)	Recovery rate of humans
γ B1, γ M1	Recovery rate of buffalos
σ i  (i = 1, 2)	Progression rate from E Hi to I Hi class
κ 1	Progression rate from E B1 to E B2 class
κ 2	Progression rate from E M1 to E M2 class
σ B2	Progression rate from E B2 to I BB class
σ M2	Progression rate from E M2 to I MB class
θ Hi  (i = 1, 2)	Exogenous reinfection rate for humans in the E Hi class
θ RB, θ RH	Exogenous reinfection rate for recovered humans
θ EB	Exogenous reinfection rate for buffalos in the exposed and recovered classes, respectively
δ H1,  δ H2	Disease-induced death rate for humans
δ B,  δ M	Disease-induced death rate for buffalos

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,

including the dynamics of early and advanced-exposed buffalos (exposed buffalo classes were not considered in the models in [8, 12, 16–18]),

allowing for BTB and MTB transmission by exposed buffalos and humans (this was not considered in [8, 12, 16–19]),

including the dynamics of humans (this was not considered in [8, 13, 18, 19]),

allowing for the reinfection of exposed and recovered buffalos and humans (this was not considered in [8, 12, 13, 16, 18]),

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 S = E = E = I = I = R = R = λ = θ = 0 in (21)), given by where, now, 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 Substituting (24) in the fifth equation in (22) shows that Similarly, by substituting (E , E ) = (0,0) into the equations for I and R in (22), it follows that Thus, the buffalo-only model reduces to the following (limited) model at steady-state:

Lemma 1 .

The following biologically feasible region of the buffalo-only model (27) is positively invariant and attracting.

Proof

Adding the equations in the buffalo-only model system (27) gives so that It follows from (30) and the Gronwall inequality that In particular, N (t) ≤ Π/μ if N (0) ≤ Π/μ . 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 data S (0) > 0, E (0) > 0, E (0) > 0, I (0) > 0, R (0) > 0. Then, the solutions S (t), E (t), E (t), I (t), and R (t) of the buffalo-only model (27) are positive for all t ≥ 0.

It is clear from the first equation of the buffalo-only model (27) that so that Using similar approach, it can be shown that E (t) > 0, E (t) > 0, I (t) > 0, and R (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 The linear stability of E 0 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 where K 1 = κ 1 + μ , K 3 = σ + μ , and K 5 = γ + μ + δ . It follows that the basic reproduction number of the buffalo-only model (27), denoted by R 0, is given by Hence, using Theorem 2 of [33], the following result is established.

Lemma 3 .

The DFE, E 0, of the buffalo-only model (27), is locally asymptotically stable (LAS) if R 0 < 1 and unstable if R 0 > 1.

The threshold quantity, R 0, 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 of R 0

The threshold quantity, R 0, 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 (E (t)), advanced-exposed buffalo with BTB (E (t)), or infected buffalo with clinical symptoms of BTB (I (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 (β η /N *) and the average duration of stay in the early-exposed class (1/K 1). Thus, the average number of BTB infections generated by early-exposed buffalos is given by

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 (β η /N *), the probability that early-exposed buffalo survived the early-exposed class and move to the advanced-exposed class (κ 1/K 1), and the average duration of stay in the advanced-exposed class (1/K 3). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by

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 (β /N *), the probability that an advanced-exposed buffalo survived the advanced-exposed class and move to the symptomatic class I (κ 1 σ /K 1 K 3), and the average duration of stay in the symptomatic class I (1/K 5). Thus, the average number of BTB infections generated by advanced-exposed buffalos is given by

The sum of the terms in (37), (38), and (39) gives R 0. That is, the average number of new infections generated by infected buffalos (early-exposed, advanced-exposed, or symptomatic) is given by (noting that S * = Π/μ and N * = Π/μ )

The epidemiological implication of Lemma 3 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 (E 0). 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 R 0 < 1 [12, 17, 25, 36]). The epidemiological implication of this phenomenon is that the classical requirement of R 0 < 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 R 0 = 1 whenever the bifurcation coefficient, a, given by (A.9) (in Appendix A), is positive.

The proof of Theorem 4, 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 θ = θ = 0), the backward bifurcation coefficient, a, given by (A.9) in Appendix A, reduces to (it should be noted, from Appendix A, that A 1 > 0, A 2 > 0, β * > 0, and, from Theorem 2, that all parameters of the buffalo-only model (22) are nonnegative): where K 2 = κ 2 + μ , K 4 = σ + μ and K 6 = γ + μ + δ . 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 dynamics Mycobacterium tuberculosis in human populations). This result is summarized below.

Lemma 5 .

The buffalo-only model (22) does not undergo backward bifurcation at R 0 = 1 in the absence of the reinfection of exposed and recovered buffalos (θ = θ = 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 (E 0) of the buffalo-only model (27) in the absence of reinfection (i.e., θ = θ = 0) below.

3.2. Global Asymptotic Stability of the DFE

Consider the buffalo-only model (27) in the absence of the reinfection of exposed (θ = 0) and recovered (θ = 0) buffalos.

Theorem 6 .

The DFE, E 0, of the buffalo-only model (27) with θ = θ = 0 is globally asymptotically stable (GAS) in Γ if R 0 ≤ 1.

Consider the buffalo-only model (27) in the absence of reinfection (θ = θ = 0). Furthermore, let R 0 ≤ 1. Consider the linear Lyapunov function F = a 0 E + a 1 E + a 2 I , where with Lyapunov derivative given by (where a dot represents differentiation with respect to time t) Since all the parameters and variables of model (27) are nonnegative (Theorem 2), it follows that for R 0 ≤ 1 with if and only if E = E = I = 0. Thus, it follows, by LaSalle's Invariance Principle [37], that Since lim⁡⁡sup⁡I (t) = 0 (from (44)), it follows that, for sufficiently small ϖ* > 0, there exists a constant M > 0, such that lim⁡⁡sup⁡I (t) ≤ ϖ* for all t > M. Hence, it follows from the fifth equation of the buffalo-only model (27) that, for . Thus, by comparison theorem [38], R = lim⁡⁡sup⁡R ≤ γ ϖ*/μ , so that, by letting, ϖ* → 0, Similarly, it can be shown that Thus, it follows from (45) and (46), that R ≥ 0 ≥ R . Hence, Furthermore, substituting (44) in the first equation of (27) shows that Thus, by combining equations (44), (47), and (48), it follows that every solution of the equations of the buffalo-only model (27), with θ = θ = 0 and initial conditions in Γ, approaches E 0 as t → ∞ (whenever R 0 ≤ 1).

Theorem 6 shows that, in the absence of the reinfection of exposed and recovered buffalos (i.e., θ = θ = 0), BTB can be eliminated from the buffalo-only population if the reproduction number of the model (R 0) 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 E 0 when R 0 < 1 (in line with Theorem 6).

Figure 4

Simulations of the buffalo-only model (27), showing the total number of infected buffalos with clinical symptoms of BTB (I (t)(t)) at time t as a function of time. Parameter values used are as given in Table 3 with (a) β = 0.00733 (so that, R 0 = 0.7036 < 1) and (b) β = 0.0733, δ = 0 (so that, ).

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., θ = θ = 0). Solving the equations of the buffalo-only model (27) at steady-state gives the following general form of the endemic equilibrium (denoted by E 1): where with the force of infection at steady-state (λ **) given by

Using (50) in the expression for λ ** in (51) shows that the nonzero equilibrium of the model (22) satisfies the linear equation: where b 1 = K 5 μ (K 3 + κ 1) + σ κ 1(μ + γ ) and b 2 = K 1 K 3 K 5 μ (1 − R 0). Clearly, the coefficient b 1 is always positive, and b 2 is positive (negative) if R 0 is less than (greater than) unity, respectively. Thus, the linear system (52) has a unique positive solution, given by λ ** = −b 2/b 1, whenever R 0 > 1. Further, the force of infection for buffalos (λ **) is negative whenever R 0 < 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 θ = θ = 0, has a unique endemic equilibrium, E 1, whenever R 0 > 1, and no endemic equilibrium otherwise.

3.3.1. Global Asymptotic Stability of Endemic Equilibrium

The global asymptotic stability of the unique endemic equilibrium (E 1) of the buffalo-only model is explored for the special case without reinfection (θ = θ = 0) and BTB-induced death in buffalos (δ = 0). It is convenient to define the stable manifold of the DFE (E 0) of the buffalo-only model (27).

Theorem 8 .

The unique endemic equilibrium (E 1) of the buffalo-only model (27), with θ = θ = δ = 0, is GAS in Γ∖Γ1 if .

Consider the buffalo-only model (27) with θ = θ = δ = 0. For this case, it follows from Theorem 7 that the buffalo-only model (27) has a unique endemic equilibrium whenever . Furthermore, setting δ = 0 in model (27) shows that N (t) → Π/μ as t → ∞. Consider the following nonlinear Lyapunov function (of Goh-Volterra type) for the subsystem of model (27) involving the state variables S , E , E , and I (noting that N (t) is now replaced by its limiting value Π/μ ): where . The Lyapunov derivative of F is given by Using the following steady-state relations (obtained from (27)), the Lyapunov derivative can be simplified to

Thus, Finally, since the arithmetic mean exceeds the geometric mean, it follows then that Furthermore, since all the model parameters are nonnegative, it follows that for . Thus, is a Lyapunov function for the subsystem of model (27) on Γ∖Γ1. Therefore, it follows, by LaSalle's Invariance Principle [37], that Since I (t) → I ** as t → ∞, it follows from the equation for dR /dt in (27) that R (t) → γ I **/μ = R ** as t → ∞. The proof is concluded using similar arguments as in the proof of Theorem 6.

The epidemiological implication of Theorem 8 is that BTB will be endemic in the buffalo population if (and θ = θ = δ = 0). Figure 4(b) depicts the solutions of model (27) for the case when and θ = θ = δ = 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 θ = θ = 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]): (S (0), E (0), E (0), E (0), E (0), I (0), I (0), R (0), R (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]
η H1	[0, 1)	0.5	Fitted
η H2	[0, 1)	0.5	Fitted
η B1	(0.4455, 0.5045)	0.45	Fitted
η B2	(0.495, 0.605)	0.55	Fitted
θ BB	[0, 1)	0.5	Fitted
θ MM	[0, 1]	0.5	Assumed
θ HH	[0, 1)	0.5	Assumed
γ i  (i = 1, 2)	(0.0000823, 0.000823)	0.000453	[14, 23]
γ B1	(0.00774, 0.00946)	0.0086	[10]
γ M1	(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]
σ B2	(0.25, 0.35)	0.3	[27]
σ M2	(0.36, 0.44)	0.4	[27]
θ Hi  (i = 1, 2)	[0, 0.1]	0.00271	[23]
θ RB, θ RH	[0.002439, 0.002981]	0.00271	[23]
θ EB	[0.002439, 0.002981]	0.00271	[23]
δ H1, δ H2	[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]

Figure 5 depicts a box plot of R 0, as a function of the number of LHS runs carried out (N = 1000), from which it is evident that the distribution of R 0 lies in the range R 0 ∈ [0.25, 0.75] (each box plot displays the upper and lower quartile ranges of R 0. a horizontal line within the box is the median value, and values of R 0 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 R 0 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 R 0 (hence the disease dynamics) are the BTB transmission rate (β ), the recovery rate of buffalos (γ ), and the BTB-induced mortality in buffalos (δ ).

Figure 5

Box plot of R 0 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 θ = θ = 0.

Figure 6

PRCC values of the parameters of the buffalo-only model (22), using R 0 as the output function. Parameter values used are as given in Table 3.

Similarly, Figure 7 depicts the box plot of the buffalo-only model (22) using total number of symptomatic buffalos (I + I ) 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 (Π), the natural (μ ), and the disease-induced (δ ) 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 (I + I ) 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 θ = θ = 0.

Figure 8

PRCC values of the parameters of the buffalo-only model (22), using total number of symptomatic buffalos (I + I ) 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 Section 3, that the following biologically feasible region 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 θ = θ = 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 The associated next generation matrices of the BTB-MTB model (21), denoted by F and V , are given, respectively, by

It follows then that the reproduction number of the BTB-MTB model (21), denoted by R , is given by where R and R are the associated reproduction numbers for humans infected with MTB and with BTB, respectively, given by where Q 1 = σ 1 + μ , Q 2 = σ 2 + μ , Q 3 = γ 1 + μ + δ , and Q 4 = γ 2 + μ + δ , and R 0 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, E 0, of model (21), with θ = θ = 0, is LAS in Ω if R < 1, and unstable if R > 1.

It can be shown, as in Section 3, 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 (θ and θ ) are set to zero. This is due to the reinfection of exposed and recovered humans (i.e., θ ≠ 0 and θ ≠ 0). To illustrate this fact, it is shown that the DFE (E 0) 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 (R ) is less than unity.

4.2. Global Asymptotic Stability of DFE

Theorem 10 .

The DFE, E 0, of the BTB-MTB model (21) with θ = θ = θ = θ = θ = θ = θ = θ = 0, is GAS in Ω if R < 1.

The proof of Theorem 10, 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 θ , for the reduced likelihood of humans acquiring BTB infection from buffalos):

no transmission of BTB from buffalos to humans: θ = 0,

low rate of transmission of BTB from buffalos to humans: θ = 0.25,

moderate rate of transmission of BTB from buffalos to humans: θ = 0.50,

high rate of transmission of BTB from buffalos to humans: θ = 0.75.

The simulation results obtained, depicted in Figure 9(a), show that the cumulative number of new MTB cases in humans decreases with increasing rate of BTB transmission to humans by buffalos (θ ).

Figure 9

Cumulative number of new cases of (a) MTB infection in humans and (b) BTB infection in buffalos. Parameter values used are as given in Table 3, with various values of θ (a) or θ (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 (θ ) are considered, namely, none (θ = 0), low (θ = 0.25), moderate (θ = 0.50), and high (θ = 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.

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.

In the absence of the reinfection of exposed and recovered buffalos (θ = θ = 0), the buffalo-only model is shown to have unique endemic equilibrium whenever its reproduction number exceeds unity (R 0 > 1). This equilibrium is shown to be globally asymptotically stable for the special case where the disease-induced mortality in buffalos is negligible (δ = 0).

Detailed uncertainty analyses of the buffalo-only model, using a reasonable set of parameter values and ranges (Table 3) 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 (R 0) 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.

The BTB-MTB model also undergoes backward bifurcation. Unlike in the buffalo-only model, this phenomenon persists even if the bovine-associated reinfection terms (θ and θ ) are set to zero. This is due to the reinfection of exposed and recovered humans (θ ≠ 0, θ ≠ 0, θ ≠ 0, θ ≠ 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.

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).

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.

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	230