Journal of Veterinary Medicine and Surgery Open Access

Mathematical Modelling of Bovine Tuberculosis Transmission Dynamics: Role of Combination of Control Measure

Sylas Oswald^1*, Theresia Crispin Marijani¹, Goodluck Mika Mlay¹ and Winifrida Benedict Kidima^{2 1}Department of Mathematics, University of Dar-es-Salaam, P.O.Box 35062, Dar-es-salaam, Tanzania
²Department of Zoology and Wildlife Conservation, University of Dar-es-Salaam, P.O.Box 35064, Dar-es-salaam, Tanzania
^*Correspondence: Sylas Oswald, Department of Mathematics, University of Dar-es-Salaam, P.O.Box 35062, Dar-es-salaam, Tanzania, Email:

Received: 12-Jun-2024, Manuscript No. IPJVMS-24-20424; Editor assigned: 14-Jun-2024, Pre QC No. IPJVMS-24-20424 (PQ); Reviewed: 28-Jun-2024, QC No. IPJVMS-24-20424; Revised: 12-Mar-2025, Manuscript No. IPJVMS-24-20424 (R); Published: 19-Mar-2025, DOI: 10.36648/2574-2868.9.1.41

Introduction

Bovine Tuberculosis (BTB) is a contagious and life-threatening infectious disease of cattle. The bacterium Mycobacterium bovis (M. bovis) that causes disease, can also afflict many other mammals such as humans, goats, pigs, cats and dogs. BTB remains a significant challenge in both human and animal health worldwide. It poses a substantial economic burden, particularly in regions heavily reliant on live-stock agriculture. In 1993, the disease was declared a global state of emergency. Despite our knowledge of how to effectively prevent and cure BTB through half a century of development and progress, more than 1.6 million people have still died from it. In 2014, BTB claimed the lives of 1.5 million people, including 890,000 men, 480,000 women, and 140,000 children. India, Indonesia, and China had the largest number of cases, accounting for 23%, 10%, and 10% of the global total, respectively. Traditional control measures, such as culling infected animals and movement restrictions, have had limited success in curbing the transmission of BTB. To address this complex issue, researchers have increasingly turned to mathematical modelling to gain insights into the dynamics of BTB transmission and assess the potential impact of various control strategies. This study delves into the mathematical modelling of BTB transmission dynamics, with a novel focus on the synergistic role of a combination of vaccination, treatment, and education campaigns as control measures [1].

The transmission route entails direct, close contact through inhaling sputum droplets (thick mucus produced in the lungs) contaminated with M. bovis bacteria exhaled by infected animals. On the other hand, indirect transmission of the disease can occur when one comes into contact with an infected animal or ingests material heavily contaminated with M. bovis, such as sputum, pus, urine, feces, and other excrement of infected animals. There are inherent practical difficulties associated with detecting M. bovis in environmental samples. However, the bacteria have been traced in diverse sources, such as soil, excrement, hay, and pasture.

Bovine Tuberculosis: A Global Challenge

Bovine tuberculosis poses a substantial threat to both animal and human populations. In cattle, it results in reduced milk and meat production, increased slaughterhouse condemnation rates, and trade restrictions, contributing significantly to economic losses in affected regions. In humans, zoonotic transmission of BTB remains a concern, particularly in areas where close contact with infected cattle occurs. Given its impact on both public health and the livestock industry, BTB demands innovative and effective control measures [2-5].

Mathematical Modelling in Infectious Disease Control

Mathematical modelling has proven to be a powerful tool in understanding and predicting the dynamics of infectious diseases. In the context of BTB, previous studies have employed compartmental models, network models, and agent-based models to simulate disease spread and evaluate control strategies. These models have provided valuable insights into the role of factors such as cattle demographics, wildlife reservoirs, and testing protocols in BTB transmission dynamics.

Individual Control Measures

Many studies have assessed the efficacy of individual control measures, such as vaccination and antimi-crobial treatment, in reducing BTB prevalence. Vaccination with Bacillus Calmette-Gu´erin (BCG) has shown promise in reducing the severity of BTB in cattle. Antimicrobial treatment can help clear infections but is often logistically challenging due to the lengthy treatment duration and concerns about antibiotic resistance.

Materials and Methods

Education Campaigns as a Control Measure

Education campaigns aimed at promoting best practices in cattle management, biosecurity measures, and early disease detection have been recognized as crucial components of BTB control. These campaigns empower farmers and stakeholders with knowledge to reduce the risk of BTB transmission within and between herds.

In Africa and Europe, bovine tuberculosis is a serious threat to the economy as well as human and animal health. For example, presented a deterministic model for examining the impact of separating two large human populations based on the probability of coming into contact with cattle for the prevention of bovine tuberculosis. The model encompasses three types of incidents: Those occurring between cattle, between humans, and between cattle and humans. The results of the study, based on the stability of the disease-free equilibrium and sensitivity analysis of the model’s parameters, showed that quarantine measures for cattle and the parameter related to medical masks significantly contributed to reducing the basic reproduction number and, consequently, decreasing the disease transmission rate. The effects of vaccines on cattle, treatments, public health education campaigns for humans, and communal water source sharing, in relation to disease transmission, were not fully considered.

The study by also presented a dynamic mathematical model for the transmission of bovine BTB. The results from the model simulations supported vaccine administration as the best strategy for reducing BTB infections. Furthermore, they suggested that coupling a public health education campaign with treatment for humans could maintain R₀ below one, significantly reducing the spread of bovine tuberculosis. However, their study did not address the impact of treatment as a measure for controlling the transmission dynamics of bovine tuberculosis.

A study by revealed that the use of slaughter and quarantine methods to reduce the number of infectious cattle was found to be the most important control measures to minimize the prevalence of BTB disease.

The study done by revealed that a relatively high proportion of BTB infections in the Ngorongoro district are due to the husbandry practices of a semi-nomadic system where the Maasai search for water and pasture during the dry season. It is recommended that implementing an education campaign for cattle owners would be the best solution.

In developing countries such as Tanzania, Bovine Tuberculosis (BTB) remains a serious threat to human health because the disease is still largely hidden in society. Reported cases of BTB in Tanzania are concentrated in specific regions, including the Northern part of the country (Arusha, Kilimanjaro, and Manyara), dairy farms in Kibaha, and some areas in Morogoro districts. The prevalence of the disease varies from one region to another, depending on the concentration of cattle herds in a particular place, ranging from 0.2 percent to 13.3 percent.

Novelty of the Study

While previous research has explored individual control measures in isolation, this study stands out for its innovative approach in combining vaccination, treatment, and education campaigns within a mathematical modelling framework. Investigating the synergistic effects of these measures can provide a more comprehensive understanding of their potential impact on BTB transmission dynamics. Moreover, it addresses the practical question of how these measures can be integrated into a cohesive control strategy to reduce both the prevalence of BTB and its economic burden [6-10].

Model Formulation

The model under consideration here entails two populations: humans and cattle, which often come into contact. The management system for grazing is a ranching system in which only cattle are confined. The model operates on the assumption that at any time t, the human population, denoted by N_h(t), is divided into five classes: Susceptible noneducated (S_nh(t)), Susceptible educated (S_eh(t)), Exposed (E_h(t)), Infectious (I_h(t)), and Recovered (R_h(t)) individuals.

The total human population N_h(t) is given by N_h(t)=S_nh(t)+ S_eh(t)+E_h(t)+I_h(t)+R_h(t).

The recruitment of humans into the susceptible non-educated class occurs at a constant rate π_h. The assumption is that the education strategy is executed at a rate of ψ only for susceptible, non-educated humans to reduce the disease’s transmissibility. Moreover, the study assumes that the education given to a targeted group does not necessarily guarantee lifelong protection. Susceptible non-educated humans can acquire infection through the consumption of cattle products and aerosols from infected cattle, as well as through inter-human transmission at the rate φ1, and then move to the exposed class. The variable φ1 is the force of infection given by

where a₁ and a₂ are probabilities that infectious cattle and human infects susceptible non educated humans, respectively and θ ∈ [0, 1] is the efficacy of the education campaign that is being implemented. Furthermore, an educated human may contract the disease by consuming contaminated cattle products and inhaling aerosols from infected cattle at the rate φ2, leading to a transition to the exposed class. The variable φ2 represents the force of infection for susceptible educated humans.

where a is the probability that infectious cattle infects susceptible educated human. If no education campaign is extended to susceptible non-educated human, then θ=0. Both educated and non-educated susceptible humans may leave their respective classes following the occurrence of natural death at a rate of μ_h. Some of the exposed individuals leave the compartment upon gaining full recovery, denoted as ρh, due to treatment and a natural or disease-induced death rate of μ_h and α_h, respectively. Considering the transmission dynamics of BTB disease, the study assumes that the treatment does not guarantee permanent protection. As such, recovered humans may either progress to the exposed class at a rate of τ_h or leave the compartment following the occurrence of natural death at a rate of μ_h. This study presupposes that at any time t, the cattle population denoted by N_c(t) is divided into four classes, of Susceptible nonvaccinated S_nv(t), Susceptible vaccinated S_v(t), Exposed E_c(t) and Infectious I_c(t) cattle. The total cattle population N_c(t) is given by N_c(t)=S_nc(t)+S_v(t)+E_c(t)+I_c(t). At any given time t, it is also presumed that cattle are recruited into the nonvaccinated cattle group at a constant rate π_c. Healthier cattle are vaccinated at a rate ω to reduce the transmission of M. bovis from the contaminated environment and inter-cattle transmission. Susceptible non-vaccinated cattle can contract infections from their respective environment, inter-cattle transmission, and infection from infected humans. The progression to the exposed class occurs at a rate σ₁. The variable σ₁ represents the force of infection for nonvaccinated cattle.

Where ε₁ is the probability of an infected human infecting susceptible, non-vaccinated cattle. Moreover, the presumption is that it is difficult to differentiate infected cattle from a contaminated environment or from inter-cattle transmission. ε₂ is the rate at which cattle get infected through nose-to-nose contact, aerosol inhalation, grazing areas, and all other environmental elements that can infect cattle. λ ∈ [0, 1] represents the vaccine efficacy. If no susceptible cattle are vaccinated, then λ=0. Furthermore, vaccinated cattle may contract infection from both a contaminated environment and inter-cattle transmission at the rate σ₂, leading them to move to the exposed class. The variable σ₂ constitutes the force of infection for susceptible vaccinated cattle, as given by 

Where ε represents the probability that an infected human will transmit the disease to susceptible vaccinated cattle. Both vaccinated and non-vaccinated susceptibles may leave their respective classes due to natural death at a rate of μ_c. Exposed cattle also exit this class, either because of natural death 5 at a rate of μ_c or due to progressive incubation into the infectious class at a rate of γ_c. Furthermore, infected cattle exit the class through natural death at a rate of μ_c and disease-induced death at a rate of αc (Figure 1).

Model Diagram

Figure 1: Schematic flow diagram for the dynamics of bovine tuberculosis disease in the presence of some intervention strategies for model.

Model Equations

The biological description and schematic flow diagram in Figure 1 together result in the following system of nine nonlinear ordinary differential equations (Table 1).

Table 1: Parameters and their descriptions.

Basic Properties of the Model

Positivity of the solutions: To ensure that the model (1) is epidemiologically meaningful and well-posed, it is necessary to demonstrate that all state variables are non-negative ∀t ≥ 0.

Lemma 1. Let {S_nh(0) ≥ 0, S_eh(0) ≥ 0, E_h(0) ≥ 0, I_h(0) ≥ 0, R_h(0) ≥ 0, S_nv(0) ≥ 0, S_v(0) ≥ 0, E_c(0) ≥ 0 and I_c(0) ≥ 0 } of themodel system (1) are satisfied then the solutions {S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c } are non-negative for all t ≥ 0 in Ω.

Proof. From the first equation of system (1),

In the absence of human recruitment πh it follows that

Upon integration gives,

In the same way, the remaining state variables give;

Therefore, the solution set {S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c} of the model system (1) is non negative for all t ≥ 0 in Ω.

Invariant Region

Lemma 2. The feasible region defined by Ω=Ωh ∪ Ωc ∈ R⁵ + × R⁴ + where Ω_h={S_nh, S_eh, E_h, I_h, R_h ∈ R⁵ +: S_nh+S_eh+E_h+I_h+R_h=N_h ≤ π_h/μ_h} and Ω_c = {S_nv, S_v, E_c, I_c ∈ R⁴ +:S_nv+S_v+E_c+I_c=N_c ≤ π_c/μ_c}, is positively invariant and attracting with regard to the model system (1).

Proof. Given S_nh (0) ≥ 0, S_eh (0) ≥ 0, E_h (0) ≥ 0, I_h (0) ≥ 0, R_h (0) ≥ 0, S_nv (0) ≥ 0, vs. (0) ≥ 0, E_c (0) ≥ 0 and I_c (0) ≥ 0, it is sufficient to prove that {S_nh(t), S_eh(t), E_h(t), I_h(t), R_h(t), S_nv(t), S_v(t), E_c(t), I_c(t)} ∈ R⁹+ is bounded.

Case 1: By summing equations for human population from model system gives;

Separating variables and applying anti-derivative of an integrating factor it gives;

The same procedure can be used to prove that;

Since the total population N_h(t) as well as N_c(t) is positive for all t ≥ 0. It is well defined that,

This proves that all solutions of the BTB model system (1) with initial conditions in Ω for all t>0.

Model Analysis

Disease free equilibrium: Disease-Free Equilibrium (DFE) of the mathematical model system (1) is the point where there is no disease. Disease-free equilibrium point is obtained by setting E_h=I_h=R_h=E_c=I_c=0 in all equations. Now, Let E₀ be disease-free equilibrium point. Therefore, the Disease Free Equilibrium Point (DFE) denoted by E₀ can be expressed as

The Effective Reproduction Number

The Effective Reproduction Number (R_e) is the expected number of secondary cases produced by a single (typical) infection in a completely susceptible population and in nonsusceptible hosts [6,7,10]. Effective reproduction number (R_e) is the the spectral radius of the next-generation matrix, denoted 10 by R_e=ρ(FV⁻¹) at DFE with F and V, respectively given by

It follows that the effective reproduction number of the model system (1), is given by

such that,

So, when there is no intervention strategy such that ψ=0 and ω=0,” which refers to the education campaign and vaccination rate,” and without any treatment for individual humans being implemented (ρ_h=0), then the basic reproduction number R0 will be

Biologically, it is meaningful to compare the threshold values: the basic reproduction number without a vaccine, denoted as R₀, and the effective reproduction numbers with and without a vaccine, denoted as R_v and R_e respectively, such that R_ev0. The basic reproduction number, often denoted as R₀, plays a crucial role in epidemiological studies by facilitating predictions of future infections of interest. When the basic reproduction number is less than one (R₀<1), it implies that, on average, an infectious individual leads to the infection of fewer than one other individual. Consequently, over time, the bTB disease could naturally die out, leading to a population free from BTB. In other words, both the human and cattle populations would remain free from the disease invasion if R₀<1 [11-15]

Conversely, if the basic reproduction number (R₀) is greater than one, it implies that each infected individual, on average, will lead to more than one newly infected individual. In this scenario, the infection persists, causing the disease to continue spreading in the population. Therefore, if R₀>1, BTB will continue to spread

Local Stability of Disease-Free Equilibrium (DFE)

The local stability analysis of DFE enables us to understand how a system behaves near an equilibrium point, but not necessarily at a specific equilibrium point. In other words, local stability investigates the nature of the system in the vicinity of the equilibrium point.

Theorem 3.1. The BTB model system at DFE (E₀) is locally asymptotically stable if R_e<1 and unstable if R_e>1.

Proof. The Jacobian matrix of the model system at disease free equilibrium point (E₀) is given by

From the Jacobian matrix J_E0, the first, fifth and seventh columns contains diagonal entries. Therefore, the diagonals −μ_h, −(μ_h + τ_h), and −μ_c are the three eigenvalues of the Jacobian J_E0. Thus, excluding these columns and the corresponding rows, the remaining eigenvalues are computed. Then the reduced 6 × 6 matrix from J_E0 becomes.

From the Jacobian matrix ξ, the first and fourth columns contains diagonal entries. Therefore, the diagonals − (μ_h + ψ) and − (μ_c + ω) are the two eigenvalues of the Jacobian ξ. Thus, excluding these columns and the corresponding rows, the remaining eigenvalues are computed. Then the reduced 4 × 4 matrix from ξ bec becomes.

Matrices G₁, G₂, G₃ and G₄ are defined as follows

Apparently, matrices G₃ and G₄ are singular matrices with | G₃| = 0 and |G₄| = 0. Now, it is sufficient to demonstrate that the disease-free equilibrium point of the model system (1) is stable only when the trace and determinant of matrices G₁ and G₂ are negative and positive, respectively. Trace and determinant of matrix G₁ denoted by Tr(G₁) and det(G₁) are respectively given by.

Trace and determinant of matrix G₂ denoted by Tr(G₂) and det(G₂) are respectively given by

Thus det(G₁)> 0 and det(G₂)> 0 if and only if R_e<1. Since the traces of matrix G₁ and G₂ are negative and the determinants are strictly greater than zero when R_e<1, then the disease-free equilibrium point E₀ is locally asymptotically stable when R_e<1 and unstable otherwise. Epidemiologically, it implies that the bTB disease can be eliminated in the endemic area when R_e<1, particularly when the vaccination rate (ω) is kept constant without limitation, and the education campaign rate (ψ) is also executed for every individual. Conversely, if R_e>1, then each infectious individual produces more than one new infected individual, implying that the disease can spread rapidly in the entire population [16-20].

Global Stability of Disease-Free Equilibrium (DFE)

This subsection analyses the global stability pertaining to DFE aimed to establish the asymptotic behaviour of the model system (1) beyond just the neighbourhood points of the model disease-free steady state.

Lemma 3. The disease-free equilibrium point is globally asymptotically stable if R_e<1 and unstable if R_e>1.

Proof. The approach of [6] is used to analyse the global stability of DFE of the model system (1). Using this approach, the model system (1) is written as follows.

where ζ_n stands for classes that do not transmit bTB disease, that is ζ_n=(S_nh, S_eh, R_h, S_nv, S_v) T and ζ_i stands for classes that can transmit bTB disease, that is ζ_i=(I_h, E_h, E_c, I_c) T. Here, T stands for the transposition of ζ_n and ζ_i, and also ζ_E0 is ζ_n at D_FE (E₀). The matrix M₁ is obtained by differentiating the nontransmitting equations of the model system (1) with respect to non-transmitting variables at E₀, whereas the matrix M₂ is obtained by differentiating the non-transmitting equations of the model system (1) with respect to the transmitting variables. The disease-free equilibrium point (E₀) is globally asymptotically stable if the eigenvalues of M₁ are real and negative, and if M₃ is a Metzler stable matrix with nonnegative off-diagonal elements. Therefore, the model system (1) yields the following.

Also,

Since the matrix M₁ is a lower triangular matrix, then the eigenvalues will be the diagonal entries presented as follows: ξ₁=−(μ_h+ψ), ξ₂=−μ_h, ξ₃=−(μ_c+ω) and ξ₄=−μ_c which are all negative and real. This shows that at the DFE the system dζ_n/d_t=M₁(ζ_n−ζ_E0)+M₂ζ_i is globally asymptotically stable. Furthermore, the matrix M₃ is obtained by differentiating the transmitting equations of the model system (1) with respect to the transmitting variables which gives

Testing whether the matrix M₃ is a Metzler stable matrix requires applying the approach as propounded by [16,11]. The following Lemma is used.

Lemma 4. Let G be a square Metzler matrix stable which can be written in block form

where A and D are square matrices. Then, G is a Metzler stable matrix if and only if matrix A and D−CA⁻¹B or matrix D and A−BD⁻¹C are Metzler stable.

By comparing the matrix M₃ and a square Metzler matrix G, then the matrices A, B, C and D are expressed as follows:

After some computations and simplifications,

A − BD⁻¹C is said to be Metzler stable matrix if and only if (μ_h +γ_h) (μ_h+α_h+ρ_h) ≥ 0. As already noted, it can be seen that matrix M₁ has all the eigenvalues that are real and negative and matrix M₃ is the Metzler matrix as its off-diagonal elements are so non-negative, that |A−BD⁻¹C| ≥ 0 =⇒ M₃(i,j) ≥ 0 for all indices i/=j. Implicitly, that the disease-free equilibrium point is globally asymptotically stable when R_e<1, otherwise it is unstable. In biological terms, the bTB disease will eventually die out in the population given R_e<1 no matter how huge the disease invaded. Otherwise, if R_e>1 the disease will spread even more as there will be many new infections in the population.

Endemic Equilibrium Point (EEP)

Existence of the equilibrium solutions: The Endemic Equilibrium (EE) denoted by E∗ such that E∗=(S∗_nh, S∗_eh, E∗_h, I∗_h, R∗_h, S∗_nv, S∗_v, E∗_c, I∗_c) is the state where the population is not free from the infection. At this point the disease persists in the population. In other words, the disease cannot be eradicated from the population. In the case of bTB, the number of infectious cases is not equal to zero such that, Eh/=0, Ih/=0, Rh/=0, Ec/=0, Ic/=0. Setting each equation in the model system (1) equal to zero, then the steady state of the system becomes.

Finally, when solving for the corresponding variable from equation (11) the Endemic Equilibrium denoted by E∗=(S∗_nh, S∗_eh, E∗_h, I∗_h, R∗_h, S∗_nv, S∗_v, E∗_c, I∗_c) is solved and given by

where,

From equation (12), the value of I∗_h and I∗_c is substituted into

and after some algebraic computations and simplifications, a polynomial of degree three results A_σ1 ∗³+B_σ1 ∗²+C_σ1 ∗=0. Upon algebraic simplification, A_σ1 ∗³+B_σ1 ∗²+C_σ1 ∗=0 results to σ1 ∗∗A_σ1 ∗²+Bσ∗ ¹+C=0, where σ1 ∗=0 or A_σ1 ∗²+Bσ∗ ¹+C=0

corresponds to Disease Free Equilibrium (DFE), whereas A_σ1 ^∗2 +Bσ^∗1+C=0, can also be written in the form:

which satisfies Endemic Equilibrium. The value of A is strictly positive. Depending on the signs of B and C, there are three cases to consider as having positive real root of the force of infection (σ^∗₁) as follows:

Case 1: If B<0 then model system (1) has a stable endemic equilibrium point when C<0. This equilibrium happens when R_e>1 as interpreted from (13). In this case backward bifurcation is not possible due to the absence of multiple equilibria.

Case 2: Exactly one endemic equilibrium point. Suppose B<0 and C=0 or B²−4AC=0. In other words, the polynomial Aσ1 ∗2 +Bσ∗ 1+C=0 has just one positive root and hence the model system (1) has unique endemic equilibrium point.

Case 3: Two endemic equilibria. If B<0, C>0 and B2−4AC>0, then the polynomial A_σ1 ∗²+Bσ∗ ¹+C=0 has two positive real roots. In other words, the model system (1) has two endemic equilibria and hence there is a possibility of backward bifurcation. These three cases are summarized under theorem 3.2.

Theorem 3.2. The number of positive endemic equilibria of bovine tuberculosis model (1) is hereunder summarised thus:

If C<0, R_e>1, then the system has a unique endemic equilibrium.

If B<0 and C=0 or B2⁻⁴AC=0, then the system has exactly one endemic equilibrium.

If B<0, C>0 and B2⁻⁴AC>0, then the system has exactly two endemic equilibria.

Otherwise there are no endemic equilibria, i.e., when AC>0 and B>0.

Stability of Endemic Equilibrium Point (EEP) of Model with Interventions

In epidemiological models, forward or backward bifurcation has vital implications for the biological control measures of infectious diseases. Bifurcation analysis of the equilibrium points reveals whether the disease is completely reducible or persistent in the afflicted population. This analysis occurs at the disease-free equilibrium using Center Manifold theory, as presented in [6]. It is done by renaming the state variables of the model system as follows.

Resulting system can be written in the form of

Let ε^∗₂ be the bifurcation parameter, then, the system is a system of model equations at disease free equilibrium point when ε₂=ε^∗₂ with R₀=1. Hence, solving for ε^∗₂ from R₀=1 in

Next, the system is transformed with ε₂=ε^∗₂, which possesses a simple zero eigenvalue. Center Manifold Theory is employed to analyze the dynamics of in the vicinity of ε₂=ε^∗₂. Consequently, the Jacobian matrix of the system at the disease-free equilibrium, denoted as J (ε^∗₂), is given by

Now, the right and left eigenvectors associated with the zero eigenvalues are calculated. The right eigenvector associated with the zero eigenvalue is given by w=(w₁, w₂, ..., w₉)T, which results in the following equations:

Solving the system (15) for w_i′_s, i=1, 2, ..., 9, gives the following right eigenvectors

furthermore, to calculate left eigenvector given by v=(v₁, v₂, ..., v₉)^T, which satisfy v. w=1, the matrix J(ε^∗₂) is transposed and becomes.

Then solving J(ε^∗₂)^T, the following are obtained

Solving the system (16) for v_i'_s, i = 1, 2, ..., 9, yields the following left eigenvectors

From [6], Theorem 4:1 is used to establish the conditions for the existence of forward or backward bifurcations of the endemic equilibrium point near R₀=1.

Lemma 5. Consider the following general system of ordinary differential equations with a parameter β such that

such that f (0, β) ≡ 0, where 0 is an equilibrium point of the system with the following conditions:

is the linearization matrix of the model system (1) around the equilibrium 0 with β evaluated at 0.

Zero is a simple eigenvalue of A, and all other eigenvalues of A have negative real parts.

Matrix A has a non-negative right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue.

The sign of a and b always determines the local dynamics of the system around the equilibrium point.

a>0, b>0. When β <0 with |β| << 1; 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0< β << 1; 0 is unstable and there exists a negative, locally asymptotically stable equilibrium.

a<0, b<0. When β <0 with |β| << 1; 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0< β << 1; 0 is asymptotically stable equilibrium, and the re exists a positive unstable equilibrium.

a<0, b<0.W hen β <0 with| β| << 1; 0 isu nstable, anda positive unstable equilibrium appears.

a<0, b>0. When β changes from positive to negative, 0 changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable. Particularly, if a>0 and b>0, then, a subcritical (or backward) bifurcation occurs at β =0.

Computation of bifurcation coefficients to determine the local dynamics of the transformed system (14), the value of an and b are computed to probe whether the model system (14) shows forward or backward bifurcation. Since v₁=v₂=v₅=v₆=v₇=0 for k=1; 2; 5; 6; 7 then k=3; 4; 8; 9 are considered. For the system (14), the associated non-zero second order partial derivatives at disease-free equilibrium and at ε₂= ε^∗₂ are given by

For the value of b, it can be shown that there exists a nonvanishing partial derivative.

Therefore, a<0 and b>0, hence the model system (1) exhibit a forward bifurcation.

Lemma 6. The unique endemic equilibrium is guaranteed and by 3.2, the model is locally asymptotically stable for R₀>1. In addition, by 3.2 item (i), the model undergoes backward bifurcation when a>0. This holds only if v₃<0 and v₈<0, otherwise it undergoes forward bifurcation.

Results and Discussion

Bifurcation Diagram

Figure 2 illustrates the forward bifurcation (transcritical bifurcation) at R₀=1 when a<0 and b>0 for a model system described by equation (1). Thus, the bifurcation diagram for the model system shows that the disease-free equilibrium and endemic equilibria exchange stability when R₀=1. This scenario indicates biologically that the model system (1) is globally asymptotically stable at the disease-free level. The equilibrium point emerges when R₀<1, and a unique endemic equilibrium point results whenever R₀>1. The unique endemic equilibrium point is locally asymptotically stable when R₀ is near one. Observably, as R₀ decreases below one (R₀<1), no endemicity exists, and the disease wanes. As R₀ rises above one (R₀>1), the disease spreads through the population.

Figure 2: Forward bifurcation diagram for the model system (1) illustrating how that the disease-free equilibrium and endemic equilibria exchange stability when R₀=1. The cyan and blue curve indicate stable equilibria whereas the dashed red curve signifies unstable equilibrium.

Global Stability of Endemic Equilibrium Point

This subsection deploys the Lyapunov function method to study the global stability of the Endemic Equilibrium Point (EEP) beyond just the issue of neighborhood equilibrium points. The Lyapunov 25 function enables the extension of the analysis to capture an enormous region, as opposed to relying only on a strip region.

Theorem 3.3. If R_e>1, then the bovine tuberculosis disease model system (1) has a unique equilibrium point E^∗ which is globally asymptotically stable.

Proof. Since a suitable Lyapunov function is developed using approach, the Lyapunov function will then be applied to analyze the global stability of the endemic equilibrium point of the model system (1). Therefore, the Lyapunov function is developed using the general form as illustrated below:

where c_i is a positive constant, xi is the population of the ith compartment and x^∗_i is the endemic equilibrium point. Thus, the following Lyapunov function is constructed:

Since the chosen function L and its constant are always differentiable and continuous, the differentiating L regarding time yields the following results:

Substituting the model equations of the model system (1) into the equation (18) to obtain

Now, at the endemic equilibrium point of the model system (1) the following are results.

Case 1: Considering the differential equations of human population, then making π_h, ψ, τ_h, γ_h and ρ_h the subject yields the following:

Case 2: Considering the differential equations of cattle population and making π_c, ω, (μ_c+γ_c) and γ_c the subject yields the following:

Now, substituting equations (20) and (21) into the equation, then dL/dt becomes.

To investigate the global stability of the endemic equilibrium using the Lyapunov function approach at the steady state, set dN∗_h/dt=0 and dN∗_c/dt=0, such that:

Substituting equation (23) and (24) into forces of infections to obtain:

After substituting equation (25) into (26), the equation becomes

Simplifying equation (26), it gives

Therefore, following the idea of the function Φ (S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) is non positive, that is Φ (S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) ≤ 0 for all (S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c). Then, this implies that d_L/d_t ≤ 0 for all (S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) and it is zero only when S_nh=S∗_nh, S_eh=S∗_eh, E_h=E∗_h, I_h=I∗_h, R_h=R∗_h, S_nv=S∗_nv, S_v=S∗_v, E_c=E∗_c, I_c=I∗_c. Therefore, the largest compact invariant set in S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c such that d_L/d_t=0 is the singleton (E∗) which is the endemic equilibrium point of the model system (1). [20] principle, indicates that (E∗) is globally asymptotically in the interior region of S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c. Now recalling the definition of Lyapunov stability, if L(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c)>0 ∀(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) ∈ Ω \{E∗}, L(E∗)=0 and L′(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) ≤ 0 ∀(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) ∈ Ω \{E∗} then, E∗=(S∗_nh, S∗_eh, E∗_h, I∗ _h, R∗ _h, S∗_nv, S∗_v , E∗_c , I∗_c) is stable. Furthermore, if L(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c)>0 ∀(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) ∈ Ω\ {E_∗}, L(E_∗)=0 and L′ (S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c)< 0 ∀(S_nh, S_eh, E_h, I_h, R_h, S_nv, S_v, E_c, I_c) ∈ Ω\{E_∗} then, E_∗=(S∗_nh, S∗_eh, E∗_h, I∗_h, R∗_h, S∗_nv, S∗_v, E∗_c, I∗_c) is asymptotically stable.

Sensitivity and Numerical Analysis

Sensitivity Analysis: The sensitivity analysis helps identify the influence of model parameters on the effective reproduction number (R_e) and disease transmission. It determines which parameters and initial conditions affect the model output. Additionally, it informs researchers about which parameters require more numerical attention. The normalized forward sensitivity index of the variable Re depends on the differentiability of a parameter p and is defined as follows:

where ψ is a parameter present in effective reproduction number R_e. For example, the sensitivity index of R_e corresponding to the parameter α is given as

Other indices are calculated using a similar approach and the results are displayed in Table 2 and Figure 3.

Table 2: Sensitivity indices of R_e using parameter values in Table 1.

Figure 3: Graph of sensitivity indices of R_e with respect to the model parameters.

Interpretation of the Sensitivity

Figure 3 displays the sensitivity profile of R_e concerning the model parameters found within R_e. Further analysis reveals that the parameters a₂, a₁, a, ε₁, ε₂, ε, γ_h, γ_c, and μ_h have positive indices, while μ_c, ρ_h, α_h, alpha_c, ψ, λ, ω, and θ have negative indices. Notably, parameters γ_c and ε₂ exhibit index values of +0.999997 and +0.999994, respectively, indicating that an increase in these parameters, while keeping other variables constant, elevates the effective reproduction number R_e.

In other words, increasing these parameters heightens the risk of a bTB outbreak in the wider population. To minimize infections in the cattle population, the rates γ_c and ε₂ must be kept sufficiently low to ensure that the removal rate from the exposed to infectious class remains minimal, while also maintaining a low contact rate in contaminated environments and inter-cattle transmission.

On the contrary, the vaccine parameter ω and the induced death rate αc emerge as the most negatively sensitive parameters, with index values of α_c=−0.999997 and ω= −0.930551, respectively. This implies that increasing the value of ω by vaccinating healthy cattle and the value of α_c by culling infected animals, while keeping other parameters constant, reduces the effective reproduction number R_e, thus alleviating the disease burden among the cattle population and promoting disease free conditions in the human population.

Numerical Simulation

The most effective approach for accurately solving an Ordinary Differential Equation (ODE) is to carefully develop an exact solution, as discussed in. The chosen method must always uphold the standard properties of the approximated solution, including consistency and convergence, and it must also preserve the qualitative properties of the solution, such as boundedness and positivity, as highlighted. In this section, we employed the Runge-Kutta fourth (RK4) order method to obtain mathematical results for the model system due to its heuristic properties within the ODE framework. This method is known for its stability and practicality, even with large time steps, as demonstrated in previous studies. We utilized the ODE45 version in MATLAB software to execute the computations using the parameter values presented in Table 1. The initial conditions for the state variables were set as follows: S_nh=15000, S_eh=1000, E_h=1000, I_h=500, R_h=300, S_nv=15000, S_v=10, E_c=6000, and I_c=10. These initial conditions were arbitrarily chosen to illustrate specific behaviors of the model system (Figure 4).

Evolution of Population against Time

Figure 4: Dynamics of human and cattle populations over time with interventions.

Figure 4 depicts the behavior of the infected human and cattle populations as time progresses. Figures 4(a) and 4(b) indicate that interventions for both the human and cattle populations reduce the number of infected individuals significantly. This result attests to the existence of the DFE point in the model system (1). When ψ=0, ρ_h=0, and ω=0, it implies that interventions are not in place, meaning the populations of humans and cattle can never reach zero, no matter how long the EE point of the model system persists.

Simulation on the effects of varying some parameter values on the model this section presents data derived from the performance of the simulation of model system to explore the behavior of certain state variables over time. This simulation involved varying selected parameter values, as presented in Table 1. The simulation results are graphically displayed in Figure 5.

Figure 5: (a) Effect of education efficacy (θ) on infected human population and (b) Effect of vaccination efficacy (λ) on infected cattle population.

Figure 5 shows that education efficacy (θ) and vaccine efficacy (λ) play pivotal roles in reducing infectious hosts in humans and cattle. Thus, as education and vaccine efficacy approach 1, the number of infectious humans and cattle decreases.

From Figure 6, it can be surmised that the application of rifampicin, isoniazid, ethambutol, and pyrazinamide as medication for infected humans at the rate ρ_h and Bacille Calmette Gu´erin (BCG) vaccination for infected cattle at the rate ω leads to a reduction in the number of infected humans and cattle when employed on a large scale and applied effectively.

Figure 6: (a) Effect of treatment (ρh) on recovered human population and (b) Effect vaccine (ω) on cattle population.

Figure 6(b) reveals that when cattle are vaccinated to at least ω ≥ 0.4 with a minimum dose of 20 mg, the number of healthy cattle increases rapidly. Additionally, the number of recovered humans also rises, as shown in Figure 6(a). Similarly, when treatment implementation reaches at least ρ_h ≥ 0.4, it triggers a significant increase in the number of recovered humans.

Simulation on Environmental Contamination

This section presents results from the simulation of the infected human and cattle populations when pastures, soil, slurry, and hay are presumably heavily contaminated with M. bovis, with potential infections stemming from sputum, pus, urine, feces, and other excretions of infectious animals (Figure 7).

Figure 7: Impacts of contaminated environment (ϵ₂) on (a) Infected human population and (b) Infected cattle population.

Figure 7(b) illustrates that the number of infected cattle rises as the grazing rate (ϵ₂) on the heavily contaminated environment increases. However, after several months, the number of infected cattle decreases due to the disease.

Simulation on Multiple Intervention Strategies

This section presents data derived from simulating the infected human and cattle populations. It explores the outcomes resulting from the execution of various intervention strategies simultaneous (Figure 8).

Figure 8: (a) Impacts of multiple interventions on infected human population and (b) Single intervention on infected cattle population.

Figure 8 shows the effects of a different combination of interventions (public health education campaigns, treatment, and vaccination) on the bTB transmission dynamics. Figure 8 shows that combining all three intervention strategies decreases disease transmission in the endemic area faster than using only two interventions.

Conclusion

This study has presented the formulation of the model for bTB disease transmission in the presence of control strategies, which include public health education campaigns, treatment, and vaccination. Specifically, it has clearly demonstrated the positivity and boundedness of the model system (1) for the model’s domain to be biologically and mathematically meaningful. The results stemming from the sensitivity analysis on all model parameters, aimed at determining their relationship with the effective reproduction number (R_e), show that the probability of cattle contracting M. bovis from the contaminated environment, inter-cattle transmission (ε₂), and the removal rate from the exposed class to the infectious class (γ_c) are the most positively sensitive parameters for the effective reproduction number. In contrast, the most negatively sensitive parameters are ω (vaccination) and α_c (induced disease death rate), indicating that all infectious cattle are subject to culling in the absence of curative medication administration. Moreover, the numerical simulation of the model reveals that the combination of all interventions has the most significant impact on the control of bTB disease. Therefore, these control measures should be implemented concurrently, especially in endemic areas, to effectively control bTB disease transmission.

Funding Statement

There is no any fund which has been given to support this study.

Declaration of Competing Interest

Authors declare that they have no conflict of interest.

Credit Authorship Contribution Statement

Sylas Oswald: Conceptualization, methodology, software, formal analysis, visualization, investigation and writing original draft

Theresia Crispin Marijan and Goodluck Mika Mlay: Conceptualization, methodology, formal analysis, investigation, writing, review, editing, resources and supervision.

Winifrida Benedict Kidima: Biology aspect visualization, investigation, writing, review, editing and supervision.

Acknowledgment

We acknowledge and thank the College of Natural and Applied Science (CoNAS) for the regular moral supports given to the all authors during writing this paper.

References

1. Aldila D, Latifah SL, Dumbela PA (2019) Dynamical analysis of mathematical model for Bovine Tuberculosis among human and cattle population. Commun Biomath Sci. 2(1): 55-64.

   [Crossref] [Google Scholar]

2. Al-Obeidi AS, Al-Azzawi SF (2022) A novel six-dimensional hyperchaotic system with selfexcited attractors and its chaos synchronisation. Int J Comput Sci Math. 15(1):72-84.

   [Crossref] [Google Scholar]

3. Ashgi R, Pratama MAA, Purwani S (2021) Comparison of Numerical Simulation of Epidemiological Model between Euler Method with 4th Order Runge Kutta Method. Int J Oper Res. 2(1):37-44.

   [Crossref] [Google Scholar]

4. Ayele WY, Neill SD, Zinsstag J, Weiss MG, Pavlik I (2004) Bovine tuberculosis: an old disease but a new threat to Africa. Int J Tuberc Lung Dis. 8(8):924-937.

   [Google Scholar] [PubMed]

5. Butcher JC (1996) A history of Runge-Kutta methods. Appl Numer Math. 20(3):247-260. [Crossref]

   [Google Scholar]

6. Castillo-Chavez C, Song B (2004) Dynamical models of tuberculosis and their applications. Math Biosci Eng. 1(2):361.

   [Google Scholar] [PubMed]

7. Mdegela RH, Kusiluka LJM, Kapaga AM, Karimuribo ED, Turuka FM, et al. (2004) Prevalence and determinants of mastitis and milk borne zoonoses in smallholder dairy farming sector in Kibaha and Morogoro districts in Eastern Tanzania. J Vet Med B. 51(3):123-128.
8. Chitnis N, Hyman JM, Cushing JM (2008) Determining Important Parameters in the Spread of Malaria through the Sensitivity Analysis of a Mathematical Model. Bull Math Biol. 70(5):1272.

   [Crossref] [Google Scholar] [PubMed] 

9. Diekmann O, Heesterbeek JAP, Metz JA (1990) On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol. 28(4):365-382.

   [Crossref] [Google Scholar] [PubMed]

10. Dumont Y, Chiroleu F, Domerg C (2008) On a temporal model for the Chikungunya disease: modelling, theory and numerics. Math Biosci. 213(1):80-91.

    [Crossref] [Google Scholar] [PubMed]

11. Durnez L, Sadiki H, Katakweba A, Machang’u RR, Kazwala RR, et al. (2009) The prevalence of Mycobacterium bovis-infection and a typical mycobacterioses in cattle in and around Morogoro, Tanzania. Trop Anim Health Prod. 41(8):1653-1659.

    [Crossref] [Google Scholar] [PubMed]

12. Ekwem D, Morrison TA, Reeve R, Enright J, Buza J, et al. (2021) Livestock Movement Informs the Risk of Disease Spread in Traditional Production Systems in East Africa. Sci Rep.

    [Crossref] [Google Scholar] [PubMed]

13. Hassan AS, Garba SM, Gumel AB, Lubuma JS (2014) Dynamics of Mycobacterium and bovine tuberculosis in a human-buffalo population. Comput Math Method Med.

    [Crossref] [Google Scholar]

14. Hyera EM, Athumani SN (2017) The status and characteristics of two populations of small East African goats of Tanzania.

    [Google Scholar]

15. Kamgang JC, Sallet G (2005) Global asymptotic stability for the disease free equilibrium for epidemiological models. C R Math. 341(7):433-438.

    [Crossref] [Google Scholar]

16. Karimuribo ED, Ngowi HA, Swai ES, Kambarage DM (2007) Prevalence of brucellosis in crossbred and indigenous cattle in Tanzania. 19(10):148-152.

    [Google Scholar]

17. Katale BZ, Mbugi EV, Matee MI, Fyumagwa RD, Keyyu JD, et al. (2012) Bovine tuberculosis at the human-livestock-wildlife interface: is it a public health problem in Tanzania? A review: proceedings. 79(2):1-8.

    [Google Scholar]

18. Korobeinikov A, Wake GC (2002) Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models. Appl Math Lett. 15(8):955-960.

    [Crossref] [Google Scholar]

19. LaSalle JP (1976) The stability of dynamical systems.

    [Google Scholar]

20. Liu S, Li A, Feng X, Zhang X , Wang K (2016) A dynamic model of human and livestock tuberculosis spread and control in Urumqi, Xinjiang, China. Comput Math Method Med.

    [Crossref] [Google Scholar] [PubMed]