Impact Factor 4.019

The world's most-cited Microbiology journal

Original Research ARTICLE

Front. Microbiol., 26 January 2016 | https://doi.org/10.3389/fmicb.2015.01530

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

  • Department of Mathematics, University of Benin, Benin City, Nigeria

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, S1, and the “educated,” low risk (high level of awareness) group, S2. The S1 class is “educated” at the per capita rate α1 and thereafter move into the S2 class. Tuberculosis infection can invade the S1 and S2 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 S2 class return to the S1 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 S1 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 α2, and are referred to a TB treatment program under DOTS for effective treatment (in Okuonghae and Omosigho, 2011, α2 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 να2 = 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 να2I will be the total number of identified cases that are tested and treated for tuberculosis under DOTS. Hence, να2 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 r2 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 r1.

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

dS1dt=Λα1S1βS1(I+ηJ)N+θS2μS1,    (1a)
dS2dt=α1S1σβS2(I+ηJ)NθS2μS2,    (1b)
 dEdt=(1p)(βS1(I+ηJ)N+σβS2(I+ηJ)N+ϵβT(I+ηJ)N)        βE(I+ηJ)N(k+μ)E,    (1c)
 dIdt=p(βS1(I+ηJ)N+σβS2(I+ηJ)N+ϵβT(I+ηJ)N)+kE       +βE(I+ηJ)N(να2+μ+d+r1)I,    (1d)
dJdt=να2Ir2JμJ,    (1e)
dTdt=nr2JμTϵβT(I+ηJ)N,    (1f)
dRdt=r1I+mr2JμR.    (1g)

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

Rc=βd+μ+r1+να2(k+pμ)k+μ(μ+r2+ηνα2)μ+r2                      [μ+θμ+α1+θ+σα1μ+α1+θ]    (2)

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., S1 and S2 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
www.frontiersin.org

Figure 1. Schematic diagram of model given in (3), where λ=β(I+η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, E1, and the “educated” latently infected, low risk (high level of awareness) group, E2. It can be stated that due to awareness, individuals in the E2 group are tested in order to know their disease status and take necessary precautionary health measures. The E1 class is made up of individuals from the S1 class and some individuals from the treated class (T) who have a low level of awareness of tuberculosis while the E2 class is made up of individuals from the S2 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 (E1) can become “educated” at the per capita rate ψ and thereafter move into the E2. 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 (E1) while the remaining fraction 1 − ω enter the E2 class.

We also assume that the education programmes produces “temporary immunity” at the per capita rate θ1 (for the susceptible groups, S1 and S2) and at the rate θ2 (for the latent groups, E1 and E2). The cases θ12) = ∞ corresponds to the situation where there is absolutely no immunity (from the education programmes) while θ12) = 0 corresponds to life-long immunity. Hence, θ12) measures the rate at which those in the S2(E2) class return to the S1(E1) 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 < p1(p2) < 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, p1p2 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 b1(b2) are modification parameters accounting for exogenous re-infection of the latently infected individuals in the E1(E2) class, with 0 ≤ b2b1 < 1 and that k1(k2) is the rate of progression of the individuals in the latent state (E1(E2)) to active tuberculosis. Assume that, in addition, to the impact of active cough identification (α2) 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:

dS1dt=Λα1S1βS1(I+ηJ)N+θ1S2μS1,    (3a)
dS2dt=α1S1σβS2(I+ηJ)Nθ1S2μS2,    (3b)
dE1dt=(1p1)(βS1(I+ηJ)N+ωϵβT(I+ηJ)N)          b1βE1(I+ηJ)N(k1+μ+ψ)E1+θ2E2,    (3c)
dE2dt=(1p2)(σβS2(I+ηJ)N+(1ω)ϵβT(I+ηJ)N)          b2βE2(I+ηJ)N(k2+μ+θ2)E2+ψE1,    (3d)
 dIdt=p1β(S1+ωϵT)(I+ηJ)N+p2β(σS2       +(1ω)ϵT)(I+ηJ)N+β(b1E1+b2E2)(I+ηJ)N       +k1E1+k2E2(να2+μ+d+π)I,    (3e)
dJdt=(π+να2)IrJμJ,    (3f)
dTdt=rJμTϵβT(I+ηJ)N,    (3g)

with N = S1 + S2 + E1 + E2 + 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 S1(0) > 0, S2(0) > 0, E1(0) > 0, E2(0) > 0, I(0) > 0, J(0) > 0 and T(0) > 0. Then, the solutions

(S1(t),S2(t),E1(t),E2(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={(S1,S2,E1,E2,I,J,T)+7;NΛμ}

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

ξ1=(S1,S2,E1,E2,I,J,T)   =(Λμμ+θ1μ+α1+θ1,Λμα1μ+α1+θ1,0,0,0,0,0)

We see that the susceptible classes, at the DFE, depend on a factor of the asymptotic population size, Λμ.

The local stability of ξ1 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=(00(1p1)βS1N(1p1)βηS1N00(1p2)βσS2N(1p2)βησS2N00β(p1S1+p2σS2)Nβη(p1S1+p2σS20)N00000)

and

V=((k1+μ+ψ)θ200ψ(k2+μ+θ2)00k1k2(να2+d+μ+π)000(π+να2)(r+μ)).

Thus, the effective reproduction number of model (3), denoted by RT, is given by

RT=β(μ+ηπ+r+ηνα2)[G1+G2+G3]G4,    (4)

where

G1=k1[(μ+σα1+θ1)(k2+θ2)+μ(μ+σα1p2+θ1)],G2=μ(μ+θ2+ψ)[σα1p2+p1(μ+θ1)],G3=k2(μ+θ1)(μp1+ψ)+k2σα1(μ+ψ),G4=(μ+r)(d+μ+π+να2)(μ+α1+θ1)[(μ+k1)            (μ+k2+θ2)+(μ+k2)ψ],

and RT is obtained from ρ(FV−1) with ρ being the spectral radius of the matrix FV−1.

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

Lemma 2.1. The DFE, ξ1, of model (3) is locally asymptotically stable if RT < 1 and unstable if RT > 1.

The threshold quantity, RT, 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 RT < 1) if the initial sizes of the subpopulation of model (3) are in the basin of attraction of the DFE, ξ1, 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

RT=RFμ(πη+μ+r+ηνα2)(πη+μ+ηνα2)(μ+r)RFA1    (5)

where

A1=μ(πη+μ+r+ηνα2)(πη+μ+ηνα2)(μ+r)

and RF = RT|r = 0 is the reproduction number of model (3), under control, without treatment (where RT|r = 0 is the effective reproduction number, with r = 0).

Since the difference of RF from RT is only in treatment, the factor A1 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 A1 will be less than unity (A1 < 1) for all 0 ≤ η, ν < 1.

Generally, if RF < 1, then we cannot expect a TB epidemic and in this case no treatment strategy is needed for control. However, when RF > 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 RF and RT as

ΔT:=RFRT=(1A1)RF.    (6)

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 A1 < 1 in Equation (6). Now, setting RT = 1 and solving for A1, we have the threshold effectiveness of treatment and cough identification taking into account the cost factor and an active case-finding programme:

A1=1RF    (7)

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 A1<A1*. Observe from Equation (7) that A1* is a decreasing function of RF; higher values of A1* results in lower values of RF, a desired outcome.

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

(i)limν1α2A1 =limπν1α2A1=limπν1α2η0A1=μμ+r,    (8)
(ii)limrA1 =limrν1A1=μμ+πη+ηα2,(iii)limrη0A1=1.    (9)

Observe that the limits of A1 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), α2 → ∞ 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 A1, 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 A1 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 π = α2 = ν = 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 π = α2 = ν = 0, RB, is written as

RB=RT|π=α2=ν=0.

Now, we have that RT = A2RB where

A2=d+μd+π+μ+να2πη+μ+r+ηνα2μ+r.

Here, A2 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 RB < 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 RB and RT i.e.,

ΔB:=RBRT=(1A2)RB

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 A2) would mean that ΔB > 0, which is satisfied if A2 < 1.

It is very interesting to note the expressions for the limits of A2, when the limits applied to A1 are used, namely:

(i)limν1α2A2=limπν1α2A2=limν1rA2=η(d+μ)μ+r,    (10)
(ii)limrA2=limrη0A2=d+μd+μ+π+να2and(iii)limπη0ν1α2A2=0.    (11)

Other than the limits of A2 in Equation (11) (iii), the results of the other limits of A2 in Equations (10) and (11) (ii) are adjusted by disease-induced death and are less than unity, when compared to the corresponding limits for A1, as discussed earlier. Interestingly, the limit in Equation (11) (ii) is zero, compared to the same limit that gave μμ+r, for A1; 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 α1 = ψ = θ1 = θ2 = 0. We observe, from model (3) that, for this case, S2 → 0 asymptotically, as t → ∞, so that the latently infected class E2 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 E2 class have life-long awareness with life long benefit (so that θ2 = 0). Hence, the effective reproduction number of the model (3), with α1 = ψ = θ1 = θ2 = 0, RA, is written as

RA=RT|α1=ψ=θ1=θ2=0.

Now, we have that RT = A3RA where

A3=B1B2,

with B1 = (μ(μ + k1)(πη + μ + r + ηνα2)(k1(μ(μ + σα1p2 + θ1) + (μ + σα1 + θ1)(k2 + θ2)) + σk2α1(μ + ψ) + k2(μ + θ1)(μp1 + ψ) + μ(σα1p2 + p1(μ + θ1))(μ + θ2 + ψ))) and B2 = [(k1 + μp1)(μ + r)(πη + μ + ηνα2)(μ + α1 + θ1)[(μ + k1)(μ + k2 + θ2) + (μ + k2)ψ]].

In this case, A3 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 RA < 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 RA and RT i.e.,

ΔA:=RART=(1A3)RB,

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

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

(i)limrα1θ10πA3=limπrψθ20A3=limν1rπθ20θ10ψα1α2A3=0,    (12)
(ii)limα1ψθ10θ20A3=μσ(μ+k1)(k2+p2μ)(πη+μ+r+ηνα2(μ+k2)(k1+p1μ)(μ+r)(πη+μ+ηνα2),    (13)

and

(iii)limν1πθ20θ10ψα1α2A3=μσ(μ+k1)(k2+p2μ)(μ+k2)(k1+p1μ)(μ+r).    (14)

Clearly, the limits of A3 that includes r (treatment rate) are zero, and this is a positive outcome. However, for the limits of A3 in Equations (13) and (14), the limiting values of A3 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 A3 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 = ν = α2 = π = α1 = ψ = θ1 = θ2 = 0. This implies that, from model (3), (S2, E2, J, T) → 0 asymptotically, as t → ∞. Hence the basic reproduction number is obtained as

R0=RT|r=ν=α2=π=α1=ψ=θ1=θ2=0.

So that RT = A4R0 where

A4=D1D2,

with D1 = ((d + μ)(μ + k1)(πη + μ + r + ηνα2)(k1(μ(μ + σα1p2 + θ1) + (μ + σα1 + θ1)(k2 + θ2)) + σk2α1(μ + ψ) + k2(μ + θ1)(μp1 + ψ) + μ(σα1p2 + p1(μ + θ1))(μ + θ2 + ψ))) and D2 = [(k1 + μp1)(μ + r)(d + π + μ + να2)(μ + α1 + θ1)[(μ + k1)(μ + k2 + θ2) + (μ + k2)ψ]].

If R0 < 1, then TB cannot develop into an epidemic and no control strategy is needed for TB control. If we take the difference between R0 and RT i.e.,

Δ0:=R0RT=(1A4)R0

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 > 0, which is satisfied if A4 < 1.

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

(i)limα1ψθ10θ20A4=(d+μ)σ(μ+k1)(k2+p2μ)(πη+μ+r+ηνα2)(μ+k2)(k1+p1μ)(μ+r)(d+π+μ+να2)    (15)

and

(ii)limν1πθ20θ10ψα1α2A4=η(d+μ)σ(μ+k1)(k2+p2μ)(μ+k2)(k1+p1μ)(μ+r)    (16)

The difference is that the limits of A4 are disease-induced death adjusted (compared to the limits of A3) and, in Equations (15) and (16), the limiting results for A4 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, A1 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, A4 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 A1 < 1, A2 < 1, A3 < 1 or A4 < 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 > 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 A1, A2, A3 and A4. However, our interests are in investigating control strategies involving the combination of the following parameters: ν, π, ψ, θ1, θ2, α2, r, α1 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, RT

Using the threshold parameter, RT, 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

limα1ψθ10θ20RT=βσ(k2+p2μ)(πη+μ+r+ηνα2)(k2+μ)(μ+r)(d+π+μ+να2)>0,    (17)
limα1ψθ10θ20rRT=βσ(k2+p2μ)(k2+μ)(d+π+μ+να2)>0,    (18)
limν1α2α1ψθ10θ20rRT=limα2ν1πrRT=0.    (19)

It seems, from the limits of RT 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 RT, 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 RT as a function of the chronic cough identification rate (α2) and the associated cost factor (ν) and as a function of the awareness rate for the susceptible group (α1) and the chronic cough identification rate (α2), respectively. Using the base parameter values in Table 1, we observe from Figures 24 that the indicated parameters have a positive effect on the effective reproduction number, RT. 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 RT, albeit at different rates.

FIGURE 2
www.frontiersin.org

Figure 2. Contour plot of RT as a function of π and r.

FIGURE 3
www.frontiersin.org

Figure 3. Contour plot of RT as a function of ν and α2.

FIGURE 4
www.frontiersin.org

Figure 4. Contour plot of RT as a function of α1 and α2.

TABLE 1
www.frontiersin.org

Table 1. Parameter information.

From the expressions in the limits of RT 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, α1 → ∞(ψ → ∞), θ1 → 0(θ2 → 0), α2 → ∞, 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, α2 → ∞, 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 RT 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 RT with respect to the key parameters (r2, α2, α1, θ1, θ2, ψ, π 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 RTr<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 RTα1>0 (RTθ1<0) if

σ>σ=k1k2+k1θ2+k2ψ+μ[k1+k2p1+p1(μ+θ2+ψ)]k1k2+k1θ2+k2ψ+μ[k1p2+k2+p2(μ+θ2+ψ)].    (20)

However, RTα1<0 (RTθ1>0) if σ < σ* and RTα1=0 (RTθ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 RT, 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., p1 and p2. Clearly, if p1 = p2 = 1, then σ* = 1 (which is also the case when p1 = p2 = p and k1 = k2 = 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, p1 and p2 and the rate of endogenous reactivation, k1 and k2 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 RTθ2<0 (RTψ>0) if

k1<k2.    (21)

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), k1, 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., k1 > k2; there will be no impact on TB dynamics if k1 = k2 and a detrimental impact when k1 < k2. 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 (θ2). 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 k1 > k2, no impact if k1 = k2 and a detrimental impact if k1 < k2. The result is reversed for the rate of loss of awareness by such latently infected: a positive impact if k1 < k2, no impact when k1 = k2 and a detrimental impact if k1 > k2.

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 p1 and p2, and probably k1 and k2) 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, k1 and k2.

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 (p1, p2, k1 and k2) 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 (p1, p2, k1, and k2) 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 RTπ<0, RTα2<0, and RTν<0 if

d<d=(1η)μ+rη,    η0    (22)

or

η<η=μ+rμ+d    (23)

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 rd. 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 RT, 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 ν, α1, α2 and others like θ1, π, ψ and θ2, 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 (RT) 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)Λμ as t → ∞. Let β^=μβΛ so that the force of infection now becomes

λ=β^(I+ηJ).    (24)

Also, let 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 ξ2=(S1**,S2**,E1**,E2**,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:

f(λe)=A1λe4+A2λe3+A3λe2+A4λe+A5=0    (25)

where the coefficients, A1, …, A5 and the endemic equilibrium, ξ2, 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 A1, of the quartic (25), is always positive, and A5 is positive (negative) if 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 b1 = 0 and R^T<1 (absence of exogenous re-infection from the E1 class).

ii has four or two endemic equilibria if A2 < 0, A3 > 0, A4 < 0 and R^T<1.

iii has two endemic equilibria if A3 is of the same sign as A2 or A4 and R^T<1.

iv has two endemic equilibria if A2 > 0, A3 < 0, A4 > 0 and R^T<1.

v no endemic equilibrium otherwise when R^T<1.

vi a unique endemic equilibrium when 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., b1 and b2, 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, b1 = 0, the system (3) will not have a backward bifurcation when R^T<1. Clearly, setting b1 = 0 in Equation (25), yields the quadratic equation

f(λe)=A^3λe2+A^4λe+A5=0    (26)

where Â3 and Â4 are now evaluated from A3 and A4, respectively, with b1 = 0. It is easy to see that Â3 > 0 and A5 > 0 (the latter when R^T<1). We want to show that, since σ ≤ 1, the coefficient Â4 > 0, so that the quadratic (26) will not have any positive roots, which implies that there is no endemic equilibrium when b1 = 0 in Equation (25). With b1 = 0 in Equation (25), we see that

A^4=Λ[(μ+π+να2)(μ+μσ+σα1+θ1)((μ+k1)           (μ+k2+θ2)+(μ+k2)ψ)]Λ[β(μ+η(π+να2))           σ((k1+μp1)(μ+k2+θ2)+(k2+μp1)ψ)]

Since σ ≤ 1, it then implies that,

A^4Λ[(μ+π+να2)(μ+μσ+σα1+θ1)((μ+k1)           (μ+k2+θ2)+(μ+k2)ψ1)]Λ[β(μ+η(π+να2))           ((k1+μp1)(μ+k2+θ2)+(k2+μp1)ψ1)].    (27)

From R^T<1, it follows that

β^G4(μ+ηπ+ηνα2)(G1+G2+G3),

where β^ is now evaluated with r = d = 0. Substituting the expression for β^ into Equation (27), we see that

A^4Λ[(μ+π+να2)(μ+μσ+σα1+θ1)((μ+k1)           (μ+k2+θ2)+(μ+k2)ψ1)]          Λ[(G4(μ+ηπ+r2+ηνα2)(G1+G2+G3))(μ+η          (π+να2))((k1+μp1)(μ+k2+θ2)+(k2+μp1)ψ1)].    (28)

Simplifying this further, we have that

A^4Λ[(μ+π+να2)[((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)            (k1(μ3+k2(μ+α1+θ1)2+μ(2μ+α1+θ1)(p2α1+θ1)        +(μ+α1+θ1)2θ2)+μα1p2(2μ+α1+θ1)(μ+θ2+ψ)        +μp1(μ2+θ1(2μ+α1+θ1))(μ+k2+θ2+ψ)        +k2(μα1(2μ+α1+θ1)+(μ+α1+θ1)2ψ))]>0,

for σ ≤ 1. Therefore, the polynomial (26) will not have a positive root, hence no endemic equilibrium, when b1 = 0 and R^T<1, ruling out the possibility of a backward bifurcation, for this case.

However, when b2 = 0 and b1 ≠ 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 E1 class i.e., latently infected individuals with low level of awareness), the quartic (25) now reduces to a cubic equation

f(λe)=A˜2λe3+A˜3λe2+A˜4λe+A5=0    (29)

where A2~,A3~ and A4~ are obtained from evaluating the corresponding coefficients of (25) with b2 = 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 R^T < 1, suggesting the presence of a backward bifurcation for the case b2 = 0 and b1 ≠ 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, b1 and b2, are present in the system i.e., b1 ≠ 0 and b2 ≠ 0.

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

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

i both exogenous re-infection parameters, b1 and b2, are absent in the system i.e., b1 = 0 and b2 = 0.

ii the exogenous re-infection parameter b1 is absent in the system (even when b2 ≠ 0) i.e., b1 = 0 and b2 ≠ 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 b1 = 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 b2 (the exogenous re-infection parameter for the latently infected individuals with high awareness level).

FIGURE 5
www.frontiersin.org

Figure 5. Backward bifurcation when b1 ≠ 0 and b2 ≠ 0.

FIGURE 6
www.frontiersin.org

Figure 6. Backward bifurcation when b1 ≠ 0 and b2 = 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 RT < 1. We claim the following:

Theorem 2.4. The model (3) exhibits backward bifurcation at RT = 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 RT = 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, b1 and b2, there could be other parameters that could determine whether the system (3) can undergo a backward bifurcation at RT = 1. These parameters, we believe, are liked to the awareness parameters (i.e., α1, ψ, θ1, and θ2) and the fractions of fast progressions, p1 and p2. 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
www.frontiersin.org

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

FIGURE 8
www.frontiersin.org

Figure 8. Simulations of model (3) showing the number of infected individuals. Here, the awareness rate for the susceptible individuals (α1) is varied from 0 to 30 with θ1 = 20, θ2 = 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
www.frontiersin.org

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

FIGURE 10
www.frontiersin.org

Figure 10. Simulations of model (3) showing the number of infected individuals. Here, the rate at which susceptibles lose awareness (θ1) is varied from 0 to 30 with α1 = 30 and α2 = 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
www.frontiersin.org

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 α1 = α2 = 5, θ1 = θ2 = 30 and ψ = 5.

FIGURE 12
www.frontiersin.org

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 α1 = α2 = 5, θ1 = θ2 = 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
www.frontiersin.org

Figure 13. Simulations of model (3) showing the number of infected individuals. Here, the cost factor (ν) is varied from 0 to 1 with α1 = α2 = 5, θ1 = θ2 = 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 = θ2 = 1).

FIGURE 14
www.frontiersin.org

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 α1 = α2 = 5, θ1 = θ2 = 1 and ψ = 30.

FIGURE 15
www.frontiersin.org

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 α1 = α2 = 20, θ1 = θ2 = 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 RT. 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, ν, α1, α2, θ1, θ2, ψ, 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 R0 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 t1 = sup{t > 0 : S1(t) > 0, S2(t) > 0, E1(t) > 0, E2(t) > 0, I(t) > 0, J(t) > 0, T(t) > 0} > 0. It follows from (3a) that

dS1dt=Λα1S1λS1μS1+θ1S2Λα1S1λS1μS1,

which can be written as

ddt{S1(t)exp[(α1+μ)t+0tλ(τ)dτ] }Λ{exp[(α1+μ)t+0tλ(τ)dτ] }.

Thus,

S1(t1)exp[(α1+μ)t1+0t1λ(τ)dτ]S1(0)0t1Λ{exp[(α1+μ)y+0yλ(τ)dτ] }dy.

so that,

S1(t1)S1(0)exp[(α1+μ)t10t1λ(τ)dτ]+         {exp[(α1+μ)t10t1λ(τ)dτ] }         0t1Λ{exp[(α1+μ)y+0yλ(τ)dτ] }dy>0.

Similarly, it can be shown that S2(t) > 0, E1(t) > 0, E2(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

dNdt=ΛμNdI.    (A1)

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

N(t)N(0)eμt+Λμ(1eμt).    (A2)

In particular, if N(0)Λμ, then N(t)Λμ. Thus, D is a positively invariant set. Furthermore, if N(0)Λμ, then either the solution enters D in finite time or N(t) approaches to Λμ 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 RT

RTr=βη(π+να2)P1((μ+r)2(d+π+μ+να2)(μ+α1+θ1)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)),    (A3)
RTα1=β(πη+μ+r+ηνα2)(μ+θ1)P2((μ+r)(d+π+μ+να2)(μ+α1+θ1)2    ((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)),    (A4)
RTθ1=βα1(πη+μ+r+ηνα2)P2((μ+r)(d+π+μ+να2)(μ+α1+θ1)2    ((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)),    (A5)
RTθ2=βμ(πη+μ+r+ηνα2)(k1k2)[(1p1)(μ+θ1)ψ+σ(1p2)α1(μ+k1+ψ)]((μ+r)(d+π+μ+να2)(μ+α1+θ1)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)2),    (A6)
RTψ=βμ(πη+μ+r+ηνα2)(k1k2)[(1p1)(μ+θ1)(μ+k2+θ2)+σ(1p2)α1θ1]((μ+r)(d+π+μ+να2)(μ+α1+θ1)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)2),    (A7)
RTπ=β(ηd(1η)μr)P1((μ+r)(d+π+μ+να2)2(μ+α1+θ1)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)),    (A8)
RTα2=βν(ηd(1η)μr)P1((μ+r)(d+π+μ+να2)2(μ+α1+θ1)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)),    (A9)
RTν=βα2(ηd(1η)μr)P1((μ+r)(d+π+μ+να2)2(μ+α1+θ1)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)),    (A10)

where

P1 = k1(μ(μ + σα1p2 + θ1 + (μ + σα1 + θ1)(k2 + θ2)) + σα1k2(μ + ψ) + k2(μ + θ1)(μp1 + ψ) + μ(σα1p2 + p1(μ + θ1))(μ + θ2 + ψ)),

and

P2 = −[(1 − σ)(k1k2 + k1θ2 + k2ψ) + k1μ(1 − σp2)] − (k2μ(p1 − σ) + μ(p1 − σp2)(μ + θ2 + ψ).

Appendix C

Existence of EEP for Special Case

The coefficients of the polynomial (25) are given as

A1=z2ϑσΛ(μ+π+να2),A2=Λz(μ+π+να2)[μ(σ+ϑ(z+σ+zσ))+ϑσk1+σk2+zϑσα1       +zϑθ1+σ(θ2+ϑψ)][β(μ+η(π+να2))z2ϑΛσ]A3=Λ(μ+π+να2)[k1(μ(σ+zϑ(1+σ))+σk2+zϑ(σα1+θ1)+σθ2)       +zθ1((1+ϑ+zϑ)μ+θ2+ϑψ)+k2(μ(z+σ+zσ)+zσα1+zθ1+σψ)       +zα1(μ(σ+ϑ(z+σ))+σθ2+ϑσψ)+μ[μ(σ+z(1+σ+ϑ(1+z+σ)))       +(z+σ+zσ)θ2(σ+zϑ(1+σ))ψ]][β(μ+η(π+να2))        zΛ(zϑμ+μσ+ϑσk1+σk2+ϑμσp1+zϑσα1+zϑθ1+σθ2+ϑσψ)]A4=Λ(μ+π+να2)[k1[k2(μ+μσ+σα1+θ1)+θ1(μ+zϑμ+θ2)+α1(μ(zϑ+σ)+σθ2)       +μ(μ(1+zϑ+σ)+(1+σ)θ2)]+k2[θ1(μ+zμ+ψ1)+α1(μ(z+σ)+σψ)       +μ(μ(1+z+σ)+(1+σ)ψ)]+μ[θ1((1+z+zϑ)μ+(1+z)θ2+(1+zϑ)ψ)       +α1(μ(z+zϑ+σ)+(z+σ)θ2+(zϑ+σ)ψ)+μ(μ(1+z+zϑ+σ)       +(1+z+σ)θ2+(1+zϑ+σ)ψ)]       β(μ+η(π+να2))Λ[k1(μ(zϑ+σ)+σk2+zϑ(σα1+θ1)+σθ2)       +k2(μσp1+z(μ+σα1+θ1)+σψ)+μp1(μ(zϑ+σ)+zϑθ1+σ(θ2+ψ))        z((μ+θ1)(μ+θ2+ϑψ)+σα1(μnp2+θ2+ϑ(μ+ψ)))]A5=Λμ(μ+α1+θ1)(μ+π+να2)((μ+k1)(μ+k2+θ2)+(μ+k2)ψ)[1R^T]

where z = b1 and ϑ = b1b2 with R^T=RT|r=d=0 being the reproduction number of the model (3) where r = d = 0. The endemic steady state solutions are

  S1∗∗=Λ(μ+λeσ+θ1)α1(μ+λeσ)+(λe+μ)(μ+λeσ+θ1),  S2∗∗=Λα1α1(μ+λeσ)+(λe+μ)(μ+λeσ+θ1),  E1∗∗=1b1λe+(k1+μ+ψ)[(1p1)S1∗∗λe+θ2E2∗∗],E2∗∗=1b2λe+(k2+μ+θ2)[(1p2)σS2∗∗λe+ψE1∗∗],  I∗∗=1να2+μ+π[p1S1∗∗λe+p2σS2∗∗λe+(b1E1∗∗+b2E2∗∗)+k1E1∗∗+k2E2∗∗], J∗∗=π+να2μI∗∗    (A11)

Appendix D

Proof of Theorem 2.4

Suppose

ξ3=(S1∗∗,S2∗∗,E1∗∗,E2∗∗,I∗∗,J∗∗,T∗∗)

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 S1 = x1, S2 = x2, E1 = x3, E2 = x4, I = x5, J = x6, and T = x7. Further, by using the vector notation X=(x1,x2,,x7)T, and noting that the force of infection is given by

λ=βx5+ηx6x1+x2+x3+x4+x5+x6+x7,

the model (3) can be written in the form dXdt=F(X), with F=(f1,f2,,f7)T, as follows

x1.f1=Λα1x1λx1+θ1x2μx1,x2.f2=α1x1σλx2θ1x2μx1,x3.f3=(1p1)λ(x1+ωϵx7)b1λx3(k1+ψ1+μ)x3+θ2x4,x4.f4=(1p2)λ(σx2+(1ω)ϵx7)b2λx4(k2+θ2+μ)x4+ψ1x3,x5.f5=p1λ(x1+ωϵx7)+p2λ(σx2+(1ω)ϵx7)+b1λx3+b2λx4              +k1x3+k2x4(να+π+d+μ)x5,x6.f6=πx5+να2x5(r2+μ)x6,x7.f7=r2x6ϵλx7μx7.    (A12)

Consider the case when β = β* is chosen as a bifurcation parameter. Solving for β = β* from RT = 1 gives

β=β=G4(μ+ηπ+r2+ηνα2)(G1+G2+G3).

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

J=J(ξ1)|β=β=(K1θ100βx1x1+x2ηβx1x1+x20α1K200σβx2x1+x2σηβx2x1+x20,00K3θ2(1p1)βx1x1+x2(1p1)ηβx1x1+x20,00ψK4(1p2)σβx2x1+x2(1p2)σηβx2x1+x20,00k1k2βp1x1+p2σx2x1+x2K5βηp1x1+p2σx2x1+x20,0000π+να2K60,00000r2μ,)

where K1 = μ + α1, K2 = θ1 + μ, K3 = k1 + ψ + μ, K4 = k2 + θ2 + μ, K5 = να2 + μ + d + π, K6 = r2 + μ, x1*=Λμμ+θ1μ+α1+θ1 and x2*=Λμα1μ+α1+θ1.

Matrix J* has a right eigenvector given by w=(ω1,ω2,,ω11)T, where

ω1=G4(μ2+θ1(2μ+σα1+θ1))ω5μ(G1+G2+G3)(μ+r)(μ+α1+θ1)2H1w5<0ω2=G4α1(μ+σμ+σα1+θ1)ω5μ(G1+G2+G3)(μ+r)(μ+α1+θ1)2H2w5<0ω3=G4[(1p1)(μ+θ1)(μ+k2+θ2)+σα1θ2(1p2)]ω5(G1+G2+G3)(μ+r)(μ+α1+θ1)[(μ+k1)(μ+k2+θ2)+(μ+k2)ψ]H3ω>0,ω4=G4[(1p1)(μ+θ1)ψ+σ(1p2)α1(μ+k1+ψ)]ω5(G1+G2+G3)(μ+r)(μ+α1+θ1)[(μ+k1)(μ+k2+θ2)+(μ+k2)ψ]H4ω>0,ω5=ω5>0,ω6=π+να2r+νω5H5ω5>0,ω7=r(π+να2)r+νω5H6ω5>0.    (A13)

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

ν1=α1(α1+μ)(θ1+μ)Z1>0,ν2=α1+μθ1α1Z2>0,ν3=k1(μ+k2+θ2)+k2ψ(μ+k1)(μ+k2+θ2)+(μ+k2)ψν5Z3ν5,ν4=k1θ2+k2(μ+k1+ψ)(μ+k1)(μ+k2+θ2)+(μ+k2)ψν5Z4ν5,ν5=ν5>0,ν6=1r+μ[ηβx1x1+x2ν1σηβx2x1+x2ν2+(1p1)ηβx1x1+x2ν3      +(1p2)σηβx2x1+x2ν4+βηp1x1+p2σx2x1+x2ν5],ν7=0.    (A14)

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, ξ1), that the associated bifurcation coefficients, a and b, defined by

a=k,i,j=1nvkwiwj2fkxixj(0,0),

and

b=k,i=1nvkwi2fkxiβ(0,0),

are computed to be

a=2βω52(1+ηH5)x1+x2{(H1+H2)[x1x1+x2Z1+σx2x1+x2Z2]   +(H3+H4+H5+H6+1)[(1p1)x1x1+x2Z3+σ(1p2)x2x1+x2Z4+p1x1x1+x2+σp2x2x1+x2]ν5   +[(1p1)H1Z3+b1H3Z3+σ(1p2)H2Z4+b2H4Z4+p12H1+σp22H2]ν5       ((H1+H2)[(1p1)x1x1+x2Z3+σ(1p2)x2x1+x2Z4+p1x1x1+x2+σp2x2x1+x2]ν5       +(H3+H4+H5+H6+1)[x1x1+x2Z1+σx2x1+x2Z2]       +H1Z1+σH2Z2+[ϵω(1p1)H6Z3+ϵ(1ω)(1p2)H6Z4       +b12H3+b22H4+ϵωp1H6+ϵ(1ω)p2H6]ν5)}    (A15)

and

b=(Z3(1p1)x1x1+x2+Z4(1p2)σx2x1+x2+p1x1x1+x2+p2x2x1+x2) v5Z1x1x1+x2+Z2σx2x1+x2,    (A16)

for v5 > 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.

Keywords: mathematical model, tuberculosis, awareness level, bifurcation, simulation

Citation: Okuonghae D and Ikhimwin BO (2016) Dynamics of a Mathematical Model for Tuberculosis with Variability in Susceptibility and Disease Progressions Due to Difference in Awareness Level. Front. Microbiol. 6:1530. doi: 10.3389/fmicb.2015.01530

Received: 19 September 2015; Accepted: 21 December 2015;
Published: 26 January 2016.

Edited by:

Cristina Vilaplana, Fundació Institut d’Investigació en Ciències de la Salut Germans Trias i Pujol, Spain

Reviewed by:

Paras Jain, Albert Einstein College of Medicine, USA
Clara Prats, Universitat Politècnica de Catalunya, Spain

Copyright © 2016 Okuonghae and Ikhimwin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Daniel Okuonghae, danny.okuonghae@corpus-christi.oxon.org; daniel.okuonghae@uniben.edu