The model consists of five states of susceptible individuals ( S) who are not yet infected, but at a high risk of acquiring TB infection or disease if exposed to an infectious individual(s); newly infected individuals ( L ₁) who have been infected once in their lifetime with Mycobacterium tuberculosis (MTB) and if exposed, they move to either reinfected or active TB state through exogenous reinfection, depending on the host immune response and virulence of the pathogen strain; reinfected individuals ( L ₂) who have been infected multiple times with MTB of either the same or different strains and remain latently infected or move to active TB state through endogenous reactivation; active TB individuals ( I) who are infectious and can transmit TB; treated or recovered individuals ( T) who were previously infectious and were cured by chemotherapy, though some are still on treatment. Note that the difference between newly infected and reinfected individuals is that newly infected individuals are latently infected who have been infected once in lifetime with MTB and they can develop active disease through either exogenous reinfection or endogenous reactivation, while reinfected individuals are latently infected individuals who are multiply infected with MTB of either the same or different virulent strains and they can only develop active disease through endogenous reactivation or remain latently reinfected, depending on the virulence of the pathogen strain and immune response of the host.

The design and development of the model was based on an actual population in Cape Town, South Africa ¹¹ , where many already infected individuals are repeatedly exposed to different infectious sources and become reinfected, and start treatment after acquiring TB disease ¹¹ .

Furthermore, since natural MTB infection and chemotherapy induce an immune response that reduces TB susceptibility ¹² , we take into consideration the fact that treated individuals move to either a newly infected or reinfected state if they maintain immunological memory after successful treatment, and if they become susceptible due to waning of the immune response, they move to a susceptible state. This is because the average duration for the waning of induced immune response by MTB and chemotherapy is not well understood. However, if treatment fails or is not completed, they can move back to active disease and become infectious. As it has been reported that TB disease incidence in TB treated individuals is several-fold higher than in latently infected ¹³ , we explored post-treatment cases returning to either the susceptible, newly infected or reinfected immunological state. This approach makes our model unique and different from the prior studies, which apply infection force for treated individuals moving to latent TB states without considering immunological memory of the host. In order to explore the population impact of effective contact rate and post treatment immune status on TB epidemiology in high TB burden settings, such as Cape Town, the model includes human immunodeficiency virus (HIV) infection, but it does not incorporate: (i) age stratification or (ii) multi-drug resistant TB (MDR-TB).

In the development of the model, we assume that the population consists of fully susceptible individuals, hence, the number of births or recruitment individuals become susceptibles at a rate π y r ^–1. When susceptible individuals become infected with MTB, a proportion p in this group develop active TB very fast and it’s complement proportion 1− p, which constitute the majority become latently infected and move to a newly infected state depending on the virulent strains of the pathogen and host immunological factors. Newly infected individuals are at a high risk of developing active disease through either endogenous reactivation or exogenous reinfection. Here, we take into consideration that one group of newly infected individuals develop active TB slowly through endogenous reactivation at a rate c y r ^–1, and a fraction q of the second group of newly infected individuals develop active TB very fast after infection through exogenous reinfection, while a remaining fraction 1-q of this group move to a reinfected state. Reinfected individuals develop active TB through endogenous reactivation either fast or slowly at a rate y y r ^–1, depending on the immunological state of the host and the virulent strains of the infecting pathogen. Individuals with active TB receive medical treatment and move to a treated state at a rate k, while some die of active TB at a rate µ _t y r ^–1. Treated and successfully cured individuals may become susceptible and move to a susceptible state at a rate r or retain immunological memory and move to either a newly infected at a rate g or a reinfected state at a rate z y r ^–1, without exposure to the force of infection ¹⁴ . This is because latently infected and treated individuals acquire an immune response that reduces TB susceptibility though it may wane and individuals become susceptible. As mentioned earlier, we only implement treatment for active TB cases because in the real world, many TB patients only start treatment after active disease development. However, if TB cases are not well treated, treatment fails or chemotherapy is incomplete, they incubate the infection and relapse back to active TB at a rate α y r ^–1. As stated earlier. it has been noted that latently infected individuals acquire immunity that reduces TB susceptibility ^{5, 15} , though they can be (re)infected and develop active disease if they become frequently exposed to infectious individuals in confined spaces. Thus, ω(0 < ω < 1) is the preventive factor due to prior natural MTB infection, which partially prevents latently infected individuals from TB susceptibility ⁵ . In each of the five states, there is a natural mortality rate µ y r ^–1 as demonstrated in Figure 1.

From the model ( Figure 1), we developed a mathematical model as follows:

S ˙ = π N − ( θ + μ ) S + r T L ˙ 1 = ( 1 − p ) θ S − ( c + ω θ + μ ) L 1 + g T L ˙ 2 = ω ( 1 − q ) θ L 1 − ( y + μ ) L 2 + z T I ˙ = p θ S + ( c + ω q θ ) L 1 + y L 2 − ( k + μ + μ t ) I + α T Τ ˙ = κ Ι − ( r + z + g + α + μ ) T } ( 1 )

where θ = β I N , which is the force of infection, β is the effective contact rate, I N is the prevalence, I is the number of infectious individuals and N is the total number of individuals in the population. Other parameters are described in Table 1.

Disease free equilibrium is defined as the point at which there is no existence of the disease in the population. By computing S, L ₁, L ₂, I and T at disease free equilibrium (i.e., L ₁ = L ₂ = I = T = 0) and assuming that birth rate is equal to death rate, we obtain the disease free equilibrium point as ( S o , L 1 o , L 2 o , I o , T o ) = ( N, 0, 0, 0, 0). This implies that at disease free equilibrium, the total number of individuals in the population is equal to that of susceptible individuals, such that N = S. However, during an epidemic of active disease at any time, t, we have N _{( t)} = S _{( t)}+ L _{1( t)}+ L _{2 ( t)}+ I _{( t)}+ T _{( t)}, implying that the population is comprised of all five states and disease free equilibrium might be unstable. The model is epidemiologically balanced. From the mathematical model ( Equation 1), we compute the basic reproduction number to determine the stability of the disease free equilibrium as discussed below.

Basic reproduction number ( R _o ) is the average number of secondary cases generated by a single case connected with the fully susceptible population ^{16, 20, 21} . We use the next generation method to determine the basic reproduction number, which is an important figure in the epidemiology of infectious diseases. It is among the most frequently estimated quantities to predict the outbreak of infectious diseases, and its value provides information for designing control interventions for the infectious disease burden ^{16, 21} . R _o thus plays an important role in the analysis of infectious disease models, such as TB ¹⁷ . If R _o < 1, the disease free equilibrium is asymptotically stable, and it is unstable if R _o > 1 ^{16, 20, 21} .

The next generation method is used extensively to compute the basic reproduction number, though numerous other methods are discussed in the literature ²¹ . It generates a next generation matrix, and the basic reproduction number is the spectral radius of this matrix. The next generation matrix is defined as a matrix that relates the numbers of newly infected individuals in the various states in consecutive generations ²¹ . The dominant eigenvalue in a next generation matrix is referred to as the basic reproduction number. The next generation matrix ( M) is defined mathematically as

M = F ( S o , L 1 o , L 2 o , I o , T o ) V − 1 ( S o , L 1 o , L 2 o , I o , T o ) ,

and the basic reproduction number is the dominant eigenvalue of M. Where F ( S o , L 1 o , L 2 o , I o , T o ) and V ( S o , L 1 o , L 2 o , I o , T o ) are the transmission and transition matrices ( n × n matrix) at disease free equilibrium, respectively. These are defined mathematically as

F ( S o , L 1 o , L 2 o , I o , T o ) = [ ∂ f i ( x 0 ) x j ] and V ( S o , L 1 o , L 2 o , I o , T o ) = [ ∂ v i ( x 0 ) x j ] for i , j = 1 , ⋯ , n

where x _o is the disease free equilibrium, f _i is the rate of appearance of the new infection in classes, v _i is the transfer of the infection from one class to another and n is the number of infective classes.

However, prior to the formation of the next generation matrix, we generate the transmission ( f) and transition ( v) subsystems from Equation (1). We then linearise subsystems f and v by applying the Jacobian ( J _c ) and obtain matrices F( S, L ₁, L ₂, I, T) and V( S, L ₁, L ₂, I, T), which are n × n transmission and transition matrices, respectively. The transmission subsystem is an infected subsystem that describes the production of newly infected individuals and changes in the states of existing infected individuals. The transition subsystem is the transfer of the infection from one class to another. On the other hand, the transmission matrix describes the production of new infections, and the transition matrix describes changes of infected states, such as the immunity acquisition or removal by death ²¹ . Considering Equation (1), the Jacobian method that is used to linearise subsystems f and v is described mathematically as

J c = ( ∂ S ˙ ∂ S ∂ S ˙ ∂ L 1 ∂ S ˙ ∂ L 2 ∂ S ˙ ∂ I ∂ S ˙ ∂ T ∂ L ˙ 1 ∂ S ∂ L ˙ 1 ∂ L ˙ 1 ∂ L ˙ 1 ∂ L ˙ 2 ∂ L ˙ 1 ∂ I ∂ L ˙ 1 ∂ T ∂ L ˙ 2 ∂ S ∂ L ˙ 2 ∂ L ˙ 1 ∂ L ˙ 2 ∂ L ˙ 2 ∂ L ˙ 2 ∂ I ∂ L ˙ 2 ∂ T ∂ I ˙ ∂ S ∂ I ˙ ∂ L 1 ∂ I ˙ ∂ L 2 ∂ I ˙ ∂ I ∂ I ˙ ∂ T ∂ T ˙ ∂ S ∂ T ˙ ∂ L 1 ∂ T ˙ ∂ L 2 ∂ T ˙ ∂ I ∂ T ˙ ∂ T ) ( 2 )

Hence, from Equation (1), we form the transmission, f and transition, v subsystems as

f = ( 0 ( 1 − p ) β I N S ω ( 1 − q ) β I N L 1 p β I N S + ω q β I N L 1 0 ) and v = ( − π N + ( β I N + μ ) S − r T ( c + ω β I N + μ ) L 1 − g T ( y + μ ) L 2 − z T c L 1 − y L 2 + ( k + μ + μ t ) I − α T − k I + ( r + z + g + α + μ ) T )

In Figure 1, there are three infected states, which are L ₁, L ₂ and I, giving a 3 ×3 matrix. Hence, we form the transmission matrix, F( S, L ₁, L ₂, I, T) by using Jacobian (see Equation (2)) in f for L ₁, L ₂ and I as

F ( S , L 1 , L 2 , I , T ) = ( 0 0 ( 1 − p ) β S N ω ( 1 − q ) β I N 0 ω ( 1 − q ) β L 1 N ω q β I N 0 p β S N + ω q β L 1 N ) ( 3 )

Substituting S, L ₁, L ₂, I and T at disease free equilibrium point ( S o , L 1 o , L 2 o , I o , T o ) = ( N, 0, 0, 0, 0) into Equation (3), we obtain the transmission matrix at disease free equilibrium, F ( S o , L 1 o , L 2 o , I o , T o ) , as

F ( S o , L 1 o , L 2 o , I o , T o ) = ( 0 0 ( 1 − p ) β 0 0 0 0 0 p β ) ( 4 )

By applying the Jacobian with respect to L ₁, L ₂ and I in v above, we obtain the transition matrix, V( S, L ₁, L ₂, I, T), as

V ( S , L 1 , L 2 , I , T ) = ( c + ω β I N + μ 0 ω β L 1 N 0 y + μ 0 − c − y k + μ + μ t ) ( 5 )

Substituting disease free equilibrium point into Equation (5), we obtain the transition matrix at disease free equilibrium, V ( S o , L 1 o , L 2 o , I o , T o ) , as

V ( S o , L 1 o , L 2 o , I o , T o ) = ( c + μ 0 0 0 y + μ 0 − c − y k + μ + μ t ) ( 6 )

Since the transition matrix consists of many variables, to be consistent we used a computational tool (Sage version 8.0; http://www.sagemath.org/) to compute the inverse. Hence the inverse of Equation (6), V ⁻¹ ( S o , L 1 o , L 2 o , I o , T o ) , is demonstrated as

V − 1 ( S o , L 1 o , L 2 o , I o , T o ) = ( 1 c + μ 0 0 0 1 y + μ 0 c ( c + μ ) ( k + μ + μ t ) y ( y + μ ) ( k + μ + μ t ) 1 k + μ + μ t ) ( 7 )

We obtain the next generation matrix, which is the product of F ( S o , L 1 o , L 2 o , I o , T o ) and V ⁻¹ ( S o , L 1 o , L 2 o , I o , T o ) as

F ( S o , L 1 o , L 2 o , I o , T o ) V − 1 ( S o , L 1 o , L 2 o , I o , T o ) = ( ( 1 − p ) β c ( c + μ ) ( k + μ + μ t ) ( 1 − p ) β y ( y + μ ) ( k + μ + μ t ) ( 1 − p ) β k + μ + μ t 0 0 0 p β c ( c + μ ) ( k + μ + μ t ) p β y ( y + μ ) ( k + μ + μ t ) p β k + μ + μ t ) ( 8 )

Thus, the basic reproduction number, which is defined as the dominant eigenvalue in the next generation matrix, was obtained from Equation (8) as

R o = β ( μ p + c ) ( k + μ + μ t ) ( μ + c )

The computed basic reproduction number has a positive real part, implying that it might be greater than 1. We further examined the stability of the disease free equilibrium by computing eigenvalues as discussed below.

Since the basic reproduction number computed in this study has a positive real part, it is likely to be greater than 1, implying that TB transmission is increasing, regardless of the implementation of medical treatment ^{16, 20, 21} . Here, we explore the stability of the disease free equilibrium by computing eigenvalues in Equation (1). Disease free equilibrium is stable if all eigenvalues have negative real parts, and unstable if any of the eigenvalues is positive ²⁰ . Hence, applying the Jacobian in Equation (1), gives

J = ( − ( β I N + μ ) 0 0 − β S N r ( 1 − p ) β I N − ( c + ω β I N + μ ) 0 ( 1 − p ) β S N − ω β L 1 N g 0 ω ( 1 − q ) β I N − ( y + μ ) ω ( 1 − q ) β L 1 N z p β I N c + ω β I N y H α 0 0 0 k − ( r + α + z + g + μ ) ) ( 9 )

where H = p β S N + ω β L 1 N − ( k + μ + μ t )

The characteristic equation of the model at disease free equilibrium is expressed as

                       [ J _{0– λI} ] = 0                      

where J ₀ is the Jacobian matrix (from the model) at disease free equilibrium, λ is the eigenvalue and I is the identity matrix.

Substituting disease free equilibrium point ( S o , L 1 o , L 2 o , I o , T o ) = ( N, 0, 0, 0, 0) into Equation (9), minus the product of eigenvalue and identity matrix diagonally, we obtain the model matrix at disease free equilibrium as follows:

( J 0 − λ I ) = ( − μ − λ 0 0 − β r 0 − c − μ − λ 0 ( 1 − p ) β g 0 0 − y − μ − λ 0 z 0 c y p β − ( k + μ + μ t ) − λ α 0 0 0 k − r − α − z − g − μ − λ ) ( 10 )

From Equation (10), we form the characteristic equation, [ J _{0 −λI} ] = 0, as

                                    (– μ – λ)[(– c – μ – λ)(– y – μ – λ)( pβ – k – μ – μ _t – λ)(– r – α – z – g – μ – λ)]                      

                          – c[(– y – μ – λ)( p β – k – μ – μ _t – λ)(– r – α – z – g – μ – λ)] = 0                      

Simplifying the characteristic equation further, gives

                      (– y – μ – λ)( pβ – k – μ – μ _t – λ)(– r – α – z – g – μ – λ)[ μ( c + μ) + λ(2 μ + c + λ) – c] = 0

From the simplified characteristic equation, we form the following four equations:

                                                                                                                  –( y + μ) – λ = 0                                                                       (11)

                                                                                                                  pβ – ( k + μ + μ _t ) – λ = 0                                                             (12)

                                                                                                           –( r + α + z + g + μ) – λ = 0                                                                (13)

                                                                                                           λ ² + (2 μ + c) λ + μ ² + cμ – c = 0                                                             (14)

Using Equation (11) to Equation (14), we can determine the eigenvalues that predict the stability of the disease free equilibrium. As mentioned earlier, disease free equilibrium is stable if all eigenvalues have negative real parts, and unstable if any of the eigenvalues is positive. Additionally, if one of the eigenvalues is zero then we cannot tell from stability analysis whether the disease free equilibrium is stable or unstable. Thus, considering Equation (11) to Equation (14), the eigenvalues are:

                                                                                                                                  λ ₁ = –( y + μ) < 0                   

                                                                                                                                λ ₂ = p β – ( k + μ + μ _t )                   

                                                                                                                              λ ₃ = –( r + α + z + g + μ) < 0                   

λ 4 , 5 = − ( 2 μ + c ) ± ( 2 μ + c ) 2 − 4 ( μ 2 + c μ − c ) 2

Here, we noted that λ ₁ and λ ₃ have negative real parts, and λ ₂ has a positive real part if pβ > ( k + µ + µ _t ). For stability of the disease free equilibrium, we assume that λ _4,5 have negative real parts (i.e., λ ₄ < 0 and λ ₅ < 0). Disease free equilibrium is stable if pβ < ( k + µ + µ _t ) and unstable if pβ > ( k + µ + µ _t ), where pβ is the effective contact rate for susceptible individuals moving directly to active disease. On the other hand, disease free equilibrium is stable if the effective contact rate for susceptible individuals moving directly to active disease is less than the summation of case detection, natural death and disease death rates. Otherwise disease free equilibrium is unstable.