Dynamics of a Mathematical Model for Tuberculosis with Variability in Susceptibility and Disease Progressions Due to Difference in Awareness Level

This work extends a mathematical model for the transmission dynamics of tuberculosis that examined the impact of certain factors on tuberculosis case detection (Okuonghae and Omosigho, 2011). The extended model now classifies the latently infected individuals by their level of tuberculosis awareness (as was done for the susceptible sub-population) and further expands the number of key factors that can positively affect the tuberculosis case detection rate. The effect of these identified factors on the associated reproduction number of the model is considered. It is shown that the system can undergo the phenomenon of backward bifurcation when the associated reproduction number of the model is less than unity; in a special case, the effect of exogenous re-infection on the backward bifurcation phenomenon is significantly dictated by the level of awareness of the latently infected individuals. Qualitative and quantitative analysis of the model showed the effect of key identified factors on the dynamics of tuberculosis while suggesting a serious concentration on tuberculosis awareness programmes, active case finding strategies and use of active cough identification for identifying likely TB cases and sustaining awareness campaigns over a long period of time.

1. Introduction

Tuberculosis (TB), an infectious disease caused by the Mycobacterium tuberculosis bacillus, remains one of the world's deadliest diseases (World Health Organization, 2014). According to the World Health Organization (WHO), in 2013, about 9 million people were infected, worldwide, with TB and 1.5 million deaths from the disease were reported, 360,000 of whom were HIV-positive (World Health Organization, 2014). Tuberculosis is seen to be declining slowly each year and an estimated 37 million lives were saved between 2000 and 2013 through effective diagnosis and treatment (World Health Organization, 2014). On the average, TB incidence fell to about 1.5% per year, between 2000 and 2013, worldwide (World Health Organization, 2014). Globally, TB mortality rate fell by an estimated 45% between 1990 and 2013 while the prevalence rate dropped by 41% (World Health Organization, 2014). Even with such positive results achieved within the last 14 years, it is still thought that deaths from tuberculosis are preventable; in fact the death toll is still considered unacceptably high. Hence, efforts are geared toward accelerating programmes that will result in a reduction in the TB burden globally [within the context of the Millennium Development Goals (MDGs)], and reach the Stop TB Partnership target of a 50% reduction by 2015 (World Health Organization, 2014).

More than half of the approximately 9 million individuals who are infected with tuberculosis in 2013 (56%) were in South-East Asia and Western Pacific Regions. A further one quarter of these infected individuals are in the African Region, which account for the highest rates of TB cases and deaths relative to population (World Health Organization, 2014).

Generally, reducing the incidence and prevalence of TB in a population hinges on successful treatment and high case detection rates (Okuonghae and Omosigho, 2010, 2011; World Health Organization, 2014; Okuonghae, 2015). It fact, it is seriously encouraged that effort should be concentrated on ensuring that all TB cases are detected, notified and commence treatment immediately (World Health Organization, 2014). It is reported that about 6.1 million TB cases were reported to the WHO in 2013 out of which about 5.7 million were individuals were newly diagnosed and another 0.4 million were already on treatment (World Health Organization, 2014). Tuberculosis notification has stabilized in recent years, with about 64% of the estimated 9 million individuals who developed TB in 2013 were notified as newly diagnosed cases (World Health Organization, 2014). This implies that about 3 million cases were either not diagnosed, or diagnosed but not reported to national TB programmes (World Health Organization, 2014) which could hinder the goal of significantly reducing the prevalence of TB in such localities. Treatment success rates (globally) have been impressive (and continue to be high) over the years, with about 86% treatment success rate reported in 2013 among all new TB cases (World Health Organization, 2014).

Tuberculosis affects family and social relationships and results in adverse health and economic consequences (Armijos et al., 2008; Chang and Cataldo, 2014). As stated earlier, improving the case detection and notification rates will result in reducing the TB burden in a population. However, several factors could be hindering efforts at improving these rates. Individuals infected with TB and their families can experience prejudice and negative attitudes, such as shame, blame and a sense of judgment as a result of the infection (Bennstam et al., 2004; Baral et al., 2007; Chang and Cataldo, 2014). Stigmatization can also be a stumbling block in improving the tuberculosis case detection rate (Johansson et al., 2000; Okuonghae and Omosigho, 2010; Chang and Cataldo, 2014) since patients and families' fears of inferiority stems from the anticipation of an adverse judgement related to a TB diagnosis (Johansson et al., 2000; Chang and Cataldo, 2014).

A survey reported in Okuonghae and Omosigho (2010) listed some factors that can adversely affect the implementation of the directly observed treatment, short-course (DOTS) strategy in Nigeria (one of the high burden countries) in reducing the incidence of TB in the country. The survey revealed that most persons do not know how TB is transmitted and the signs and symptoms of tuberculosis; several individuals are not even aware of the government's health policies on tuberculosis and TB treatment. Further, the survey revealed that this lack of awareness can lead to delays in reporting likely TB cases for treatment (Okuonghae and Omosigho, 2010), increasing the likelihood of disease transmission.

Analysis of the results from the survey (Okuonghae and Omosigho, 2010) identified four key factors that must be combined for an effective control of tuberculosis: “effective awareness programme, active cough identification, associated cost factor for treatment of identified cases and effective treatment” (Okuonghae and Omosigho, 2010).

A mathematical model for TB dynamics that incorporated the identified factors (as parameters) gleaned from the work in Okuonghae and Omosigho (2010) was formulated and analyzed in Okuonghae and Omosigho (2011). Control strategies, based on the identified parameters, that can lead to minimizing the incidence (as well as the prevalence) of the disease in a population were proposed. In summary, the qualitative and quantitative analysis of the model in Okuonghae and Omosigho (2011) showed that a serious concentration on tuberculosis awareness programmes and active cough identification as a marker for identifying potential TB cases can be significant in minimizing the severity of the disease in a population, with effective treatment.

The purpose of this article is to present and analyze a new mathematical model for TB dynamics; this model is an extension of that presented and analyzed in Okuonghae and Omosigho (2011). The aim of this work is to further study the effects of additional heterogeneities based on the level of awareness of tuberculosis within the population and active-case finding, on the dynamics of the disease. In the work in Okuonghae and Omosigho (2011), only the susceptible subpopulation was stratified by their level of tuberculosis awareness. In this work, we will now stratify both the susceptible and latently infected sub-populations by their level of awareness of TB (symptoms and signs of TB as well as testing and treatment programmes available by the government). Note that the measure of the case detection and notification rates will be by the number of infectious individuals detected, notified, and treated for tuberculosis.

2. Materials and Methods

2.1. Basic Model

This section briefly describes the model in Okuonghae and Omosigho (2011). In Okuonghae and Omosigho (2011), we assumed that susceptible individuals are divided into two groups depending on their level of awareness of the disease (and any treatment policy): the high risk (low level of awareness) group, S₁, and the “educated,” low risk (high level of awareness) group, S₂. The S₁ class is “educated” at the per capita rate α₁ and thereafter move into the S₂ class. Tuberculosis infection can invade the S₁ and S₂ classes, depending on the “efficacy” of the education programme. The programme is assumed to reduce the likelihood of infection by a factor of σ (0 ≤ σ ≤ 1). The case σ = 0 signifies a completely effective education program, while σ = 1 models the situation where the program is totally ineffective.

It was further assumed that the “vaccine” (education program) produces temporary immunity at the per capita rate θ. The case θ = ∞ corresponds to the case where there is absolutely no immunity while θ = 0 corresponds to life-long immunity. Hence, θ measures the rate at which those in the S₂ class return to the S₁ class due to forgetfulness caused by lack of continuous exposure to the enlightenment program while the disease persists in the community. We assumed β to be the disease transmission rate.

We also assumed that the variables E, I, J, and T represented the “primary” latent, infectious, identified infectious (for treatment under DOTS) and the effectively treated individuals, respectively (by “primary” latency, we are referring to susceptible individuals who are infected for the first time as well as treated individuals who now get reinfected after recovering from a previous infection). In addition to these groups, a separate class (R(t)) was added to account for individuals who become latent due to failed treatment or self cure (we refer to this situation as “secondary” latency).

The parameter η was used as a modification parameter (0 ≤ η ≤ 1) to account for the relative infectiousness of infectious individuals in the J class. We also assumed that ϵ is the reduced likelihood of reinfection of effectively treated individuals, where 0 ≤ ϵ ≤ 1. Also, we took 0 < p < 1 as the fraction of individuals with new infections who develop TB fast per unit of time, with Λ being the rate of recruitment of uninfected newborns and immigrants into the low risk susceptible class. For simplifications, we assumed that all entrants into the population move into the S₁ group. We also assumed that μ is the natural death rate while d is the tuberculosis-induced death rate.

To account for exogenous reinfection of latently infected individuals, we assumed that β* be the transmission rate amongst this group in infected persons. We also assumed that k is the rate of progression of infected individuals in the latent stage to active tuberculosis.

Individuals with active TB can be identified using chronic cough lasting more than 2 weeks as a marker, at the rate α₂, and are referred to a TB treatment program under DOTS for effective treatment (in Okuonghae and Omosigho, 2011, α₂ was known as the cough identification rate). However, a fraction of these identified cases will eventually get into the treatment program when we consider the cost factor. Hence, a cost improvement factor (ν : 0 < ν ≤ 1) will affect the actual number of identified cases that commences treatment. The cost factor considers the effect the actual cost of medical tests and treatment will have on the care givers when presenting the infectious individual for treatment. If ν = 0, then the cost of medical tests and treatment is prohibitively high and να₂ = 0 implies that the TB case will not get into a TB treatment program due to the financial cost on the care givers or family members. However, if ν = 1, it means that the cost of medical tests and treatment is totally free and να₂I will be the total number of identified cases that are tested and treated for tuberculosis under DOTS. Hence, να₂ was taken to be the proportion of identified TB cases that commences treatment under DOTS. Since in most developing countries, like Nigeria, it was observed that medical tests for TB still cost money especially when the infectious person is about commencing treatment (Okuonghae and Omosigho, 2010), it then implies that ν ≠ 1; it lies in the range 0 < ν < 1. In all, ν → 1 implies that the cost of testing and treating TB becomes affordable.

We assumed that r₂ is the treatment rate for the identified infectious individuals under the DOTS scheme while the fraction of the detected cases who were successfully treated under the DOTS was n, with m = 1 − n being the fraction of those whose treatment were unsuccessful and, thereafter, moved to the “secondary” latency group.Tuberculosis cases that are not detected either die at the rate d, or self-cure and revert to the “secondary” latent state (in R) at the rate r₁.

The mathematical model was then given by the following system of non-linear ordinary differential equations (Okuonghae and Omosigho, 2011):

$\begin{array}{cc}\frac{d{S}_1}{dt}=\Lambda -{\alpha }_1{S}_1-\beta S_1\frac{(I+\eta J)}{N}+\theta S_2-\mu S_1,\hfill & (1\text{a})\hfill \end{array}$

$\begin{array}{cc}\frac{d{S}_2}{dt}={\alpha }_1{S}_1-\sigma \beta S_2\frac{(I+\eta J)}{N}-\theta S_2-\mu S_2,\hfill & \left(1\text{b}\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{dE}{dt}=(1-p)(\beta S_1\frac{(I+\eta J)}{N}+\sigma \beta S_2\frac{(I+\eta J)}{N}+ϵ\beta T\frac{(I+\eta J)}{N})\hfill \\ -{\beta }^{\ast }E\frac{(I+\eta J)}{N}-(k+\mu )E,\hfill \end{array}& \left(1\text{c}\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{dI}{dt}=p(\beta S_1\frac{(I+\eta J)}{N}+\sigma \beta S_2\frac{(I+\eta J)}{N}+ϵ\beta T\frac{(I+\eta J)}{N})+kE\hfill \\ +{\beta }^{\ast }E\frac{(I+\eta J)}{N}-(\nu {\alpha }_2+\mu +d+r_1)I,\hfill \end{array}& \left(1\text{d}\right)\hfill \end{array}$

$\begin{array}{cc}\frac{dJ}{dt}=\nu {\alpha }_{2}I-r_{2}J-\mu J,\hfill & \left(1\text{e}\right)\hfill \end{array}$

$\begin{array}{cc}\frac{dT}{dt}=n{r}_{2}J-\mu T-ϵ\beta T\frac{(I+\eta J)}{N},\hfill & \left(1\text{f}\right)\hfill \end{array}$

$\begin{array}{cc}\frac{dR}{dt}=r_{1}I+m{r}_{2}J-\mu R.\hfill & \left(1\text{g}\right)\hfill \end{array}$

The effective reproduction number of model (1) is given as

$\begin{array}{cc}\begin{array}{l}{{R}}_c=\frac{\beta }{d+\mu +r_1+\nu {\alpha }_2}\frac{(k+p\mu )}{k+\mu }\frac{(\mu +r_2+\eta \nu {\alpha }_2)}{\mu +r_2}\\ \left[\frac{\mu +\theta }{\mu +{\alpha }_1+\theta }+\frac{\sigma {\alpha }_1}{\mu +{\alpha }_1+\theta }\right]\end{array}& \left(2\right)\end{array}$

See Okuonghae and Omosigho (2011) and Okuonghae (2015) for the qualitative and quantitative analysis of model (1) and the effect of the key parameters, gleaned from the survey in Okuonghae and Omosigho (2010), on the dynamics of TB in a population.

2.2. Modified Mathematical Model

Instead of stratifying only the susceptible subpopulation by their level of awareness (i.e., S₁ and S₂ as described in Section 2.1), we will also stratify the latently infected class by their level of awareness. This is reasonable since latently infected individuals do not transmits TB (they do not show the signs and symptoms of TB) and, in most cases, will become aware of their disease status only when tuberculin tests are carried out on them; for example, the Mantoux tuberculin skin test (TST). Figure 1 shows the schematic diagram of the modified mathematical model given in (3).

FIGURE 1

Figure 1. Schematic diagram of model given in (3), where $\lambda =\frac{\beta (I+\eta J)}{N}$.

In addition to some of the epidemiological classes in Section 2.1, we will break up the latent class (E) into two classes, depending on their level of awareness of the disease: the latently infected, high risk (low level of awareness) group, E₁, and the “educated” latently infected, low risk (high level of awareness) group, E₂. It can be stated that due to awareness, individuals in the E₂ group are tested in order to know their disease status and take necessary precautionary health measures. The E₁ class is made up of individuals from the S₁ class and some individuals from the treated class (T) who have a low level of awareness of tuberculosis while the E₂ class is made up of individuals from the S₂ class and some persons from the treated class (T) who retained their high level of awareness following their recovery. The latently infected individuals with low awareness (E₁) can become “educated” at the per capita rate ψ and thereafter move into the E₂. Individuals who recover (following effective treatment) can get reinfected, with a fraction 0 ≤ ω ≤ 1 of such persons entering the class of latent individuals with low level of awareness (E₁) while the remaining fraction 1 − ω enter the E₂ class.

We also assume that the education programmes produces “temporary immunity” at the per capita rate θ₁ (for the susceptible groups, S₁ and S₂) and at the rate θ₂ (for the latent groups, E₁ and E₂). The cases θ₁(θ₂) = ∞ corresponds to the situation where there is absolutely no immunity (from the education programmes) while θ₁(θ₂) = 0 corresponds to life-long immunity. Hence, θ₁(θ₂) measures the rate at which those in the S₂(E₂) class return to the S₁(E₁) class due to a lack of continuous exposure to the awareness programmes while the disease persists in the community. We further assume that β is the disease transmission rate.

Also, we assume that 0 < p₁(p₂) < 1 represent the fraction of persons with new infections who develop TB fast per unit of time from the class of infected individuals with low level of awareness (high level of awareness). We expect that, by virtue of the benefits of awareness, p₁ ≥ p₂ as new infections are quickly detected and fewer cases of fast progressions to active TB are recorded amongst individuals with a high level of awareness; also, note that the period of fast latency can span between 1 and 5 years (Styblo, 1980; Flynn and Chan, 2001; Colijn et al., 2007) so that early detection can prevent more cases of active TB.

We further assume that b₁(b₂) are modification parameters accounting for exogenous re-infection of the latently infected individuals in the E₁(E₂) class, with 0 ≤ b₂ ≤ b₁ < 1 and that k₁(k₂) is the rate of progression of the individuals in the latent state (E₁(E₂)) to active tuberculosis. Assume that, in addition, to the impact of active cough identification (α₂) and the cost factor (ν) on improving the case detection (and notification) rates, evident by the number of infectious individuals in the J class (and treatment), let π be the rate at which an active case-finding strategy is used in searching for infectious TB cases for onward treatment with the treatment rate now represented as r. The remaining parameters in the modified model are as defined in Section 2.1. Since there is significant improvement in treatment success rate for TB (worldwide) (World Health Organization, 2014), especially in countries and communities, where there is an efficient treatment strategy put in place (for example, DOTS), we assume an insignificant number of failed treatments and self cure and will omit the R class from the modified model.

On the basis of the above assumptions, the modified model is now given by the following system of non-linear ordinary differential equations:

$\begin{array}{cc}\frac{d{S}_1}{dt}=\Lambda -{\alpha }_1{S}_1-\beta S_1\frac{(I+\eta J)}{N}+{\theta }_1{S}_2-\mu S_1,\hfill & \left(3\text{a}\right)\hfill \end{array}$

$\begin{array}{cc}\frac{d{S}_2}{dt}={\alpha }_1{S}_1-\sigma \beta S_2\frac{(I+\eta J)}{N}-{\theta }_1{S}_2-\mu S_2,\hfill & \left(3\text{b}\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{d{E}_1}{dt}=(1-p_1)\left(\beta S_1\frac{(I+\eta J)}{N}+\omega ϵ\beta T\frac{(I+\eta J)}{N}\right)\hfill \\ -b_1\beta E_1\frac{(I+\eta J)}{N}-(k_1+\mu +\psi )E_1+{\theta }_2{E}_2,\hfill \end{array}\hfill & \left(3\text{c}\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{d{E}_2}{dt}=(1-p_2)\left(\sigma \beta S_2\frac{(I+\eta J)}{N}+(1-\omega )ϵ\beta T\frac{(I+\eta J)}{N}\right)\hfill \\ -b_2\beta E_2\frac{(I+\eta J)}{N}-(k_2+\mu +{\theta }_2)E_2+\psi E_1,\hfill \end{array}\hfill & \left(3\text{d}\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{dI}{dt}=p_1\beta (S_1+\omega ϵT)\frac{(I+\eta J)}{N}+p_2\beta (\sigma S_2\hfill \\ +(1-\omega )ϵT)\frac{(I+\eta J)}{N}+\beta (b_1{E}_1+b_2{E}_2)\frac{(I+\eta J)}{N}\hfill \\ +k_1{E}_1+k_2{E}_2-(\nu {\alpha }_2+\mu +d+\pi )I,\hfill \end{array}& \left(3\text{e}\right)\hfill \end{array}$

$\begin{array}{cc}\frac{dJ}{dt}=(\pi +\nu {\alpha }_2)I-rJ-\mu J,\hfill & \left(3\text{f}\right)\hfill \end{array}$

$\begin{array}{cc}\frac{dT}{dt}=rJ-\mu T-ϵ\beta T\frac{(I+\eta J)}{N},\hfill & \left(3\text{g}\right)\hfill \end{array}$

with N = S₁ + S₂ + E₁ + E₂ + I + J + T.

2.2.1. Basic Properties of Model (3)

For model (3) to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all time t. In other words, the solutions of model (3) with positive initial data will remain positive for all time t ≥ 0.

Theorem 2.1. Let the initial data for model (3) be S₁(0) > 0, S₂(0) > 0, E₁(0) > 0, E₂(0) > 0, I(0) > 0, J(0) > 0 and T(0) > 0. Then, the solutions

$(S_1(t),S_2(t),E_1(t),E_2(t),I(t),J(t),T(t))$

of model (3), with positive initial data, will remain positive for all time t > 0.

Theorem 2.2. The closed set

${D}=\left\{(S_1,S_2,E_1,E_2,I,J,T)\in {ℝ}_{+}^7;N\le \frac{\Lambda }{\mu }\right\}$

is positively invariant and attracts all positive solutions of model (3).

See Appendix A for the proofs of Theorems 2.1 and 2.2.

Since the region ${D}$ is positively invariant, the unique solution of model (3) exists and depends continuously on the initial data of the model (hence, it is sufficient to study its asymptotic dynamics in the region ${D}$, Hethcote, 2000).

2.2.2. Local Asymptotic Stability of Disease-free Equilibrium

Model (3) has a disease-free equilibrium (DFE), obtained by setting the right hand sides of the equations in model (3) to zero and solve for the state variables with no infections, given by

$\begin{array}{l}{\xi }_1=(S_1^{\ast },S_2^{\ast },E_1^{\ast },E_2^{\ast },I^{\ast },J^{\ast },T^{\ast })\\ =\left(\frac{\Lambda }{\mu }\frac{\mu +{\theta }_1}{\mu +{\alpha }_1+{\theta }_1},\frac{\Lambda }{\mu }\frac{{\alpha }_1}{\mu +{\alpha }_1+{\theta }_1},0,0,0,0,0\right)\end{array}$

We see that the susceptible classes, at the DFE, depend on a factor of the asymptotic population size, $\frac{\Lambda }{\mu }$.

The local stability of ξ₁ can be established with the next generation operator method on system (3) (Diekmann et al., 1990; van den Driessche and Watmough, 2002). Using the notations in van den Driessche and Watmough (2002), it follows that matrices F and V, for the new infection terms and the remaining transition terms, respectively, are given by

$F=\left(\begin{array}{ccccc}0& 0& (1-p_1)\beta \frac{S_1^{\ast }}{N^{\ast }}& (1-p_1)\beta \eta \frac{S_1^{\ast }}{N^{\ast }}& \\ 0& 0& (1-p_2)\beta \sigma \frac{S_2^{\ast }}{N^{\ast }}& (1-p_2)\beta \eta \sigma \frac{S_2^{\ast }}{N^{\ast }}& \\ 0& 0& \beta \frac{(p_1{S}_1^{\ast }+p_2\sigma S_2^{\ast })}{N^{\ast }}& \beta \eta \frac{(p_1{S}_1^{\ast }+p_2\sigma S_2^0)}{N^{\ast }}& 0\\ 0& 0& 0& 0& \end{array}\right)$

and

$V=\left(\begin{array}{cccc}(k_1+\mu +\psi )& \text{}\text{}-{\theta }_2& 0& 0\\ \text{}\text{}-\psi & \text{}\text{}(k_2+\mu +{\theta }_2)& 0& 0\\ \text{}\text{}{k}_1& -k_2& \text{}\text{}(\nu {\alpha }_2+d+\mu +\pi )& 0\\ \text{}\text{}0& 0& -(\pi +\nu {\alpha }_2)& \text{}\text{}(r+\mu )\text{}\text{}\end{array}\right)\text{}.$

Thus, the effective reproduction number of model (3), denoted by ${R}$_T, is given by

$\begin{array}{cc}{{R}}_T=\frac{\beta (\mu +\eta \pi +r+\eta \nu {\alpha }_2)[G_1+G_2+G_3]}{G_4},& \left(4\right)\end{array}$

where

$\begin{array}{l}{G}_1=k_1[(\mu +\sigma {\alpha }_1+{\theta }_1)(k_2+{\theta }_2)+\mu (\mu +\sigma {\alpha }_1{p}_2+{\theta }_1)],\\ G_2=\mu (\mu +{\theta }_2+\psi )[\sigma {\alpha }_1{p}_2+p_1(\mu +{\theta }_1)],\\ G_3=k_2(\mu +{\theta }_1)(\mu p_1+\psi )+k_2\sigma {\alpha }_1(\mu +\psi ),\\ G_4=(\mu +r)(d+\mu +\pi +\nu {\alpha }_2)(\mu +{\alpha }_1+{\theta }_1)[(\mu +k_1)\\ (\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ],\end{array}$

and ${R}$_T is obtained from ρ(FV⁻¹) with ρ being the spectral radius of the matrix FV⁻¹.

The following result follows from Theorem 2 in (van den Driessche and Watmough, 2002):

Lemma 2.1. The DFE, ξ₁, of model (3) is locally asymptotically stable if ${R}$_T < 1 and unstable if ${R}$_T > 1.

The threshold quantity, ${R}$_T, represents the average number of secondary tuberculosis infections generated by a typical infectious individual in a completely susceptible population with already existing controls such as treatment (Hethcote, 2000; van den Driessche and Watmough, 2002). The epidemiological implication of Lemma 2.1 is that tuberculosis can be effectively controlled in the community (when ${R}$_T < 1) if the initial sizes of the subpopulation of model (3) are in the basin of attraction of the DFE, ξ₁, in the presence of control strategies including chronic cough identification, active-case finding, effective cost factor, effective awareness programmes and treatment. Hence, a small influx of individuals with active TB into the community will not generate large TB outbreaks, and the disease will die out with time.

2.2.3. Comparing Reproduction Numbers Under Different Control Scenarios

It is important to compare different scenarios involving the presence of control measures based on the reproduction number for the different situations.

Case 1:

Re-write the effective reproduction number as

$\begin{array}{cc}{{R}}_T={{R}}_F\frac{\mu (\pi \eta +\mu +r+\eta \nu {\alpha }_2)}{(\pi \eta +\mu +\eta \nu {\alpha }_2)(\mu +r)}\equiv {{R}}_F{A}_1& \left(5\right)\end{array}$

where

$A_1=\frac{\mu (\pi \eta +\mu +r+\eta \nu {\alpha }_2)}{(\pi \eta +\mu +\eta \nu {\alpha }_2)(\mu +r)}$

and ${R}$_F = ${R}$_T|_{r = 0} is the reproduction number of model (3), under control, without treatment (where ${R}$_T|_{r = 0} is the effective reproduction number, with r = 0).

Since the difference of ${R}$_F from ${R}$_T is only in treatment, the factor A₁ compares a population with and without treatment; however other control strategies such as chronic cough identification (with an effective cost factor) and active-case finding are present in the community. We observe that A₁ will be less than unity (A₁ < 1) for all 0 ≤ η, ν < 1.

Generally, if ${R}$_F < 1, then we cannot expect a TB epidemic and in this case no treatment strategy is needed for control. However, when ${R}$_F > 1, we need to determine the necessary condition for slowing the development of tuberculosis in the community. Following Mukandavire et al. (2007); Sharomi et al. (2008); Okuonghae and Omosigho (2011), we have the difference between ${R}$_F and ${R}$_T as

$\begin{array}{cc}{\Delta }_T:={{R}}_F-{{R}}_T=(1-A_1){{R}}_F.& \left(6\right)\end{array}$

For effective treatment, use of chronic cough as a marker for potential TB cases (with an effective cost factor) and an effective active case-finding strategy to slow down the spread of tuberculosis in a population, we expect that Δ_T > 0, and this condition is satisfied if A₁ < 1 in Equation (6). Now, setting ${R}$_T = 1 and solving for A₁, we have the threshold effectiveness of treatment and cough identification taking into account the cost factor and an active case-finding programme:

$\begin{array}{cc}{A}_1^{\ast }=\frac{1}{{{R}}_F}& (7)\end{array}$

Hence, tuberculosis can be eradicated from the population if the control measures include effective treatment, active chronic cough identification (taking into cognisance an effective cost factor) and an active case-finding strategy if $A_1<A_1^{*}$. Observe from Equation (7) that $A_1^{*}$ is a decreasing function of ${R}$_F; higher values of $A_1^{*}$ results in lower values of ${R}$_F, a desired outcome.

Taking the following limits of A₁ provides further insight into possible ways of reducing the TB burden in a community:

$\begin{array}{cc}(i)\underset{\begin{array}{c}\nu \to 1\\ {\alpha }_2\to \infty \end{array}}{\lim }{A}_1=\underset{\begin{array}{c}\pi \to \infty \\ \nu \to 1\\ {\alpha }_2\to \infty \end{array}}{\lim }{A}_1=\underset{\begin{array}{c}\pi \to \infty \\ \nu \to 1\\ {\alpha }_2\to \infty \\ \eta \to 0\end{array}}{\lim }{A}_1=\frac{\mu }{\mu +r},& \left(8\right)\end{array}$

$\begin{array}{cc}(ii)\underset{r\to \infty }{\lim }{A}_1=\underset{\begin{array}{c}r\to \infty \\ \nu \to 1\end{array}}{\lim }{A}_1=\frac{\mu }{\mu +\pi \eta +\eta {\alpha }_2},(iii)\underset{\begin{array}{c}r\to \infty \\ \eta \to 0\end{array}}{\lim }{A}_1=1.& \left(9\right)\end{array}$

Observe that the limits of A₁ in Equations (8) and (9)(ii) are less than one; hence control strategies for tuberculosis gleaned from these results can be pursued, if it is feasible and practicable economically. In practice, ν → 1 implies near total elimination of costs (medical tests and treatment), α₂ → ∞ implies high rate of identifying “potential” TB cases using chronic cough as a marker, r → ∞ implies high treatment rate, π → ∞ implies a high active case-finding rate for infectious TB individuals and η → 0 implies that those identified for treatment have a significantly reduced likelihood of transmitting the disease.

Using the factor A₁, one observes that an effective combination of chronic cough identification (with the associated cost factor), an active case-finding strategy with effective treatment and taking care to prevent detected infectious TB cases, receiving treatment, from causing new TB infections, will result in reducing the burden of TB in the population.

It is very interesting to observe that more than one combination of control strategies could yield the same positive, desired result. For example, as seen in the limits of A₁ in Equation (8), concentrating on using chronic cough as a marker for a potential TB case with little cost involved in testing and treatment is as effective as including an effective active case-finding strategy to the aforementioned strategy. Also, an effective combination of treatment rate and cost factor can also lead to a reduction in the number of new TB infections in the population.

Case 2:

Consider the situation where there are no control strategies put in place for detecting active TB cases via the use of chronic cough identification and an active case-finding strategy so that π = α₂ = ν = 0. Clearly then, we observe that (J, T) → 0 asymptotically, as t → ∞, in model (3). Therefore, the effective reproduction number of the model (3), with π = α₂ = ν = 0, ${R}$_B, is written as

${{R}}_B={{R}}_T{|}_{\pi ={\alpha }_2=\nu =0}.$

Now, we have that ${R}$_T = A₂${R}$_B where

$A_2=\frac{d+\mu }{d+\pi +\mu +\nu {\alpha }_2}\frac{\pi \eta +\mu +r+\eta \nu {\alpha }_2}{\mu +r}.$

Here, A₂ is known as the effectiveness of treatment, chronic cough identification (with an effective cost factor), an active case-finding strategy, and reduced infectivity of isolated infectious cases undergoing treatment on the control of tuberculosis.

If ${R}$_B < 1, then TB cannot develop into an epidemic and no control strategy involving treatment, chronic cough identification (with a cost factor) and an active case-finding strategy is needed for TB control. If we take the difference between ${R}$_B and ${R}$_T i.e.,

${\Delta }_B:={{R}}_B-{{R}}_T=(1-A_2){{R}}_B$

then an effective combination of the above-mentioned controls for slowing down the spread of TB in the population (as stated above in the expression for A₂) would mean that Δ_B > 0, which is satisfied if A₂ < 1.

It is very interesting to note the expressions for the limits of A₂, when the limits applied to A₁ are used, namely:

$\begin{array}{cc}(i)\underset{\begin{array}{c}\nu \to 1\hfill \\ {\alpha }_2\to \infty \hfill \end{array}}{\lim }{A}_2=\underset{\begin{array}{c}\pi \to \infty \hfill \\ \nu \to 1\hfill \\ {\alpha }_2\to \infty \hfill \end{array}}{\lim }{A}_2=\underset{\begin{array}{c}\nu \to 1\hfill \\ r\to \infty \hfill \end{array}}{\lim }{A}_2=\frac{\eta (d+\mu )}{\mu +r},\hfill & \left(10\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}(ii)\underset{r\to \infty }{\lim }{A}_2=\underset{\begin{array}{c}r\to \infty \hfill \\ \eta \to 0\hfill \end{array}}{\lim }{A}_2=\frac{d+\mu }{d+\mu +\pi +\nu {\alpha }_2}\text{and}\hfill \\ (iii)\underset{\begin{array}{c}\pi \to \infty \hfill \\ \eta \to 0\hfill \\ \nu \to 1\hfill \\ {\alpha }_2\to \infty \hfill \end{array}}{\lim }{A}_2=0.\hfill \end{array}\hfill & \left(11\right)\hfill \end{array}$

Other than the limits of A₂ in Equation (11) (iii), the results of the other limits of A₂ in Equations (10) and (11) (ii) are adjusted by disease-induced death and are less than unity, when compared to the corresponding limits for A₁, as discussed earlier. Interestingly, the limit in Equation (11) (ii) is zero, compared to the same limit that gave $\frac{\mu }{\mu +r}$, for A₁; hence, the former have a lesser value in the limit compared to the latter.

Case 3

In this case, we are considering the situation where there are no awareness programmes for both susceptible and latently infected individuals. However, we assume that treated individuals now have a measure of awareness following their treatment and recovery from TB. Hence, we have that α₁ = ψ = θ₁ = θ₂ = 0. We observe, from model (3) that, for this case, S₂ → 0 asymptotically, as t → ∞, so that the latently infected class E₂ is now populated by treated individuals who now have some level of awareness, probably due to the experience they went through during the course of their treatment. We assume that those who are now populating the E₂ class have life-long awareness with life long benefit (so that θ₂ = 0). Hence, the effective reproduction number of the model (3), with α₁ = ψ = θ₁ = θ₂ = 0, ${R}$_A, is written as

${{R}}_A={{R}}_T{|}_{{\alpha }_1=\psi ={\theta }_1={\theta }_2=0}.$

Now, we have that ${R}$_T = A₃${R}$_A where

$A_3=\frac{B_1}{B_2},$

with B₁ = (μ(μ + k₁)(πη + μ + r + ηνα₂)(k₁(μ(μ + σα₁p₂ + θ₁) + (μ + σα₁ + θ₁)(k₂ + θ₂)) + σk₂α₁(μ + ψ) + k₂(μ + θ₁)(μp₁ + ψ) + μ(σα₁p₂ + p₁(μ + θ₁))(μ + θ₂ + ψ))) and B₂ = [(k₁ + μp₁)(μ + r)(πη + μ + ηνα₂)(μ + α₁ + θ₁)[(μ + k₁)(μ + k₂ + θ₂) + (μ + k₂)ψ]].

In this case, A₃ is known as the effectiveness of treatment, chronic cough identification in the presence of a cost improvement factor, an active case-finding strategy, reduced infectivity of isolated infectious cases (for treatment) and awareness program for both susceptible and latent individuals on the control of tuberculosis.

If ${R}$_A < 1, then TB cannot develop into an epidemic and no control strategy involving the aforementioned factors (as stated in the preceding paragraph) is needed for TB control.

If we take the difference between ${R}$_A and ${R}$_T i.e.,

${\Delta }_A:={{R}}_A-{{R}}_T=(1-A_3){{R}}_B,$

then an effective combination of the controls (as mentioned above and expressed in A₃) for slowing down the spread of TB in the population would mean that Δ_A > 0, which is satisfied if A₃ < 1.

Taking some limits of A₃ threw up some interesting conclusions as to effectively controlling tuberculosis, in this case:

$\begin{array}{cc}(i)\underset{\begin{array}{c}r\to \infty \\ {\alpha }_1\to \infty \\ {\theta }_1\to 0\\ \pi \to \infty \end{array}}{\lim }{A}_3=\underset{\begin{array}{c}\pi \to \infty \\ r\to \infty \\ \psi \to \infty \\ {\theta }_2\to 0\end{array}}{\lim }{A}_3=\underset{\begin{array}{c}\nu \to 1\\ r\to \infty \\ \pi \to \infty \\ {\theta }_2\to 0\\ {\theta }_1\to 0\\ \psi \to \infty \\ {\alpha }_1\to \infty \\ {\alpha }_2\to \infty \end{array}}{\lim }{A}_3=0,& \left(12\right)\end{array}$

$\begin{array}{cc}(ii)\underset{\begin{array}{c}{\alpha }_1\to \infty \\ \psi \to \infty \\ {\theta }_1\to 0\\ {\theta }_2\to 0\end{array}}{\lim }{A}_3=\frac{\mu \sigma (\mu +k_1)(k_2+p_2\mu )(\pi \eta +\mu +r+\eta \nu {\alpha }_2}{(\mu +k_2)(k_1+p_1\mu )(\mu +r)(\pi \eta +\mu +\eta \nu {\alpha }_2)},& \left(13\right)\end{array}$

and

$\begin{array}{cc}(iii)\underset{\begin{array}{c}\nu \to 1\\ \pi \to \infty \\ {\theta }_2\to 0\\ {\theta }_1\to 0\\ \psi \to \infty \\ {\alpha }_1\to \infty \\ {\alpha }_2\to \infty \end{array}}{\lim }{A}_3=\frac{\mu \sigma (\mu +k_1)(k_2+p_2\mu )}{(\mu +k_2)(k_1+p_1\mu )(\mu +r)}.& \left(14\right)\end{array}$

Clearly, the limits of A₃ that includes r (treatment rate) are zero, and this is a positive outcome. However, for the limits of A₃ in Equations (13) and (14), the limiting values of A₃ will be zero if σ → 0 i.e., the situation where susceptible individuals who have a high awareness level have very reduced likelihood of getting infected with tuberculosis due to their level of awareness and the benefits that accrue from such awareness. Clearly, effective control strategies for tuberculosis that utilizes strategies associated with the expression for A₃ will assist in reducing the TB burden in a community as can be seen from the results in Equation (12).

Case 4:

We will now discuss the situation where there are no disease controls in the system. In this worse case scenario, we will set r = ν = α₂ = π = α₁ = ψ = θ₁ = θ₂ = 0. This implies that, from model (3), (S₂, E₂, J, T) → 0 asymptotically, as t → ∞. Hence the basic reproduction number is obtained as

${{R}}_0={{R}}_T{|}_{r=\nu ={\alpha }_2=\pi ={\alpha }_1=\psi ={\theta }_1={\theta }_2=0}.$

So that ${R}$_T = A₄${R}$₀ where

$A_4=\frac{D_1}{D_2},$

with D₁ = ((d + μ)(μ + k₁)(πη + μ + r + ηνα₂)(k₁(μ(μ + σα₁p₂ + θ₁) + (μ + σα₁ + θ₁)(k₂ + θ₂)) + σk₂α₁(μ + ψ) + k₂(μ + θ₁)(μp₁ + ψ) + μ(σα₁p₂ + p₁(μ + θ₁))(μ + θ₂ + ψ))) and D₂ = [(k₁ + μp₁)(μ + r)(d + π + μ + να₂)(μ + α₁ + θ₁)[(μ + k₁)(μ + k₂ + θ₂) + (μ + k₂)ψ]].

If ${R}$₀ < 1, then TB cannot develop into an epidemic and no control strategy is needed for TB control. If we take the difference between ${R}$₀ and ${R}$_T i.e.,

${\Delta }_0:={{R}}_0-{{R}}_T=(1-A_4){{R}}_0$

then an effective combination of the controls discussed in this case will slow down the spread of TB in the population, which would mean that Δ₀ > 0, which is satisfied if A₄ < 1.

Applying the limits in Case 3 on A₄, we observe that the results are the same except for

$\begin{array}{cc}(i)\underset{\begin{array}{c}{\alpha }_1\to \infty \\ \psi \to \infty \\ {\theta }_1\to 0\\ {\theta }_2\to 0\end{array}}{\lim }{A}_4=\frac{(d+\mu )\sigma (\mu +k_1)(k_2+p_2\mu )(\pi \eta +\mu +r+\eta \nu {\alpha }_2)}{(\mu +k_2)(k_1+p_1\mu )(\mu +r)(d+\pi +\mu +\nu {\alpha }_2)}& \left(15\right)\end{array}$

and

$\begin{array}{cc}(ii)\underset{\begin{array}{c}\nu \to 1\\ \pi \to \infty \\ {\theta }_2\to 0\\ {\theta }_1\to 0\\ \psi \to \infty \\ {\alpha }_1\to \infty \\ {\alpha }_2\to \infty \end{array}}{\lim }{A}_4=\frac{\eta (d+\mu )\sigma (\mu +k_1)(k_2+p_2\mu )}{(\mu +k_2)(k_1+p_1\mu )(\mu +r)}& \left(16\right)\end{array}$

The difference is that the limits of A₄ are disease-induced death adjusted (compared to the limits of A₃) and, in Equations (15) and (16), the limiting results for A₄ will be zero if σ → 0 or η → 0 [η → 0 applies just for Equation (16)] i.e., susceptible individuals who have a high awareness level have very reduced likelihood of getting infected with tuberculosis due to their level of awareness (and the benefits that accrue therefrom) and infectious individuals detected for treatment also have a reduced likelihood of infecting susceptible and latently infected individuals with tuberculosis during the course of their treatment.

Of course, the above results (discussed in Cases 1 to 4, in this section) are dependent on the fact that the baseline populations are not subjected to the same interventions. For example, in Case 1, A₁ measures treatment effectiveness for a population where education, chronic cough identification and an active case-finding strategy are already implemented (without treatment) while in Case 4, A₄ measures the effort to reduce the reproduction number for a population that had no intervention whatsoever; hence in the latter case (Case 4), the effort to curtail the disease must be higher than in the former (Case 1).

In summary, parameter values that would make A₁ < 1, A₂ < 1, A₃ < 1 or A₄ < 1 could yield control strategies with the capacity of reducing the number of secondary TB infections and slow down the spread of tuberculosis in the population. Hence, to determine the necessary condition for slowing down the development of TB at the population level, we will require that Δ_F > 0, Δ_B > 0, Δ_A > 0 or Δ₀ > 0, as the case may be and, appropriately, determine control scenarios that will give better results in reducing the incidence (and prevalence) of tuberculosis in the population.

It is imperative to state that there are several limits that could be considered while studying A₁, A₂, A₃ and A₄. However, our interests are in investigating control strategies involving the combination of the following parameters: ν, π, ψ, θ₁, θ₂, α₂, r, α₁ and η. The combined effect of these parameters on control measures and their effect on on the dynamics of tuberculosis is worth examining and tested in improving the case detection rate as well as reduce the tuberculosis burden in a population.

2.2.4. Analysis of Effective Reproduction Number Under Controls, ${R}$_T

Using the threshold parameter, ${R}$_T, we want to study the effect of treatment, active cough identification and the associated cost factor, tuberculosis awareness, an active case-finding technique as well as a combination of some of these factors on the dynamics of TB in the population.

It is evident from (4) that

$\begin{array}{cc}\underset{\begin{array}{c}{\alpha }_1\to \infty \\ \psi \to \infty \\ {\theta }_1\to 0\\ {\theta }_2\to 0\end{array}}{\lim }{{R}}_T=\frac{\beta \sigma (k_2+p_2\mu )(\pi \eta +\mu +r+\eta \nu {\alpha }_2)}{(k_2+\mu )(\mu +r)(d+\pi +\mu +\nu {\alpha }_2)}>0,& \left(17\right)\end{array}$

$\begin{array}{cc}\underset{\begin{array}{c}{\alpha }_1\to \infty \\ \psi \to \infty \\ {\theta }_1\to 0\\ {\theta }_2\to 0\\ r\to \infty \end{array}}{\lim }{{R}}_T=\frac{\beta \sigma (k_2+p_2\mu )}{(k_2+\mu )(d+\pi +\mu +\nu {\alpha }_2)}>0,& \left(18\right)\end{array}$

$\begin{array}{cc}\underset{\begin{array}{c}\nu \to 1\\ {\alpha }_2\to \infty \\ {\alpha }_1\to \infty \\ \psi \to \infty \\ {\theta }_1\to 0\\ {\theta }_2\to 0\\ r\to \infty \end{array}}{\lim }{{R}}_T=\underset{\begin{array}{c}{\alpha }_2\to \infty \\ \nu \to 1\\ \pi \to \infty \\ r\to \infty \end{array}}{\lim }{{R}}_T=0.& \left(19\right)\end{array}$

It seems, from the limits of ${R}$_T in Equations (17)–(19), that a high treatment rate is very effective in the control of tuberculosis only when the effect of the other critical parameters are properly harnessed e.g., high awareness rates, less loss of awareness (forgetfulness), increased active case-finding and chronic cough identification (with minimal costs involved in testing and treatments).

The contour plot of ${R}$_T, as a function of the active case-finding rate (π) and the treatment rate (r) is shown in Figure 2 while Figures 3, 4 depicts the contour plots of ${R}$_T as a function of the chronic cough identification rate (α₂) and the associated cost factor (ν) and as a function of the awareness rate for the susceptible group (α₁) and the chronic cough identification rate (α₂), respectively. Using the base parameter values in Table 1, we observe from Figures 2–4 that the indicated parameters have a positive effect on the effective reproduction number, ${R}$_T. Hence, improving on the awareness rate for the susceptible individuals, the associated cost factor, active case-finding strategy, the use of chronic cough as a TB marker (for likely cases) and effective treatment, will reduce the value of ${R}$_T, albeit at different rates.

FIGURE 2

Figure 2. Contour plot of ${R}$_T as a function of π and r.

FIGURE 3

Figure 3. Contour plot of ${R}$_T as a function of ν and α₂.

FIGURE 4

Figure 4. Contour plot of ${R}$_T as a function of α₁ and α₂.

TABLE 1

Table 1. Parameter information.

From the expressions in the limits of ${R}$_T in Equation (19), it is seen that near total eradication of tuberculosis is achievable. One strategy is to effectively combine making tuberculosis medical tests and treatment free, having high awareness rates (for both susceptible and latently infected persons), continuous enlightenment programmes to reduce the likelihood of loss of awareness (or forgetfulness), high reportage of chronic cough (to improve the detection of likely TB cases) and effective and high treatment rate (i.e., ν → 1, α₁ → ∞(ψ → ∞), θ₁ → 0(θ₂ → 0), α₂ → ∞, r → ∞). It is remarkable that a similar conclusion is reached when a control strategy concentrates on high chronic cough identification rate, making tuberculosis tests and treatment free, high active case-finding rate and a high and effective treatment rate (i.e., ν → 1, α₂ → ∞, r → ∞, π → ∞).

Of course, the above mentioned control strategies could be cost prohibitive and may seem unrealistic, especially in developing countries. However, the alternative strategies that can be gleaned from the expressions in Equations (17) and (18) can be applied to the target community, in reducing the severity of tuberculosis in the community; note that, for example, the limit of ${R}$_T in Equation (17) does not necessarily involve having a high treatment rate especially when such treatment level is lacking in the community.

Computing the partial derivatives of ${R}$_T with respect to the key parameters (r₂, α₂, α₁, θ₁, θ₂, ψ, π and ν) further reveals the effect of these parameters on tuberculosis control in the population. The derivatives are given in Appendix B.

Clearly, it follows from Equation (A3) (Appendix B) that $\frac{\partial {{R}}_T}{\partial r}<0$, unconditionally. Therefore, effective treatment rates of tuberculosis will have a positive impact in reducing the disease burden in the population. This result is stated in the following lemma

Lemma 2.2. Effective treatment will have a positive impact on the TB burden in a community by reducing the incidence of the disease in the population regardless of the values of the other parameters in the expression for the effective reproduction number.

Also, from Appendix B i.e., Equations (A4), (A5), we see that $\frac{\partial {{R}}_T}{\partial {\alpha }_1}>0$ $(\frac{\partial {{R}}_T}{\partial {\theta }_1}<0)$ if

$\begin{array}{cc}\sigma >{\sigma }^{\ast }=\frac{k_1{k}_2+k_1{\theta }_2+k_2\psi +\mu [k_1+k_2{p}_1+p_1(\mu +{\theta }_2+\psi )]}{k_1{k}_2+k_1{\theta }_2+k_2\psi +\mu [k_1{p}_2+k_2+p_2(\mu +{\theta }_2+\psi )]}.& \left(20\right)\end{array}$

However, $\frac{\partial {{R}}_T}{\partial {\alpha }_1}<0$ $(\frac{\partial {{R}}_T}{\partial {\theta }_1}>0)$ if σ < σ* and $\frac{\partial {{R}}_T}{\partial {\alpha }_1}=0$ $(\frac{\partial {{R}}_T}{\partial {\theta }_1}=0)$ when σ = σ*. This implies that a high awareness level by some of the susceptible individuals will have a positive impact on the TB burden in the population if σ < σ* i.e., the likelihood of such susceptible individuals getting infected with TB should be less than the computed σ*. The strategy of using awareness (and its efficiency on TB control) to combat the disease will fail to reduce the burden of tuberculosis in a community if σ = σ*, and, in fact, will have a detrimental impact on the community, by increasing the value of ${R}$_T, if σ > σ*. This result is reversed if there is a loss of awareness over time by susceptible individuals who previously had a high level of awareness (forgetfulness). Hence, such loss of awareness will bring about a reduction in the value of the effective reproduction number if the likelihood of infection of these susceptible (with previously high level of awareness) is greater than the computed σ*, it will have no impact if σ = σ* and will have a detrimental effect on the community if σ < σ*. The result is summarized thus:

Lemma 2.3. A high awareness level for susceptible individuals will have a positive impact on the reduction of the TB burden in a community if σ < σ*, no impact if σ = σ* and a detrimental impact if σ > σ*. The result is reversed for the rate of loss of awareness by these susceptibles: a positive impact if σ > σ*, no impact when σ = σ* and a detrimental impact if σ < σ*.

This result highlights the need for sustaining the tuberculosis awareness programmes so that forgetfulness could be minimized and the power of the knowledge of TB could help in reducing the likelihood of new tuberculosis infections in the population.

A very close look at the expression for σ*, giving in Equation (20), shows the significance of the fraction of fast progressions in the system i.e., p₁ and p₂. Clearly, if p₁ = p₂ = 1, then σ* = 1 (which is also the case when p₁ = p₂ = p and k₁ = k₂ = k); recall that one major assumption in the model (3) is that, due to the effectiveness of awareness, σ ≤ 1. Therefore, we conjecture that these fractions of fast progressions, p₁ and p₂ and the rate of endogenous reactivation, k₁ and k₂ are very crucial in determining whether σ* is less than, equal to or greater than one. Hence, the impact of awareness or loss of it over time, on the dynamics of tuberculosis, is tied to the critical value of the reduced likelihood of susceptible individuals, with high awareness level, getting infected with TB (σ) and to the fractions of fast TB progressions from new infections and the rates of endogenous reactivations (based on the levels of awareness in the population).

Taking a look at Appendix B, i.e., from Equations (A6), (A7), we observe that $\frac{\partial {{R}}_T}{\partial {\theta }_2}<0$ $(\frac{\partial {{R}}_T}{\partial \psi }>0)$ if

$\begin{array}{cc}{k}_1<k_2.& \left(21\right)\end{array}$

This implies that the high awareness level of some latently infected individuals will bring about a reduction in the value of the effective reproduction number (a positive impact) if the progression rate (endogenous reactivation), k₁, of latent individuals with low awareness level, to the active stage of TB is greater than that for latently infected individuals with high awareness level i.e., k₁ > k₂; there will be no impact on TB dynamics if k₁ = k₂ and a detrimental impact when k₁ < k₂. This result is reversed when we investigate the relationship between the effective reproduction number and the rate of loss of awareness by the latently infected individuals who previously had a high awareness level (θ₂). This leads to the following lemma:

Lemma 2.4. A high awareness rate for some of the latently infected individuals will have a positive impact on the reduction of the TB burden in a community if k₁ > k₂, no impact if k₁ = k₂ and a detrimental impact if k₁ < k₂. The result is reversed for the rate of loss of awareness by such latently infected: a positive impact if k₁ < k₂, no impact when k₁ = k₂ and a detrimental impact if k₁ > k₂.

This result significantly demonstrates the connection between awareness levels (and loss of awareness over time) of latently infected individuals and their progression rates (endogenous reactivation) to active tuberculosis, on the dynamics of the disease in a community. We therefore conjecture that, from Lemmas 2.3 and 2.4, the effect of the awareness levels (and loss of previously high awareness level) of susceptible and latently infected individuals is tied to some critical parameters: for the former (susceptible), the effect of high awareness level (and lack of it) on the effective reproduction number depends heavily on σ (which further depends on p₁ and p₂, and probably k₁ and k₂) while for the latter (latently infected individuals), the effect of high awareness level (and lack of it) on the effective reproduction number depends only on the progression rates to active TB, k₁ and k₂.

Basically, these Lemmas (2.3 and 2.4) demonstrates that, even when awareness levels wanes probably due to stoppages in enlightenment programmes in several media, then the progression rates (p₁, p₂, k₁ and k₂) plays a significant role in the dynamics of tuberculosis in the population and their effect should be harnessed by health officials for proper control of the disease.

From the analysis and results stated in Lemmas 2.3 and 2.4, it seems that the progression and reactivation rates (p₁, p₂, k₁, and k₂) can vary from individual to individual or from community to community. Generally, tuberculosis is most likely to occur in the first year following infection, with stepwise reduction year on year over the following 5–10 years, by which time incidence approaches that of uninfected contacts (Esmail et al., 2014). Also, reactivation several decades after initial infection occurs, but beyond 10 years, it is difficult to assess how common reactivation is (Esmail et al., 2014). Also, we earlier stated that the period of fast latency can span between 1 and 5 years (Styblo, 1980; Flynn and Chan, 2001; Colijn et al., 2007). Hence, the time line for the risk of disease over time is, generally, not fixed, as this can be determined by several biological and environmental factors, so that we can assume that the parameters that describe disease progressions can vary within some reasonable bounds.

Looking at Appendix B, i.e., Equations (A8), (A9), and (A10), we observe that $\frac{\partial {{R}}_T}{\partial \pi }<0$, $\frac{\partial {{R}}_T}{\partial {\alpha }_2}<0$, and $\frac{\partial {{R}}_T}{\partial \nu }<0$ if

$\begin{array}{cc}d<d^{\ast }=\frac{(1-\eta )\mu +r}{\eta },\eta \ne 0& \left(22\right)\end{array}$

or

$\begin{array}{cc}\eta <{\eta }^{\ast }=\frac{\mu +r}{\mu +d}& \left(23\right)\end{array}$

The results in Equations (22) and (23) reveals that the active case-finding strategy, the use of chronic cough as a marker for identifying potential TB cases, together with an efficient cost factor, will bring about a reduction in the value of the effective reproduction number if the disease-induced death rate is less than the computed d* (22) or for the relative infectiousness of individuals with active TB who are detected for immediate treatment to be less than the computed η* (23). The aforementioned parameters will have no impact on TB control if d = d* or η = η* but will have a detrimental effect on tuberculosis control and the TB burden in the population if d > d* or η > η*. We can state the following lemma:

Lemma 2.5. A high active case-finding and chronic cough identification rates with a low cost factor will have a positive impact on the reduction of the TB burden in a community if d < d*(η < η*), no impact if d = d*(η = η*) and a detrimental impact if d > d*(η > η*).

Recall that we assume that the modification parameter η ≤ 1. Therefore, we expect that η* ≤ 1, which implies that r ≤ d. Hence, if the disease-induced death rate is greater than the treatment rate, then the active case-finding strategy, coupled with the effective use of chronic cough identification and an effective cost factor will have a positive impact on the dynamics of tuberculosis and bring down the TB burden in the population.

From the analysis carried out on ${R}$_T, one observes that taking treatment alone without considering the effect of other parameters on the dynamics of tuberculosis, may not be enough in reducing the TB burden in the population. Clearly, a critical combination of two or more parameters from the set of key parameters, notably ν, α₁, α₂ and others like θ₁, π, ψ and θ₂, can significantly affect the value of the effective reproduction number. These parameters also affect the number of detected cases (which improves the case detection and notification rates) and, with effective treatment, will reduce the tuberculosis burden in the population.

2.2.5. Backward Bifurcation Analysis: Special Case

The phenomenon of backward bifurcation, which has been observed in several disease transmission models (see, e.g., Hadeler and van den Driessche, 1997; Castillo-Chavez and Song, 2004), is typically characterized by the coexistence of a stable DFE and a stable endemic equilibrium when the associated reproduction number of the model is less than unity. The public health implication of the backward bifurcation phenomenon of model (3) is that the classical epidemiological requirement of having the reproduction number (${R}$_T) be less than one, although necessary, is no longer sufficient for effective control of the disease in the population.

We can determine the possibility of the existence of a backward bifurcation in model (3); this is possible when we check for the existence of multiple endemic equilibria when the reproduction number of model (3) is less than one.

Consider the case when there is no treatment in the model (3) (this is fitting in communities where treatment facilities are not available or not enough to cater for the affected individuals) and insignificant disease-induced death rate i.e., ϵ = r = d = 0. It should be noted that setting d = 0 in (3) gives $N(t)\to \frac{\Lambda }{\mu }$ as t → ∞. Let $\widehat{\beta }=\frac{\mu \beta }{\Lambda }$ so that the force of infection now becomes

$\begin{array}{cc}\lambda =\widehat{\beta }(I+\eta J).& \left(24\right)\end{array}$

Also, let ${\widehat{{R}}}_T$ be the effective reproduction number of model (3) with ϵ = r = d = 0.

To find the conditions for the existence of the endemic equilibrium for the model (3) with ϵ = r = d = 0, denoted by ${\xi }_2=(S_1^{**},S_2^{**},E_1^{**},E_2^{**},I^{**},J^{**})$, the Equation in (3), with ϵ = r = d = 0, are solved in terms of the force of infection, at steady state [using Equation (24)] (where the force of infection at steady state will be denoted by λ^e), and this must satisfy the following polynomial:

$\begin{array}{cc}f({\lambda }^e)=A_1{\lambda }^e{}^4+A_2{\lambda }^e{}^3+A_3{\lambda }^e{}^2+A_4{\lambda }^e+A_5=0& \left(25\right)\end{array}$

where the coefficients, A₁, …, A₅ and the endemic equilibrium, ξ₂, are given in Appendix C.

The components of the endemic equilibrium, given in Equation (A11) (Appendix C), are then obtained by solving for λ^e from the quartic (25), and substituting the positive values of λ^e into the expressions of the endemic steady states given in Equation (A11). Furthermore, it follows that the coefficient A₁, of the quartic (25), is always positive, and A₅ is positive (negative) if ${\widehat{{R}}}_T$ is less (greater) than one. The following can be deduced:

Theorem 2.3. The treatment-free model of (3), with d = 0, has

i no endemic equilibrium if b₁ = 0 and ${\widehat{{R}}}_T<1$ (absence of exogenous re-infection from the E₁ class).

ii has four or two endemic equilibria if A₂ < 0, A₃ > 0, A₄ < 0 and ${\widehat{{R}}}_T<1$.

iii has two endemic equilibria if A₃ is of the same sign as A₂ or A₄ and ${\widehat{{R}}}_T<1$.

iv has two endemic equilibria if A₂ > 0, A₃ < 0, A₄ > 0 and ${\widehat{{R}}}_T<1$.

v no endemic equilibrium otherwise when ${\widehat{{R}}}_T<1$.

vi a unique endemic equilibrium when ${\widehat{{R}}}_T>1$.

Items (ii)—(iv) of Theorem 2.3 suggests the possibility of backward bifurcation in the treatment-free model, with d = 0. Determining specific parameters that could be the cause of the phenomenon, for the system (3), is quite challenging due to the number of parameters in the model and the degree of the polynomial (25). However, since it has been established that exogenous re-infection can trigger the backward bifurcation phenomenon in the dynamics of tuberculosis (Hadeler and van den Driessche, 1997; Castillo-Chavez and Song, 2004; Okuonghae and Omosigho, 2011), it will be interesting to observe the effects of both exogenous re-infection parameters in the system i.e., b₁ and b₂, since these parameters are based on the levels of awareness of the latently infected individuals.

First of all, it is important to show that in the absence of exogenous re-infection from the latently infected class of individuals with low awareness level, b₁ = 0, the system (3) will not have a backward bifurcation when ${\widehat{{R}}}_T<1$. Clearly, setting b₁ = 0 in Equation (25), yields the quadratic equation

$\begin{array}{cc}f({\lambda }^e)={\widehat{A}}_3{\lambda }^e{}^2+{\widehat{A}}_4{\lambda }^e+A_5=0& \left(26\right)\end{array}$

where Â₃ and Â₄ are now evaluated from A₃ and A₄, respectively, with b₁ = 0. It is easy to see that Â₃ > 0 and A₅ > 0 (the latter when ${\widehat{{R}}}_T<1$). We want to show that, since σ ≤ 1, the coefficient Â₄ > 0, so that the quadratic (26) will not have any positive roots, which implies that there is no endemic equilibrium when b₁ = 0 in Equation (25). With b₁ = 0 in Equation (25), we see that

$\begin{array}{l}{\widehat{A}}_4=\Lambda [(\mu +\pi +\nu {\alpha }_2)(\mu +\mu \sigma +\sigma {\alpha }_1+{\theta }_1)((\mu +k_1)\hfill \\ (\mu +k_2+{\theta }_2)+(\mu +k_2)\psi )]-\Lambda [\beta (\mu +\eta (\pi +\nu {\alpha }_2))\hfill \\ \sigma ((k_1+\mu p_1)(\mu +k_2+{\theta }_2)+(k_2+\mu p_1)\psi )]\hfill \end{array}$

Since σ ≤ 1, it then implies that,

$\begin{array}{cc}\begin{array}{l}{\widehat{A}}_4\ge \Lambda [(\mu +\pi +\nu {\alpha }_2)(\mu +\mu \sigma +\sigma {\alpha }_1+{\theta }_1)((\mu +k_1)\hfill \\ (\mu +k_2+{\theta }_2)+(\mu +k_2){\psi }_1)]-\Lambda [\beta (\mu +\eta (\pi +\nu {\alpha }_2))\hfill \\ ((k_1+\mu p_1)(\mu +k_2+{\theta }_2)+(k_2+\mu p_1){\psi }_1)].\hfill \end{array}\hfill & \left(27\right)\hfill \end{array}$

From ${\widehat{{R}}}_T<1$, it follows that

$\widehat{\beta }\le \frac{G_4}{(\mu +\eta \pi +\eta \nu {\alpha }_2)(G_1+G_2+G_3)},$

where $\widehat{\beta }$ is now evaluated with r = d = 0. Substituting the expression for $\widehat{\beta }$ into Equation (27), we see that

$\begin{array}{cc}\begin{array}{l}{\widehat{A}}_4\ge \Lambda [(\mu +\pi +\nu {\alpha }_2)(\mu +\mu \sigma +\sigma {\alpha }_1+{\theta }_1)((\mu +k_1)\hfill \\ (\mu +k_2+{\theta }_2)+(\mu +k_2){\psi }_1)]\hfill \\ -\Lambda [(\frac{G_4}{(\mu +\eta \pi +r_2+\eta \nu {\alpha }_2)(G_1+G_2+G_3)})(\mu +\eta \hfill \\ (\pi +\nu {\alpha }_2))((k_1+\mu p_1)(\mu +k_2+{\theta }_2)+(k_2+\mu p_1){\psi }_1)].\hfill \end{array}\hfill & \left(28\right)\hfill \end{array}$

Simplifying this further, we have that

$\begin{array}{l}{\widehat{A}}_4\ge \Lambda [(\mu +\pi +\nu {\alpha }_2)[((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi )\hfill \\ (k_1({\mu }^3+k_2{\left(\mu +{\alpha }_1+{\theta }_1\right)}^2+\mu (2\mu +{\alpha }_1+{\theta }_1)(p_2{\alpha }_1+{\theta }_1)\hfill \\ +{\left(\mu +{\alpha }_1+{\theta }_1\right)}^2{\theta }_2)+\mu {\alpha }_1{p}_2(2\mu +{\alpha }_1+{\theta }_1)(\mu +{\theta }_2+\psi )\hfill \\ +\mu p_1({\mu }^2+{\theta }_1(2\mu +{\alpha }_1+{\theta }_1))(\mu +k_2+{\theta }_2+\psi )\hfill \\ +k_2(\mu {\alpha }_1(2\mu +{\alpha }_1+{\theta }_1)+{\left(\mu +{\alpha }_1+{\theta }_1\right)}^2\psi ))]>0,\hfill \end{array}$

for σ ≤ 1. Therefore, the polynomial (26) will not have a positive root, hence no endemic equilibrium, when b₁ = 0 and ${\widehat{{R}}}_T<1$, ruling out the possibility of a backward bifurcation, for this case.

However, when b₂ = 0 and b₁ ≠ 0 (in this case, there is no exogenous re-infection of latently infected individuals with high awareness level but there is exogenous re-infection in the E₁ class i.e., latently infected individuals with low level of awareness), the quartic (25) now reduces to a cubic equation

$\begin{array}{cc}f({\lambda }^e)={\tilde{A}}_2{\lambda }^e{}^3+{\tilde{A}}_3{\lambda }^e{}^2+{\tilde{A}}_4{\lambda }^e+A_5=0& \left(29\right)\end{array}$

where $\widetilde{A_2},\widetilde{A_3}$ and $\widetilde{A_4}$ are obtained from evaluating the corresponding coefficients of (25) with b₂ = 0.

Following the analysis above, it is easy to show that the polynomial (29) can have more than one positive root (hence more than one endemic equilibria) when ${\widehat{{R}}}_T$ < 1, suggesting the presence of a backward bifurcation for the case b₂ = 0 and b₁ ≠ 0. We can now summarise the above discussions thus:

Lemma 2.6. The treatment-free model of (3) with d = 0 can undergo the backward bifurcation phenomenon when

i both exogenous re-infection parameters, b₁ and b₂, are present in the system i.e., b₁ ≠ 0 and b₂ ≠ 0.

ii only the exogenous re-infection parameter b₁ is present in the system i.e., b₁ ≠ 0 and b₂ = 0

The treatment-free model of (3) with d = 0 cannot undergo the backward bifurcation phenomenon when

i both exogenous re-infection parameters, b₁ and b₂, are absent in the system i.e., b₁ = 0 and b₂ = 0.

ii the exogenous re-infection parameter b₁ is absent in the system (even when b₂ ≠ 0) i.e., b₁ = 0 and b₂ ≠ 0.

This study has, to the best of our knowledge, showed for the first time the critical relationship between TB awareness, heterogeneity in exogenous re-infection (by virtue of the level of awareness of the latently infected individuals) and the tuberculosis burden in a community. We can see, from Lemma 2.6, that exogenous re-infection of latently infected individuals with low awareness level will more negatively impact on TB control than incidences of exogenous re-infection in the latently infected group with high awareness level, due to the influence of the former in allowing for the existence of a backward bifurcation in the system, when the associated reproduction number is less than one. This can be seen from the backward bifurcation diagrams in Figures 5, 6; clearly, in Figure 5, the backward bifurcation range is larger than what we have in Figure 6. Interestingly, when we set b₁ = 0 (the exogenous re-infection parameter for the latently infected individuals with low awareness level), there is no backward bifurcation in the system when the associated reproduction number is less than unity, regardless of the value of b₂ (the exogenous re-infection parameter for the latently infected individuals with high awareness level).

FIGURE 5

Figure 5. Backward bifurcation when b₁ ≠ 0 and b₂ ≠ 0.

FIGURE 6

Figure 6. Backward bifurcation when b₁ ≠ 0 and b₂ = 0.

Characterizing the backward bifurcation phenomenon for the complete model (3) is not trivial. However, using the Center Manifold Theorem (Carr, 1981; Castillo-Chavez and Song, 2004; Okuonghae and Omosigho, 2011), we can show the condition required for the system (3) to undergo the backward bifurcation phenomenon when ${R}$_T < 1. We claim the following:

Theorem 2.4. The model (3) exhibits backward bifurcation at ${R}$_T = 1 whenever the bifurcation coefficients, denoted by a and b (and given by (A15) and (A16) in Appendix D, respectively), are positive. However, if the coefficient a is negative, then the system (3) will not undergo a backward bifurcation at ${R}$_T = 0.

See Appendix D for proof of Theorem 2.4.

Of course determining, by inspection, which parameter(s) could cause the phenomenon of backward bifurcation in the model Equation (3) [by checking for specific parameter(s) that will make the bifurcation coefficients a, given in Equation (A15), to be negative and b, given in Equation (A16), to be positive] is non-trivial.

However, we conjecture that, in addition the exogenous re-infection parameters, b₁ and b₂, there could be other parameters that could determine whether the system (3) can undergo a backward bifurcation at ${R}$_T = 1. These parameters, we believe, are liked to the awareness parameters (i.e., α₁, ψ, θ₁, and θ₂) and the fractions of fast progressions, p₁ and p₂. Their overall effect on the backward bifurcation phenomenon for the model (3) is left for future work.

3. Results

Model (3) is now numerically simulated with the parameter estimates in Table 1 to gain insight into some of its quantitative features. The system of equations in model (3) was solved numerically using MATLAB; we used the ode45 solver, which is based on the Runge-Kutta method. The parameter values used were based on what exists in literature such as the authors previous works (Okuonghae and Omosigho, 2011) and other references cited in Okuonghae and Omosigho (2011). Few parameter values were assumed, as seen (Okuonghae and Omosigho, 2011) and other relevant literature cited in Okuonghae and Omosigho (2011). The numerical simulations are performed to illustrate various dynamical regimen characteristics by varying some of the key parameters in the model. The effects of varying some of the key parameters on the infected classes are presented and we will examine practicable preventive measures characterized by these variations. Parameter values that are different from those stated in Table 1 are shown in the caption of the respective figure.

Figures 7, 8 shows the effect of varying the awareness rate of the susceptible population on the infected classes, vis à viz increasing the active case-finding rate. Clearly, we observe that, even with loss of awareness on the part of susceptible individuals (who previously had a high awareness rate), increasing the active case-finding rate affected (positively) the dynamics of the system with a reduction in the proportion of individuals in the infected classes, and not just the infectious class (I), only.

FIGURE 7

Figure 7. Simulations of model (3) showing the number of infected individuals. Here, the awareness rate for the susceptible individuals (α₁) is varied from 0 to 30 with θ₁ = 20, θ₂ = 1and π = 5.

FIGURE 8

Figure 8. Simulations of model (3) showing the number of infected individuals. Here, the awareness rate for the susceptible individuals (α₁) is varied from 0 to 30 with θ₁ = 20, θ₂ = 1 and π = 30.

Figure 9 shows that increasing the chronic cough identification rate (and by implication improving the case detection rate) significantly reduced the proportion of infected individuals, even with a “high” loss of awareness in the susceptible group while Figure 10 shows a worse case scenario whereby an increase in loss of awareness in the susceptible group leads to an increase in the proportion of infected individuals.

FIGURE 9

Figure 9. Simulations of model (3) showing the number of infected individuals. Here, the cough identification rate (α₂) is varied from 0 to 30 with θ₁ = 20 and θ₂ = 1.

FIGURE 10

Figure 10. Simulations of model (3) showing the number of infected individuals. Here, the rate at which susceptibles lose awareness (θ₁) is varied from 0 to 30 with α₁ = 30 and α₂ = 5.

Figures 11, 12 shows the proportion of individuals in the infected classes when we vary the reduced likelihood of infection by detected infectious TB cases receiving treatment, η, between 0 and 1. Of course, as expected, reducing the value of η, from 1 to 0, brings down the proportion of infected individuals, albeit at different rates, depending on the values of the awareness parameters for the susceptible and latently infected populations.

FIGURE 11

Figure 11. Simulations of model (3) showing the number of infected individuals. Here, the disease transmission modification parameter (η) is varied from 0 to 1 with α₁ = α₂ = 5, θ₁ = θ₂ = 30 and ψ = 5.

FIGURE 12

Figure 12. Simulations of model (3) showing the number of infected individuals. Here, the disease transmission modification parameter (η) is varied from 0 to 1 with α₁ = α₂ = 5, θ₁ = θ₂ = 1 and ψ = 30.

Remarkably, Figure 13 shows that there is little influence of the cost factor (ν) on the proportion of infected individuals, in the presence of awareness, chronic cough identification, less loss of awareness and an impressive awareness rate for the latently infected individuals and active case-finding rate.

FIGURE 13

Figure 13. Simulations of model (3) showing the number of infected individuals. Here, the cost factor (ν) is varied from 0 to 1 with α₁ = α₂ = 5, θ₁ = θ₂ = 1, ψ = 30 and π = 20.

Taking a look at Figures 14, 15 reveals the effect of awareness in preventing tuberculosis infection amongst the group of susceptibles with high awareness level as we vary the infection modification parameter σ, from 0 to 1. Of course, a powerful awareness campaign with high awareness rate and with σ = 0, brings about a drop in the proportion of infected individuals especially with an impressive awareness level (e.g., ψ = 30) and less loss of awareness (θ₁ = θ₂ = 1).

FIGURE 14

Figure 14. Simulations of model (3) showing the number of infected individuals. Here, the reduced likelihood of infection due to awareness (σ) is varied from 0 to 1 with α₁ = α₂ = 5, θ₁ = θ₂ = 1 and ψ = 30.

FIGURE 15

Figure 15. Simulations of model (3) showing the number of infected individuals. Here, the reduced likelihood of infection due to awareness (σ) is varied from 0 to 1 with α₁ = α₂ = 20, θ₁ = θ₂ = 20 and ψ = 30.

Clearly, focusing on awareness programmes (for both susceptible and latently infected individuals) is very beneficial for TB control. This should be taken in collaboration with other interventions, like the use of chronic cough as a marker for identifying potential TB cases, improved active case-finding strategy and reduced cost of disease management as well as prolonged awareness programmes to prevent loss of awareness over time.

It is important to state here that the parameter values used in the simulations can be sensitive to the model (3) and there could be uncertainties in their values. However, the use of these parameter values is to demonstrate their effect on the system and gain some insight into the quantitative (especially asymptotic) behaviors of the model.

Remark: Note the difference in the time window for the different figures, some having 10, 20, or 100 years of simulation. The graphs having a 10 or 20 year time window was used to observe the pattern of the dynamics of the diseases for the first few years of simulations or disease outbreak. We observe that the results could easily be drowned out if the graphs should have a longer time window. However, the figures with 100 year time window was used to investigate the existence of equilibria based on the parameter values used. It is important to state that tuberculosis dynamics are slow, and consequently, TB epidemics unfold over several decades (Aparicio and Castillo-Chavez, 2009). Hence determining asymptotic behaviors of the disease model could require running simulations spanning a long period of time.

4. Discussions and Conclusions

A modified and realistic deterministic mathematical model for the transmission dynamics of tuberculosis in a population has been designed and mathematically analyzed. The model (3) extends that formulated and studied in Okuonghae and Omosigho (2011) by classifying both susceptible and latently infected individuals by their level of awareness of tuberculosis (what the disease is, signs and symptoms of tuberculosis and government policy on testing and medical treatment) and included an active case-finding parameter in addition to the use of chronic cough as a marker for identifying potential TB cases, for the control of tuberculosis in a population.

The model (3) has a locally asymptotically stable DFE whenever the associated effective reproduction number is less than one. However, for the treatment free model with insignificant disease-induced death rate, it was shown that the system will undergo the phenomenon of backward bifurcation where the stable disease-free equilibrium will coexist with a stable endemic equilibrium when the effective reproduction number is less than one. For this special case, this phenomenon was caused by the exogenous re-infection of latently infected individuals. We then showed that the effect of exogenous re-infection on the backward bifurcation phenomenon depended strongly on the level of awareness of the latently infected individuals; backward bifurcation due to exogenous re-infection was largely sustained by the latently infected individuals with a low awareness level.

The study further showed that concentrating on TB treatment alone may not significantly reduce the value of the reproduction number if attention is not paid to other key control parameters such as awareness levels and active case-finding rate; in fact, it is possible for the value of the reproduction number to still be greater than one when we concentrate mainly on treatment rate. However, if we take two or more key parameters at the same time and glean out effective control strategies from their combination, then it is possible that, in the long run, the disease burden in the community can be reduced.

Qualitative study of the effective reproduction number showed how different control scenarios involving awareness levels, loss of awareness over time, active case-finding strategy, chronic cough identification and minimal (or no) cost incurred for TB management could lead to a reduction in the value of the effective reproduction number ${R}$_T. The analyses suggested that concentrating on increasing tuberculosis awareness campaign for both susceptible and latently infected individuals and also increasing chronic cough identification and active case-finding rates will result in reducing the incidence of TB in the population in the presence of an effective cost factor whereby the cost of conducting TB tests and commencing treatment is very small, if not free.

Numerical simulations suggested practical preventive measures, represented by changing the value of some parameters, notably, ν, α₁, α₂, θ₁, θ₂, ψ, and π. These simulations showed that improving on the cost factor, increased awareness programmes (especially as it affects susceptible and latently infected individuals), cough identification rate as well as minimizing new TB infections caused by fairly isolated TB infectious individuals that are undergoing treatment, will help in reducing the prevalence of TB in the community, together with minimal loss of awareness from both susceptible and latently infected individuals.

In summary, this work have shown, that the prospect of effectively controlling the spread of tuberculosis in a population is very bright. Preventive measures through the use of control strategies should concentrate on awareness programmes. The enlightenment programmes should include helping the general public make use of simple signs in quickly identifying a likely TB case such as chronic cough lasting at least 2 weeks. This will not only quicken the treatment of the infectious individuals, it can also reduce the likelihood of disease transmission. Also awareness programmes should be sustained over a long period of time (especially in places with high endemic levels of tuberculosis) as this work has demonstrated the effect of loss of awareness on the dynamics of tuberculosis in the population.

Author Contributions

Both authors were involved in the discussions and formulation of the mathematical model. BI was involved in the derivation of the effective reproduction number and the backward bifurcation analysis, for the special case considered. DO carried out the analysis of the effective reproduction number, bifurcation analysis for the complete model, as well as the numerical simulations in the manuscript. Both authors were involved in the write-up as well as proof reading of the final draft of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

DO will like to thank Prof. Denise Kirschner for the meaningful discussions, during the 45th Union World Conference on Lung Health, held in Barcelona, Spain in 2014, that led to the extensions of the basic model resulting in the formulation and analysis of the model presented in this article.

References

Aparicio, J. P., and Castillo-Chavez, C. (2009). Mathematical modelling of tuberculosis epidemics. Math. Biosci. Eng. 6, 209–237. doi: 10.3934/mbe.2009.6.209

PubMed Abstract | CrossRef Full Text | Google Scholar

Armijos, J. W., Weigel, M. M., Qincha, M., and Ulloa, B. (2008). The meaning and consequences of tuberculosis for an at-risk urban group in Ecuador. Rev. Panam. Salud. Public. 23, 188–197. doi: 10.1590/S1020-49892008000300006

PubMed Abstract | CrossRef Full Text | Google Scholar

Baral, S. C., Karki D. K., and Newell, J. N. (2007). Causes of stigma and discrimination accosiated with tuberculosis in Nepal: a qualitative study. BMC Public Health 7:211. doi: 10.1186/1471-2458-7-211

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennstam, A. L., Strandmark, M., and Diwan, V. K. (2004). Perception of tuberculosis in the Denogratic Republic of Congo: wali ya nkumu in the Mai Ndombe District. Qual. Health Res. 14, 299–312. doi: 10.1177/1049732303261822

PubMed Abstract | CrossRef Full Text | Google Scholar

Carr, J. (1981). Applications of Centre Manifold Theory. New York, NY: Springer-Verlag.

Google Scholar

Castillo-Chavez, C., and Song, B. (2004). Dynamical models of tuberculosis and their applications. Math. Biosci. Eng. 2, 361–404. doi: 10.3934/mbe.2004.1.361

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, S.-H., and Cataldo J. K. (2014). A systematic review of global cultural variations in knowledge, attitudes and health responses to tuberculosis stigma. Int. J. Tuberc. Lung. Dis. 18, 168–173. doi: 10.5588/ijtld.13.0181

PubMed Abstract | CrossRef Full Text | Google Scholar

Colijn, C., Cohen, T., and Murray, M. (2007). Emergent heterogeneity in declining tuberculosis epidemics. J. Theor. Biol. 247, 765–774. doi: 10.1016/j.jtbi.2007.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Diekmann, O., Heesterbeek, J. A., and Metz, J. A. (1990). On the definition and the computation of the basic reproduction ratio R₀ in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28, 365–382. doi: 10.1007/BF00178324

PubMed Abstract | CrossRef Full Text | Google Scholar

Esmail, H., Barry, C. E., Young, D. B. III, and Wilkinson, R. J. (2014). The ongoing challenge of latent tuberculosis. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 369:20130437. doi: 10.1098/rstb.2013.0437

PubMed Abstract | CrossRef Full Text | Google Scholar

Flynn, J. L., and Chan, J. (2001). Tuberculosis: latency and reactivation. Infect. Immun. 69, 4195–4201. doi: 10.1128/IAI.69.7.4195-4201.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadeler, K. P., and van den Driessche, P. (1997). Backward bifurcation in epidemic control. Math. Bio. 146, 15-35. doi: 10.1016/S0025-5564(97)00027-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Hethcote, H. W. (2000). The mathematics of infectious diseases. SIAM Rev. 42, 599-653. doi: 10.1137/S0036144500371907

CrossRef Full Text | Google Scholar

Johansson, E., Long, N. H., Diwan, V. K., and Winkvist, A. (2000). Gender and tuberculosis control: perspectives on health-seeking behaviour among men and women in Viet Nam. Health Policy 52, 33–51. doi: 10.1016/S0168-8510(00)00062-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakshmikantham, S., Leela, S., and Martynyuk, A. A. (1989). Stability Analysis of Nonlinear Systems. New York; Basel: Marcel Dekker, Inc.

Google Scholar

Mukandavire, Z., Bowa, K., and Garira, W. (2007). Modelling circumcision and condom use as HIV/AIDS preventive control strategies. Math. Comput. Modell. 46, 1353–1372. doi: 10.1016/j.mcm.2007.01.001

CrossRef Full Text | Google Scholar

Okuonghae, D. (2015). Lyapunov functions and global properties of some tuberculosis models. J. Appl. Maths Comput. 48, 421–439. doi: 10.1007/s12190-014-0811-4

CrossRef Full Text | Google Scholar

Okuonghae, D., and Omosigho, S. E. (2010). Determinants of TB case detection in nigeria: a survey. Global J. Health Sci. 2, 123–128. doi: 10.5539/gjhs.v2n2p123

CrossRef Full Text | Google Scholar

Okuonghae, D., and Omosigho, S. E. (2011). Analysis of a mathematical model for tuberculosis: what could be done to increase case detection. J. Theor. Biol. 269, 31–45. doi: 10.1016/j.jtbi.2010.09.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharomi, O., Podder, C. N., Gumel A. B., and Song, B. (2008). Mathematical analysis of the transmission dynamics of HIV/TB co-infection in the presence of treatment. Math. Biosci. Eng. 5, 145–174. doi: 10.3934/mbe.2008.5.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Styblo, K. (1980). Recent advances in epidemiological research in tuberculosis. Adv. Tuberc. Res. 20, 1–63.

PubMed Abstract | Google Scholar

van den Driessche, P., and Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Bios. 180, 29-48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

World Health Organization (2014). Global Tuberculosis Report. Geneva.

Appendix A

Proofs of Theorems 2.1 and 2.2

Proof of Theorem 2.1:

Proof. Let t₁ = sup{t > 0 : S₁(t) > 0, S₂(t) > 0, E₁(t) > 0, E₂(t) > 0, I(t) > 0, J(t) > 0, T(t) > 0} > 0. It follows from (3a) that

$\frac{d{S}_1}{dt}=\Lambda -{\alpha }_1{S}_1-\lambda S_1-\mu S_1+{\theta }_1{S}_2\ge \Lambda -{\alpha }_1{S}_1-\lambda S_1-\mu S_1,$

which can be written as

$\frac{d}{dt}\{S_1(t)exp[({\alpha }_1+\mu )t+{\displaystyle {\int }_0^t\lambda }(\tau )d\tau \left]\right\}\ge \Lambda \left\{exp\right[({\alpha }_1+\mu )t+{\displaystyle {\int }_0^t\lambda }(\tau )d\tau \left]\right\}.$

Thus,

$S_1(t_1)exp[({\alpha }_1+\mu )t_1+{\displaystyle {\int }_0^{t_1}\lambda }(\tau )d\tau ]-S_1(0)\ge {\displaystyle {\int }_0^{t_1}\Lambda }\left\{exp\right[({\alpha }_1+\mu )y+{\displaystyle {\int }_0^y\lambda }(\tau )d\tau \left]\right\}dy.$

so that,

$\begin{array}{l}{S}_1(t_1)\ge S_1(0)exp\left[-({\alpha }_1+\mu )t_1-{\displaystyle {\int }_0^{t_1}\lambda }(\tau )d\tau \right]+\hfill \\ -\left\{exp\left[-({\alpha }_1+\mu )t_1-{\displaystyle {\int }_0^{t_1}\lambda }(\tau )d\tau \right]\right\}\hfill \\ {\displaystyle {\int }_0^{t_1}\Lambda }\left\{exp\left[({\alpha }_1+\mu )y+{\displaystyle {\int }_0^y\lambda }(\tau )d\tau \right]\right\}dy>0.\hfill \end{array}$

Similarly, it can be shown that S₂(t) > 0, E₁(t) > 0, E₂(t) > 0, I(t) > 0, J(t) > 0 and T(t) > 0 for all time t > 0. Hence all solutions of the model (3) remain positive for all non-negative initial conditions, as required. □

Proof of Theorem 2.2:

Proof. Adding the equations of the model (3) gives

$\begin{array}{cc}\frac{dN}{dt}=\Lambda -\mu N-dI.& \left(\text{A}1\right)\end{array}$

Since the right hand side of the above equality is bounded by Λ − μN, a standard comparison theorem (Lakshmikantham et al., 1989) can be used to show that

$\begin{array}{cc}N(t)\le N(0)e^{-\mu t}+\frac{\Lambda }{\mu }(1-e^{-\mu t}).& \left(\text{A}2\right)\end{array}$

In particular, if $N(0)\le \frac{\Lambda }{\mu }$, then $N(t)\le \frac{\Lambda }{\mu }$. Thus, ${D}$ is a positively invariant set. Furthermore, if $N(0)\ge \frac{\Lambda }{\mu }$, then either the solution enters ${D}$ in finite time or N(t) approaches to $\frac{\Lambda }{\mu }$ as t → ∞. Hence no solution path leave through any boundary of ${D}$ and the region ${D}$ attracts all solutions in ${ℝ}_{+}^7$. □

Appendix B

Partial derivatives of ${R}$_T

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial r}=-\frac{\beta \eta (\pi +\nu {\alpha }_2)P_1}{({(\mu +r)}^2(d+\pi +\mu +\nu {\alpha }_2)(\mu +{\alpha }_1+{\theta }_1)((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ))},& \left(\text{A}3\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial {\alpha }_1}=\frac{\beta (\pi \eta +\mu +r+\eta \nu {\alpha }_2)(\mu +{\theta }_1)P_2}{((\mu +r)(d+\pi +\mu +\nu {\alpha }_2){(\mu +{\alpha }_1+{\theta }_1)}^2((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ))},& \left(\text{A}4\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial {\theta }_1}=-\frac{\beta {\alpha }_1(\pi \eta +\mu +r+\eta \nu {\alpha }_2)P_2}{((\mu +r)(d+\pi +\mu +\nu {\alpha }_2){(\mu +{\alpha }_1+{\theta }_1)}^2((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ))},& \left(\text{A}5\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial {\theta }_2}=\frac{\beta \mu (\pi \eta +\mu +r+\eta \nu {\alpha }_2)(k_1-k_2)[(1-p_1)(\mu +{\theta }_1)\psi +\sigma (1-p_2){\alpha }_1(\mu +k_1+\psi )]}{((\mu +r)(d+\pi +\mu +\nu {\alpha }_2)(\mu +{\alpha }_1+{\theta }_1){((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi )}^2)},& \left(\text{A}6\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial \psi }=\frac{\beta \mu (\pi \eta +\mu +r+\eta \nu {\alpha }_2)(k_1-k_2)[(1-p_1)(\mu +{\theta }_1)(\mu +k_2+{\theta }_2)+\sigma (1-p_2){\alpha }_1{\theta }_1]}{((\mu +r)(d+\pi +\mu +\nu {\alpha }_2)(\mu +{\alpha }_1+{\theta }_1){((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi )}^2)},& \left(\text{A}7\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial \pi }=\frac{\beta (\eta d-(1-\eta )\mu -r)P_1}{((\mu +r){(d+\pi +\mu +\nu {\alpha }_2)}^2(\mu +{\alpha }_1+{\theta }_1)((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ))},& \left(\text{A}8\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial {\alpha }_2}=\frac{\beta \nu (\eta d-(1-\eta )\mu -r)P_1}{((\mu +r){(d+\pi +\mu +\nu {\alpha }_2)}^2(\mu +{\alpha }_1+{\theta }_1)((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ))},& \left(\text{A}9\right)\end{array}$

$\begin{array}{cc}\frac{\partial {{R}}_T}{\partial \nu }=\frac{\beta {\alpha }_2(\eta d-(1-\eta )\mu -r)P_1}{((\mu +r){(d+\pi +\mu +\nu {\alpha }_2)}^2(\mu +{\alpha }_1+{\theta }_1)((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ))},& \left(\text{A}10\right)\end{array}$

where

P₁ = k₁(μ(μ + σα₁p₂ + θ₁ + (μ + σα₁ + θ₁)(k₂ + θ₂)) + σα₁k₂(μ + ψ) + k₂(μ + θ₁)(μp₁ + ψ) + μ(σα₁p₂ + p₁(μ + θ₁))(μ + θ₂ + ψ)),

and

P₂ = −[(1 − σ)(k₁k₂ + k₁θ₂ + k₂ψ) + k₁μ(1 − σp₂)] − (k₂μ(p₁ − σ) + μ(p₁ − σp₂)(μ + θ₂ + ψ).

Appendix C

Existence of EEP for Special Case

The coefficients of the polynomial (25) are given as

$\begin{array}{l}{A}_1=z^2\vartheta \sigma \Lambda (\mu +\pi +\nu {\alpha }_2),\hfill \\ A_2=\Lambda z(\mu +\pi +\nu {\alpha }_2)[\mu (\sigma +\vartheta (z+\sigma +z\sigma ))+\vartheta \sigma k_1+\sigma k_2+z\vartheta \sigma {\alpha }_1\hfill \\ +z\vartheta {\theta }_1+\sigma ({\theta }_2+\vartheta \psi )]-[\beta (\mu +\eta (\pi +\nu {\alpha }_2))z^2\vartheta \Lambda \sigma ]\hfill \\ A_3=\Lambda (\mu +\pi +\nu {\alpha }_2)[k_1(\mu (\sigma +z\vartheta (1+\sigma ))+\sigma k_2+z\vartheta (\sigma {\alpha }_1+{\theta }_1)+\sigma {\theta }_2)\hfill \\ +z{\theta }_1((1+\vartheta +z\vartheta )\mu +{\theta }_2+\vartheta \psi )+k_2(\mu (z+\sigma +z\sigma )+z\sigma {\alpha }_1+z{\theta }_1+\sigma \psi )\hfill \\ +z{\alpha }_1(\mu (\sigma +\vartheta (z+\sigma ))+\sigma {\theta }_2+\vartheta \sigma \psi )+\mu [\mu (\sigma +z(1+\sigma +\vartheta (1+z+\sigma )))\hfill \\ +(z+\sigma +z\sigma ){\theta }_2(\sigma +z\vartheta (1+\sigma ))\psi ]]-[\beta (\mu +\eta (\pi +\nu {\alpha }_2))\hfill \\ z\Lambda (z\vartheta \mu +\mu \sigma +\vartheta \sigma k_1+\sigma k_2+\vartheta \mu \sigma p_1+z\vartheta \sigma {\alpha }_1+z\vartheta {\theta }_1+\sigma {\theta }_2+\vartheta \sigma \psi )]\hfill \\ A_4=\Lambda (\mu +\pi +\nu {\alpha }_2)[k_1[k_2(\mu +\mu \sigma +\sigma {\alpha }_1+{\theta }_1)+{\theta }_1(\mu +z\vartheta \mu +{\theta }_2)+{\alpha }_1(\mu (z\vartheta +\sigma )+\sigma {\theta }_2)\hfill \\ +\mu (\mu (1+z\vartheta +\sigma )+(1+\sigma ){\theta }_2)]+k_2[{\theta }_1(\mu +z\mu +{\psi }_1)+{\alpha }_1(\mu (z+\sigma )+\sigma \psi )\hfill \\ +\mu (\mu (1+z+\sigma )+(1+\sigma )\psi )]+\mu [{\theta }_1((1+z+z\vartheta )\mu +(1+z){\theta }_2+(1+z\vartheta )\psi )\hfill \\ +{\alpha }_1(\mu (z+z\vartheta +\sigma )+(z+\sigma ){\theta }_2+(z\vartheta +\sigma )\psi )+\mu (\mu (1+z+z\vartheta +\sigma )\hfill \\ +(1+z+\sigma ){\theta }_2+(1+z\vartheta +\sigma )\psi )]\hfill \\ -\beta (\mu +\eta (\pi +\nu {\alpha }_2))\Lambda [k_1(\mu (z\vartheta +\sigma )+\sigma k_2+z\vartheta (\sigma {\alpha }_1+{\theta }_1)+\sigma {\theta }_2)\hfill \\ +k_2(\mu \sigma p_1+z(\mu +\sigma {\alpha }_1+{\theta }_1)+\sigma \psi )+\mu p_1(\mu (z\vartheta +\sigma )+z\vartheta {\theta }_1+\sigma ({\theta }_2+\psi ))\hfill \\ z((\mu +{\theta }_1)(\mu +{\theta }_2+\vartheta \psi )+\sigma {\alpha }_1(\mu n{p}_2+{\theta }_2+\vartheta (\mu +\psi )))]\hfill \\ A_5=\Lambda \mu (\mu +{\alpha }_1+{\theta }_1)(\mu +\pi +\nu {\alpha }_2)((\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi )[1-{\widehat{{R}}}_T]\hfill \end{array}$

where z = b₁ and ϑ = b₁b₂ with ${\widehat{{R}}}_T={{R}}_T{|}_{r=d=0}$ being the reproduction number of the model (3) where r = d = 0. The endemic steady state solutions are

$\begin{array}{cc}\begin{array}{l}{S}_1^{\ast \ast }=\frac{\Lambda (\mu +{\lambda }^e\sigma +{\theta }_1)}{{\alpha }_1(\mu +{\lambda }^e\sigma )+({\lambda }^e+\mu )(\mu +{\lambda }^e\sigma +{\theta }_1)},\hfill \\ S_2^{\ast \ast }=\frac{\Lambda {\alpha }_1}{{\alpha }_1(\mu +{\lambda }^e\sigma )+({\lambda }^e+\mu )(\mu +{\lambda }^e\sigma +{\theta }_1)},\hfill \\ E_1^{\ast \ast }=\frac{1}{b_1{\lambda }^e+(k_1+\mu +\psi )}[(1-p_1)S_1^{\ast \ast }{\lambda }^e+{\theta }_2{E}_2^{\ast \ast }],\hfill \\ E_2\ast \ast =\frac{1}{b_2{\lambda }^e+(k_2+\mu +{\theta }_2)}[(1-p_2)\sigma S_2^{\ast \ast }{\lambda }^e+\psi E_1^{\ast \ast }],\hfill \\ I^{\ast \ast }=\frac{1}{\nu {\alpha }_2+\mu +\pi }[p_1{S}_1^{\ast \ast }{\lambda }^e+p_2\sigma S_2^{\ast \ast }{\lambda }^e+(b_1{E}_1^{\ast \ast }+b_2{E}_2^{\ast \ast })+k_1{E}_1^{\ast \ast }+k_2{E}_2^{\ast \ast }],\hfill \\ J^{\ast \ast }=\frac{\pi +\nu {\alpha }_2}{\mu }{I}^{\ast \ast }\hfill \end{array}& \left(\text{A}11\right)\hfill \end{array}$

Appendix D

Proof of Theorem 2.4

Suppose

${\xi }_3=(S_1^{\ast \ast },S_2^{\ast \ast },E_1^{\ast \ast },E_2^{\ast \ast },I^{\ast \ast },J^{\ast \ast },T^{\ast \ast })$

represents any arbitrary endemic equilibrium of model (3) (i.e., an equilibrium in which at least one of the infected components is non-zero). We will explore the existence of backward bifurcation using the center manifold theory (Carr, 1981; Castillo-Chavez and Song, 2004). To apply this theory, it is convenient to perform the following change of variables. Let S₁ = x₁, S₂ = x₂, E₁ = x₃, E₂ = x₄, I = x₅, J = x₆, and T = x₇. Further, by using the vector notation $X={(x_1,x_2,\dots ,x_7)}^T$, and noting that the force of infection is given by

$\lambda =\beta \frac{x_5+\eta x_6}{x_1+x_2+x_3+x_4+x_5+x_6+x_7},$

the model (3) can be written in the form $\frac{dX}{dt}=F(X)$, with $F={(f_1,f_2,\dots ,f_7)}^T$, as follows

$\begin{array}{cc}\begin{array}{l}\stackrel{.}{x_1}\equiv f_1=\Lambda -{\alpha }_1{x}_1-\lambda x_1+{\theta }_1{x}_2-\mu x_1,\hfill \\ \stackrel{.}{x_2}\equiv f_2={\alpha }_1{x}_1-\sigma \lambda x_2-{\theta }_1{x}_2-\mu x_1,\hfill \\ \stackrel{.}{x_3}\equiv f_3=(1-p_1)\lambda (x_1+\omega ϵx_7)-b_1\lambda x_3-(k_1+{\psi }_1+\mu )x_3+{\theta }_2{x}_4,\hfill \\ \stackrel{.}{x_4}\equiv f_4=(1-p_2)\lambda (\sigma x_2+(1-\omega )ϵx_7)-b_2\lambda x_4-(k_2+{\theta }_2+\mu )x_4+{\psi }_1{x}_3,\hfill \\ \stackrel{.}{x_5}\equiv f_5=p_1\lambda (x_1+\omega ϵx_7)+p_2\lambda (\sigma x_2+(1-\omega )ϵx_7)+b_1\lambda x_3+b_2\lambda x_4\hfill \\ +k_1{x}_3+k_2{x}_4-(\nu \alpha +\pi +d+\mu )x_5,\hfill \\ \stackrel{.}{x_6}\equiv f_6=\pi x_5+\nu {\alpha }_2{x}_5-(r_2+\mu )x_6,\hfill \\ \stackrel{.}{x_7}\equiv f_7=r_2{x}_6-ϵ\lambda x_7-\mu x_7.\hfill \end{array}& \left(\text{A}12\right)\hfill \end{array}$

Consider the case when β = β* is chosen as a bifurcation parameter. Solving for β = β* from ${R}$_T = 1 gives

$\beta ={\beta }^{\ast }=\frac{G_4}{(\mu +\eta \pi +r_2+\eta \nu {\alpha }_2)(G_1+G_2+G_3)}.$

The Jacobian of the transformed system (3), evaluated at the disease-free equilibrium (ξ₁) with β = β*, is given by

$J^{\ast }=J({\xi }_1){|}_{\beta ={\beta }^{\ast }}=\left(\begin{array}{ccccccc}-K_1& {\theta }_1& 0& 0& -\frac{{\beta }^{\ast }{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}& -\frac{\eta {\beta }^{\ast }{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}& 0\\ {\alpha }_1& -K_2& 0& 0& -\frac{\sigma {\beta }^{\ast }{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}& -\frac{\sigma \eta {\beta }^{\ast }{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}& 0,\\ 0& 0& -K_3& {\theta }_2& \frac{(1-p_1){\beta }^{\ast }{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}& \frac{(1-p_1)\eta {\beta }^{\ast }{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}& 0,\\ 0& 0& \psi & -K_4& \frac{(1-p_2)\sigma {\beta }^{\ast }{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}& \frac{(1-p_2)\sigma \eta {\beta }^{\ast }{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}& 0,\\ 0& 0& k_1& k_2& \beta \frac{p_1{x}_1^{\ast }+p_2\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}-K_5& \beta \eta \frac{p_1{x}_1^{\ast }+p_2\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}& 0,\\ 0& 0& 0& 0& \pi +\nu {\alpha }_2& -K_6& 0,\\ 0& 0& 0& 0& 0& r_2& -\mu ,\end{array}\right)$

where K₁ = μ + α₁, K₂ = θ₁ + μ, K₃ = k₁ + ψ + μ, K₄ = k₂ + θ₂ + μ, K₅ = να₂ + μ + d + π, K₆ = r₂ + μ, $x_1^{*}=\frac{\Lambda }{\mu }\frac{\mu +{\theta }_1}{\mu +{\alpha }_1+{\theta }_1}$ and $x_2^{*}=\frac{\Lambda }{\mu }\frac{{\alpha }_1}{\mu +{\alpha }_1+{\theta }_1}$.

Matrix J* has a right eigenvector given by $w={({\omega }_1,{\omega }_2,\dots ,{\omega }_{11})}^T$, where

$\begin{array}{cc}\begin{array}{l}{\omega }_1=-\frac{G_4({\mu }^2+{\theta }_1(2\mu +\sigma {\alpha }_1+{\theta }_1)){\omega }_5}{\mu (G_1+G_2+G_3)(\mu +r){(\mu +{\alpha }_1+{\theta }_1)}^2}\equiv -H_1{w}_5<0\hfill \\ {\omega }_2=-\frac{G_4{\alpha }_1(\mu +\sigma \mu +\sigma {\alpha }_1+{\theta }_1){\omega }_5}{\mu (G_1+G_2+G_3)(\mu +r){(\mu +{\alpha }_1+{\theta }_1)}^2}\equiv -H_2{w}_5<0\hfill \\ {\omega }_3=\frac{G_4[(1-p_1)(\mu +{\theta }_1)(\mu +k_2+{\theta }_2)+\sigma {\alpha }_1{\theta }_2(1-p_2)]{\omega }_5}{(G_1+G_2+G_3)(\mu +r)(\mu +{\alpha }_1+{\theta }_1)[(\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ]}\equiv H_3\omega >0,\hfill \\ {\omega }_4=\frac{G_4[(1-p_1)(\mu +{\theta }_1)\psi +\sigma (1-p_2){\alpha }_1(\mu +k_1+\psi )]{\omega }_5}{(G_1+G_2+G_3)(\mu +r)(\mu +{\alpha }_1+{\theta }_1)[(\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi ]}\equiv H_4\omega >0,\hfill \\ {\omega }_5={\omega }_5>0,\hfill \\ {\omega }_6=\frac{\pi +\nu {\alpha }_2}{r+\nu }{\omega }_5\equiv H_5{\omega }_5>0,\hfill \\ {\omega }_7=\frac{r(\pi +\nu {\alpha }_2)}{r+\nu }{\omega }_5\equiv H_6{\omega }_5>0.\hfill \end{array}& \left(\text{A}13\right)\hfill \end{array}$

Furthermore, J* has a left eigenvector v = (ν₁, ν₂, …., ν₁₁) satisfying v.w = 1, with

$\begin{array}{cc}\begin{array}{l}{\nu }_1=\frac{{\alpha }_1}{({\alpha }_1+\mu )({\theta }_1+\mu )}\equiv Z_1>0,\hfill \\ {\nu }_2=\frac{{\alpha }_1+\mu }{{\theta }_1{\alpha }_1}\equiv Z_2>0,\hfill \\ {\nu }_3=\frac{k_1(\mu +k_2+{\theta }_2)+k_2\psi }{(\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi }{\nu }_5\equiv Z_3{\nu }_5,\hfill \\ {\nu }_4=\frac{k_1{\theta }_2+k_2(\mu +k_1+\psi )}{(\mu +k_1)(\mu +k_2+{\theta }_2)+(\mu +k_2)\psi }{\nu }_5\equiv Z_4{\nu }_5,\hfill \\ {\nu }_5={\nu }_5>0,\hfill \\ {\nu }_6=\frac{1}{r+\mu }[-\frac{\eta {\beta }^{\ast }{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}{\nu }_1-\frac{\sigma \eta {\beta }^{\ast }{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{\nu }_2+\frac{(1-p_1)\eta {\beta }^{\ast }{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}{\nu }_3\hfill \\ +\frac{(1-p_2)\sigma \eta {\beta }^{\ast }{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{\nu }_4+\beta \eta \frac{p_1{x}_1^{\ast }+p_2\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{\nu }_5],\hfill \\ {\nu }_7=0.\hfill \end{array}& \left(\text{A}14\right)\hfill \end{array}$

It follows from Theorem 4.1 in Castillo-Chavez and Song (2004), if we compute the associated non-zero partial derivatives of F(x) (evaluated at the DFE, ξ₁), that the associated bifurcation coefficients, a and b, defined by

$a={\displaystyle \sum _{k,i,j=1}^n{v}_k}{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}(0,0),$

and

$b={\displaystyle \sum _{k,i=1}^n{v}_k}{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial {\beta }^{\ast }}(0,0),$

are computed to be

$\begin{array}{cc}\begin{array}{l}a=-\frac{2{\beta }^{\ast }{\omega }_5^2(1+\eta H_5)}{x_1^{\ast }+x_2^{\ast }}\{(H_1+H_2)\left[\frac{x_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_1+\frac{\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_2\right]\hfill \\ +(H_3+H_4+H_5+H_6+1)\left[\frac{(1-p_1)x_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_3+\frac{\sigma (1-p_2)x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_4+\frac{p_1{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}+\frac{\sigma p_2{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}\right]{\nu }_5\hfill \\ +\left[(1-p_1)H_1{Z}_3+b_1{H}_3{Z}_3+\sigma (1-p_2)H_2{Z}_4+b_2{H}_4{Z}_4+\frac{p_1}{2}{H}_1+\frac{\sigma p_2}{2}{H}_2\right]{\nu }_5\hfill \\ -((H_1+H_2)\left[\frac{(1-p_1)x_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_3+\frac{\sigma (1-p_2)x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_4+\frac{p_1{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}+\frac{\sigma p_2{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}\right]{\nu }_5\hfill \\ +(H_3+H_4+H_5+H_6+1)\left[\frac{x_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_1+\frac{\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}{Z}_2\right]\hfill \\ +H_1{Z}_1+\sigma H_2{Z}_2+[ϵ\omega (1-p_1)H_6{Z}_3+ϵ(1-\omega )(1-p_2)H_6{Z}_4\hfill \\ +\frac{b_1}{2}{H}_3+\frac{b_2}{2}{H}_4+ϵ\omega p_1{H}_6+ϵ(1-\omega )p_2{H}_6]{\nu }_5)\}\hfill \end{array}\hfill & (\text{A}15)\hfill \end{array}$

and

$\begin{array}{cc}b=\left(Z_3\frac{(1-p_1)x_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}+Z_4\frac{(1-p_2)\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}+\frac{p_1{x}_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}+\frac{p_2{x}_2^{\ast }}{x_1^{\ast }+x_2^{\ast }}\right)v_5-Z_1\frac{x_1^{\ast }}{x_1^{\ast }+x_2^{\ast }}+Z_2\frac{\sigma x_2^{\ast }}{x_1^{\ast }+x_2^{\ast }},& \left(\text{A}16\right)\end{array}$

for v₅ > 0.

It follows from Theorem 4.1 in Castillo-Chavez and Song (2004) that the model (3), or the transformed model (A12), will undergo a backward bifurcation if the backward bifurcation coefficients a and b, given by (A15) and (A16), are positive. However, if a is negative, with b still remaining positive, then the system [either (3) or (A12)] will not undergo a backward bifurcation. □