Dynamics of tuberculosis transmission with exogenous reinfections and endogenous reactivation

Accepted Manuscript Dynamics of tuberculosis transmission with exogenous reinfections and endogenous reactivation Subhas Khajanchi, Dhiraj Kumar Das, Tapan Kumar Kar

PII: DOI: Reference:

S0378-4371(18)30014-1 https://doi.org/10.1016/j.physa.2018.01.014 PHYSA 19094

To appear in:

Physica A

Received date : 16 October 2017 Revised date : 2 January 2018 Please cite this article as: S. Khajanchi, D.K. Das, T.K. Kar, Dynamics of tuberculosis transmission with exogenous reinfections and endogenous reactivation, Physica A (2018), https://doi.org/10.1016/j.physa.2018.01.014 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Highlights

Highlights:  A TB transmission model with exogenous reinfection and endogenous reactivation.  Backward bifurcation arises as the exogenous reinfection crosses a critical value.  As the transmission rate exceeds a critical value a bifurcating limit cycle occurs.  Stability of the bifurcating periodic solution is studied.  Simulations covering the breadth of the feasible parameter space are conducted.

*Manuscript Click here to view linked References

Dynamics of tuberculosis transmission with exogenous reinfections and endogenous reactivation Subhas Khajanchia , Dhiraj Kumar Dasb , Tapan Kumar Karb a

b

Department of Mathematics, Bankura University, Bankura - 722155, India Department of Mathematics, Indian Institute of Engineering Science and Technology Shibpur, Howrah 711103, India

Abstract We propose and analyze a mathematical model for tuberculosis (TB) transmission to study the role of exogenous reinfection and endogenous reactivation. The model exhibits two equilibria: a disease free and an endemic equilibria. We observe that the TB model exhibits transcritical bifurcation when basic reproduction number R0 = 1. Our results demonstrate that the disease transmission rate β and exogenous reinfection rate α plays an important role to change the qualitative dynamics of TB. The disease transmission rate β give rises to the possibility of backward bifurcation for R0 < 1, and hence the existence of multiple endemic equilibria one of which is stable and another one is unstable. Our analysis suggests that R0 < 1 may not be sufficient to completely eliminate the disease. We also investigate that our TB transmission model undergoes Hopf-bifurcation with respect to the contact rate β and the exogenous reinfection rate α. We conducted some numerical simulations to support our analytical findings. Keywords: Tuberculosis, Basic reproduction number, Backward bifurcation, Transcritical bifurcations, Hopf-bfurcation 1. Introduction Tuberculosis (TB) is an airborne-transmitted bacterial disease caused by infection with bacterium Mycobacterium tuberculosis (Mtb) with almost one-third of the world human population as its reservoir [1]. Mtb customarily attacks the lungs, but can also affect other organ of the human body. The progression of TB is relatively slow dynamics compared to other viral infections: Mtb doubles approximately every 18 - 48 hours, while other pathogens have doubling times on the order of minutes. The primary TB progression may have a time course 1 - 5 years with reactivation an even longer process, at least 33 years [2, 6]. Thus, the epidemics of TB must be studied and estimated over acutely long windows in time. Individuals who are latently infected with TB transmission do not have enough symptoms or not capable of transmitting TB [3]. In spite of that, what happens in the primary stages of the interaction between host-pathogen is likely relevant in understanding different infection outcomes. Recent research has explored the pivotal implications of individual behavior and heterogeneous contact patterns in networked populations, as well as the many feedback loops that exist between vaccinating behavior and disease Preprint submitted to Physica A

January 2, 2018

propagation [29, 30, 31]. In addition to these different vaccination strategies, network structure also plays an important role in the transmission routes and infection levels [30]. The major difference between TB transmission and other communicable diseases is that at the primary infection, nearly 5 - 10% of latently infected individuals develop the progressive disease (active TB), that is, approximately 90 - 95% remain latently infected. Meanwhile, exposed individuals may persist in this latent phase for long period of time [17]. The age of infection and chronological age are crucial factors in disease development [5]. The development towards progressive disease (active TB) may influence with re-exposure to TB bacilli through repeated contacts with individuals with progressive disease (active TB). Thus, we must not only look at TB transmission as the development from primary infection but also incorporate the possibility of exogenous reinfection or endogenous reactivation. Exogenous reinfection plays a critical role in the development of active TB. Nevertheless, the risk factor of developing TB transmission as a upshot of exogenous reinfection is lower than that of developing the primary case for most age groups [4]. Exogenous reinfection has been studied in immunocomponent [7, 17, 25, 26, 27] and immunosuppressed individuals [8]. The significance of exogenous reinfection when compared with the effect of primary infection and endogenous reactivation was illustrated using epidemiological modeling [9]. The importance of exogenous reinfection in the progression of active TB from latent period for populations who face high contact rates needs to be appraised. In our study, we introduced exogenous reinfections into our typical epidemiological model for the dynamics of TB transmission. We investigate that the exogenous reinfection plays an elementary role in the dynamics of TB transmission. Indeed, our mathematical analysis demonstrates that reinfection may - theoretically - increase the number of active TB, may allow less predictable dynamics, and thus may decrease the effectiveness of public health measures. Uys et al. [26] studied a mathematical model for TB transmission with exogenous reinfection. In their study, they assumed that the reinfection rate is a multiple of the rate of first-time infection. Feng et al. [17] proposed a TB model with exogenous reinfection and showed that an individual can be infected by one infectious individual per contact per unit time. The word exogenous reinfection is known as super-infection was first introduced by Martcheva and Thieme [28]. Using the same term super-infection, the authors Yang and Raimundo [32] investigate the role of multiple infections and long latency periods in the dynamics of TB transmission. They investigate that for short latent period, that is, below a threshold value the backward bifurcation disappears whereas the rates of super-infection and reinfection decreased their system undergoes forward bifurcation. The basic goal for studying the infectious diseases is to gain understanding the process, and to design better treatment policies, or improve existing policies to eradicate disease or at least to improve patient’s quality of life. Mathematical modeling is a viable tool in understanding the communicable diseases such as TB, ebola, rabies, influenza, West Nile virus, Zika, measles, pertussis and HIV/AIDS etc. [1, 2, 17]. These communicable diseases cause repeated epidemic outbreaks, and their associated transmission rates solely depends on agespecific contact rates. The etiological agents of these infectious diseases are different viruses, 2

but all are efficient of producing similar symptoms at the individual level. Common responses comprises short latent periods, followed by comparatively short infectious periods, and permanent immunity after recovery. In contrast, the investigation of the control and spread of TB transmission using mathematical and statistical models has not gained sufficient attention. In spite of that, a series of mathematical models have already been developed, and each one contributes in their own way to a better understanding of TB transmission and the qualitative dynamics which determine the patient’s outcome [1, 10, 15, 22, 23, 24]. The dynamics of TB transmission shows long and changeable periods of latency. In [10], the authors showed that the incorporation of a distributed delay to model the long and different periods of latency do not change the qualitative behavior of their TB model. Still, the authors studied a TB model with reinfection as long as variable periods of latency [17]. In their article [17], the authors conclude that the results will hold in TB model with realistic distributed delay. Bowong and Tewa [22] introduced a SEI type of TB model with a general contact rate and derived the global stability of equilibria by using suitable Lyapunov stability theory and LaSalles invariant principle. Huo and Zou [24] constructed a TB model with two different kinds of treatment, namely, treatment at home and treatment in hospital. They prove the global stability for R0 < 1 and R0 > 1 in case of disease free equilibrium and endemic equilibrium point respectively, by constructing suitable Lyapunov function. Optimal treatment strategy is an important tool to control the disease transmission [10, 19, 20, 21]. In [19], the authors used the concept of optimal control theory to eradicate TB transmission. In their study they showed that the number of infectious individuals decreased very rapidly in presence of treatment control whereas the number of treated individuals increase throughout the treatment period. The qualitative dynamics of most epidemiological models depends on the threshold value known as basic reproduction number R0 [33, 34, 36]. In general, if R0 < 1, the disease will die out in time whereas if R0 > 1, the disease will persist and thus transmit through the individuals. Such type of characterization is familiar in most of the epidemiological models. Recently, many researchers investigate backward bifurcations [17, 18, 19, 23, 28, 37, 38] due to nonlinear incidence rates, macro parasite infection and age-structure and so on. Backward bifurcations can be characterized by multiple endemic equilibria, may give rise to a disease persist even though R0 < 1. Specifically, one can observe a stable equilibrium point and one unstable equilibrium point for R0 < 1. The appearance of backward bifurcation is important in a practical sense as control programs must reduce R0 further than below one to eradicate a disease. Singer and Kirschner [37] investigate a TB epidemic model with exogenous reinfection to study the influence of backward bifurcation. They showed that their model exhibits two different kinds of bifurcation: one is transcritical and another is backward bifurcation. Feng et al. [17] observed that their model undergoes backward bifurcation and can sustain multiple endemic equilibria when R0 < 1 provided that the probability of reinfection exceeds a critical value. The main aim of this paper is to study the role of exogenous reinfections and the contact rate on the progression of active TB. The introduction of exogenous reinfection on the transmission dynamics of our classical TB model exhibits interesting outcomes including backward 3

bifurcation. This backward bifurcation indicates that the model system can exhibits multiple endemic equilibrium points for R0 < 1, which implies that the probability of reinfection α crossed the threshold value αcr . Backward bifurcations, can be characterized as multiple endemic equilibria coexists, may give rise to TB transmission persists even though R0 < 1. We also observed that our TB transmission model undergoes Hopf-bifurcation with respect to the exogenous reinfection rate α as well as the contact rate β. The existence of Hof-bifurcation has been counterchecked by the software Matcont (MATLAB package for numerical bifurcation analysis). The importance of Hopf-bifurcation in the context of TB transmission is that at the bifurcation point a limit cycle is formed around the endemic equilibrium point, hence resulting in stable periodic solutions. The significance of periodic solutions is relevant in epidemic models, as it indicates that the disease levels may oscillate around the equilibrium points even in the absence of any treatment strategy. The framework of this paper is as follows. In section 2, we introduce the TB transmission model with exogenous reinfection. Section 3 introduces the qualitative behavior of the model including positivity of solutions, boundedness of the system and the existence of equilibria and their stability analysis. In this section, we calculate the basic reproduction number R0 and the existence of multiple equilibria whenever R0 < 1, which can be characterized as subcritical bifurcation or backward bifurcation. In the section 4, we investigate that our system undergoes Hopf-bifurcation. Also, we find the direction of Hopf-bifurcation and the stability of bifurcating periodic solutions in the Section 4. Section 5 numerically confirms the analytical results for hypothetical set of parameters, along with the effect of exogenous reinfection and the contact rate β. Section 6 discusses the biological significance of our results as well as various concluding remarks. 2. The model In this section, we propose a mathematical model for TB transmission with multiple reinfections, based on epidemiological status. The total host population is divided into four epidemiological classes: susceptibles (S), exposed (E, infected but not infectious), infectious (I) and recovered (R, still susceptible). The total host population is denoted by N, that is, N = S + E + I + R, at any time t. We assume that an individual can be infected only when contacts with infectious individuals. The dynamics of TB infection can be described by the following system of nonlinear ordinary differential equations: dS = dt dE = dt dI = dt dR = dt with initial conditions:

Λ − βIS − δS, (1 − c)βIS − αβEI − kE − δE + γβRI,

(1)

αβEI + kE − δI − hI − dI + cβIS, hI − δR − γβRI,

S(0) > 0, E(0) > 0, I(0) > 0, R(0) > 0. 4

(2)

The total number of susceptible individuals increased by recruitment through births and immigration at a constant rate Λ. Susceptible individuals who come into contact with infected individuals and move to exposed class but they still not infectious. Thus the susceptible individuals decreased due to contact with infected class at a rate β; δ represents the percapita natural death rate. A fraction c (0 < c < 1) of newly infected individuals are assumed to undergo a fast progression due to low resistant, while the remainder are latently infected and moves to the exposed class. Further, a latently infected individual can accrue a new infection due to long latency period k1 . To model this scenario we have taken the term αβEI, where α denotes the level of exogenous re-infection. A value of α ∈ (0, 1) implies that a new infection is likely to be re-infection. Again a recovered individual is not permanently immune from the disease, a fraction γ has taken to model this situation. Here, γ = 0 means the immunity is life long and γ = 1 means no immunity. The infected individuals decreased when individuals recover from TB transmission due to treatment at a rate h and d (> δ) is the per-capita disease induced death rate. Figure 1 describes the schematic representation of the above biological scenarios.

Figure 1: The figure represents the schematic diagram, explaining the dynamics of Tuberculosis transmission.

3. Qualitative behavior of the model 3.1. Positivity Theorem 3.1. Every solution of (1) with respect to the initial conditions (2), defined in [0, +∞), remain positive ∀ t > 0. Proof. The system (1) can be written in vector form as ˙ X(t) = AX(t), where, X(t) = col(S, E, I, R), X(0) = col(S(0), E(0), I(0), R(0), with    A1 (X(t)) Λ − βIS − δS  A2 (X(t))   (1 − c)βIS − αβEI − (k + δ)E + γβRI     A3 (X(t))  =  αβEI + kE − (δ + h + d)I + cβIS A4 (X(t)) hI − δR − γβRI 5

(3) 

 , 

where, A : R4 → R4+ and A ∈ C ∞ (R4 ). It is obvious that in the system (1), Ai (Xi )|Xi =0 ≥ 0, for i = 1, 2, 3, 4. Due to Nagumo’s theorem [11], we can conclude that, the solution of (3) with initial conditions A0 ∈ R4+ , say A(t) = A(t, A0 ) such that A(t) ∈ R4+ ∀ t > 0, that is, for all finite time. 3.2. Boundedness Now, we examine the boundedness of (1), as the solution ’live’ in biological space where all the variables must be nonnegative. Theorem 3.2. Every solution of (1) which initiate in R4+ are uniformly bounded. Proof. Adding all the equations of the model system (1), we obtain dN dt

= Λ − δN − dI, + δN ≤ Λ. ⇒ dN dt

Now, integrating both sides of the above inequality and using the theory of differential inequality [12, 13], we get 0 < N (S, E, I, R) ≤

Λ δ

+ N (0)e−δt ,

where, N (0) = S(0) + E(0) + I(0) + R(0). Taking, t → +∞, we get 0 < N ≤ Λδ . Therefore, all the solution of the system (1) that initiating at {R4+ \ 0} are confined in the region Ω = {(S, E, I, R) ∈ R4+ : 0 < S + E + I + R ≤

Λ + }, δ

for any  > 0 and for t → ∞. Hence the theorem. 3.3. Equilibria and their stability The model system (1) always has the disease free equilibrium point at P0 ( λδ , 0, 0, 0) and the I ∗ {(I ∗ β+δ)(δ+d+h)−cβΛ} λ ∗ endemic equilibrium point P ∗ (S ∗ , E ∗ , I ∗ , R∗ ), where S ∗ = δ+βI , ∗, E = (k+I ∗ αβ)(I ∗ β+δ) ∗ hI R∗ = I ∗ βγ+δ and I ∗ is a positive root(s) of the equation ρ1 I 3 + ρ2 I 2 + ρ3 I + ρ4 = 0,

(4)

where ρ1 ρ2 ρ3 ρ4

= = = =

αβ 3 γ(d + δ), β 2 {(δ + d + h)(αδ + kδ + γδ + αγδ) − γ(αβΛ + αδh + hk)}, δβ(δ + d + h){αδ + (1 + δ)(k + δ)} − β 2 Λ(αδ + kγ + cγδ) − βγδkh, δ 2 (k + δ)(δ + h + d)(1 − R0 ),

and the basic reproduction number R0 is given by    Λ(cδ + k) β R0 = . δ(k + δ) δ+h+d

(5)

This basic reproduction number R0 represents the average number of secondary infectious cases produced by an infectious individual during his/her effective infectious period in a 6

whole susceptible individuals. Clearly, the basic reproduction number R0 is an increasing function of disease transmission rate β and primary progression rate c. It can be noticed that the basic reproduction number R0 is independent of exogenous reinfection rate α. The equation (4) can be written in the standard form as z 3 + 3Hz + G = 0,

(6)

where z = I+

ρ2 3ρ1 ρ3 − ρ22 27ρ21 ρ4 − 9ρ1 ρ2 ρ3 + 2ρ32 , H = , G = . 3ρ1 9ρ21 27ρ31

The strum’s functions can be expressed as p(z) = z 3 + 3Hz + G, p1 (z) = z 2 + H, p2 (z) = − 2Hz + G, p3 (z) = − (G2 + 4H 3 ). Using well-known Strum’s theorem one can easily verify that all the roots of (6) are real and distinct if G2 + 4H 3 < 0. Thus, all the positive roots of (6) yield a nonnegative roots of (4) if, ρ2 ≤ 0.

Now, G2 + 4H 3 < 0 implies that H < 0. Therefore, p(∞) > 0 and pi (∞) > 0, for i=1,2,3. Also it can be noticed that p(0) = G, p1 (0) < 0, p2 (0) = −G and p3 (0) > 0. Therefore, the number of roots of (4) depends on the sign of G. If G < 0 then the cubic equation has exactly one nonnegative root and it has exactly two positive roots if G > 0. Now,ρ4 < 0, that is, R0 > 1 and ρ3 > 0, ρ2 ≤ 0 together implies that G < 0 but note that R0 < 1, ρ2 ≤ 0 and ρ3 > 0 gives G > 0. We can summarize the above results in the following theorem: Theorem 3.3. (a) The system (1) has at least one endemic equilibrium point if R0 > 1. (b) The system (1) has exactly two endemic equilibrium points if G2 +4H 3 > 0, R0 < 1, ρ2 ≤ 0 and ρ3 > 0. Otherwise, the system (1) has exactly one endemic equilibrium point if G2 + 4H 3 < 0, R0 > 1 and ρ2 ≤ 0. (c) The system (1) has no endemic equilibrium point if R0 < 1, ρ2 > 0 and ρ3 > 0. Note: It can be observed from the Theorem 3.3, that the system of equations (1) admits multiple interior equilibrium points when the basic reproduction number is less than 1 and this scenario reveals that the system undergoes a subcritical bifurcation, which is known as backward bifurcation, that is, there is a branch of interior equilibrium points bifurcating backwards from the disease free equilibrium at R0 = 1 [17, 18]. Thus, our model system (1) demonstrates that R0 < 1, may not be sufficient to eradicate the disease. 3.4. Local stability analysis To perform local stability analysis of (1) around each equilibrium point first we have to compute the variational matrix corresponding to each of the equilibrium point(s). The variational matrix of the system (1) at any point (S, E, I, R) is given by the following matrix:   −βI − δ 0 −βS 0  (1 − c)βI −αβI − (k + δ) (1 − c)βS − αβE + γβR  γβI , J =   cβI αβI + k αβE + cβS − (δ + h + d) 0 0 0 h − γβR −δ − γβI 7

3.4.1. Transcritical bifurcation Theorem 3.4. The disease free equilibrium (DFE) point P0 ( Λδ , 0, 0, 0) is locally asymptotically stable for R0 < 1 and unstable for R0 > 1. Proof. The characteristic equation of the Jacobian matrix of (1) at DFE is given by (λ + δ)2 (λ2 + a1 λ + a2 ),

(7)

where a1 = 2δ + h + d + k −

cβΛ , δ

β a2 = (k + δ)(d + h + δ) − Λ(k + cδ) . δ Note that a2 > 0 implies a1 > 0. Hence, all the eigenvalues of the Jacobian matrix of the system (1) at disease free equilibrium point has negative real parts if a2 > 0, that is, if R0 < 1. Therefore, the disease free equilibrium P0 is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. Note: From the above theorem we see that the disease free equilibrium P0 changes its stability from stable to unstable as the basic reproduction number R0 increases through 1, and at R0 = 1 our system (1) experiences a bifurcation around P0 , known as transcritical bifurcation, which has been verified in the next theorem. Theorem 3.5. The system (1) undergoes a transcritical bifurcation around the disease free equilibrium point P0 when the bifurcation parameter R0 = 1. Proof. At R0 = 1, the Jacobian matrix JP0 has one geometrically simple zero eigenvalue with T    (1−c)βΛ h k+δ , , 1, and left eigenvector v = 0, 1, , 0 . Now right eigenvector ω = − βΛ 2 δ δ(k+δ) δ k 

  Λ − βSI − δS f1  (1 − c)βSI − αβEI − (k + δ)E + γβRI   f2 = f (X, β) =   αβEI + kE − (δ + h + d)I + cβSI   f3 hI − δR − γβRI f4

where X = (S, E, I, R) and we have 

 −SI  (1 − c)SI − αEI + γRI  . Dβ f =    αEI + cSI −γRI

From the above expressions, we have

(v(DX Dβ f )w)P0

cδ  Λ 1+ 6= 0, = δ k 8



 , 

(8)

and v((DXX f )(w, w)) = (v

i=4 X

(ei wT DX (DX fi )T w))P0

i=1

=

2hγβ 2β 2 Λ h (1 − c)αδ  1 i + − 1+ 6 0. = δ δ k(k + δ) δ

Hence the system (1) undergoes a transcritical bifurcation at the disease free equilibrium P0 [14]. 3.5. Backward bifurcation In the subsection 3.3, we have seen that our model system (1) exhibits two endemic equilibrium points among which one is stable and other is unstable even though R0 < 1. Backward bifurcations, which are characterized by more than one endemic equilibrium points, allow a disease to persist even though R0 < 1. Thus, it is worthy to prove the backward bifurcation, and in order to do that we use the result of Castillo-Chavez and Song [15]. For simplification and understanding the P4 theory behind backward bifurcation, Twe write S = x1 , E = x2 , I = x3 , R = x4 and N = j=1 xj . By setting X = (x1 , x2 , x3 , x4 ) , the system = F (X) with F = (f1 , f2 , f3 , f4 )T . Thus, we have (1) can be written as dX dt     Λ − βx1 x3 − δx1 f1  (1 − c)βx1 x3 − αβx2 x3 − (k + δ)x2 + γβx3 x4   f2  =  F (X) =  (9)    f3  . αβx2 x3 + kx2 − (δ + h + d)x3 + cβx1 x3 hx3 − δx4 − γβx3 x4 f4 The variational matrix of the system (1) at the disease free equilibrium P0 ( λδ , 0, 0, 0) is   0 −δ 0 − −βΛ δ  0 −(k + δ) (1 − c) βΛ 0  . δ JP0 =  cβΛ  0 k − (δ + h + d) 0  δ 0 0 h −δ

1, we choose β as our bifurcating parameter, then R0 = 1 yields β ∗ = At the threshold value β = β ∗ , the transformed system (9) has one simple zero eigenvalue and other eigenvalues are non-positive. Now, we compute the right and left eigenvectors of JP0 |β=β ∗ , corresponding to the simple zero eigenvalue and these are  T  T (1−c)βΛ h k+δ w = − βΛ , , 1, and v = 0, 1, , 0 respectively. δ2 δ(k+δ) δ δ At R0 = δ(k+δ)(δ+h+d) . Λ(k+cδ)

Now we proceed to obtain the associated backward bifurcation coefficients respectively denoted by a and b as reported in Theorem 4.1 of [15]. Hence, a =

4 X

k,i,j=1

vk wi wj

∂ 2 fk (P0 , β) ∂xi ∂xj

= v2 w1 w3 (1 − c)β + v2 w2 w3 (−αβ) + v2 w4 w3 γβ + v3 w2 w3 αβ + v3 w1 w3 cβ 9

β 2 Λ h α(1 − c) cδ + k i hγβ − + , k k+δ δ2 δ 4 X ∂ 2 fk b = (P0 , β ∗ ) vk wi ∂x i ∂β k,i=1 =

= v2 w3 (1 − c) =

Λ cΛ + v3 w3 δ δ

Λ(k + cδ) . δk

Due to Theorem 4.1 of [15], if both the bifurcation coefficients a and b are positive, then the system (1) demonstrates backward bifurcation in which a stable and an unstable nontrivial equilibrium points coexist. Clearly, b is always positive for all the biologically feasible parameters. Therefore, the existence of backward bifurcation at β = β ∗ , equivalently R0 = 1, . Hence, our system exhibits backward arises if a > 0, that is if α > αmin = (k+cδ)(k+δ) (1−c)δ 2 bifurcation where a stable and a unstable endemic equilibrium point coexist with the stable disease free equilibrium P0 even though R0 < 1. Epidemiologically, it says that only reducing R0 to less than one, can not eradicate the disease from the host.

4. Analysis of Hopf - bifurcation In this section, we shall study the Hopf - bifurcation of the system (1) around the endemic equilibrium point P ∗ (S ∗ , E ∗ , I ∗ , R∗ ). The variational matrix at P ∗ is given by 

 −βI ∗ − δ 0 −βS ∗ 0  (1 − c)βI ∗ −αβI ∗ − k − δ (1 − c)βS ∗ − αβE ∗ + γβR∗  γβI ∗ . A= ∗ ∗ ∗ ∗   cβI αβI + k cβS + αβE − (δ + d + h) 0 ∗ ∗ 0 0 h − γβR −δ − γβI

The characteristic equation of the variational matrix A is given by x4 + A1 x3 + A2 x2 + A3 x + A4 = 0, where A1 A2 A3 A4

= = = =

βm1 − m2 , cS ∗ I ∗ β 2 − ηa32 β + m3 , β 2 m4 − βm5 − m6 , β 2 m7 − βm8 − m9 ,

with m1 = β{(1 + α + γ)I ∗ − cS ∗ − αE ∗ }, m2 = 4δ + k + d + h, m3 = a11 a22 + a11 a44 + a22 a44 − a33 (a11 + a22 + a44 ), 10

(10)

m4 m5 m6 m7 m8 η a11 a22 a32 a33 a43 a44

= = = = = = = = = = = =

S ∗ I ∗ {ca22 + (1 − c)a32 }a44 , a32 {(1 + a44 )η1 + γa43 I ∗ }, a11 (a22 a33 + a33 a44 − a22 a44 ) − a22 a33 a44 , S ∗ I ∗ {ca22 + (1 − c)a32 }a44 , a11 a32 (ηa44 + γa43 a32 I ∗ ) − a11 a22 a33 a44 , (1 − c)S ∗ − αE ∗ + γR∗ , βI ∗ + δ, αβI ∗ + k + δ, αβI ∗ + k, cβS ∗ + αβE ∗ − (δ + d + h), h − γβI ∗ , δ + γβI ∗ .

The Routh-Hurtwitz criteria gives a set of necessary and sufficient conditions for the roots of the equation (10) to have negative real part are as follows: A1 > 0, A3 > 0, A4 > 0, A1 A2 > A3 and A1 A2 A3 − A23 > A21 A4 . From these expressions it is very difficult to interpret the outcomes in terms of TB transmission. So, we will explore the local stability criteria by using numerical simulations. Now we shall derive the transversality conditions for Hopf - bifurcation around the endemic equilibrium point P ∗ by considering β as bifurcating parameter. Theorem 4.1. The endemic equilibrium point P ∗ (S ∗ , E ∗ , I ∗ , R∗ ) undergoes Hopf - bifurcation as β varies over R. Let φ : (0, ∞) → R, be a continuously differentiable function of β can be defined as φ(β) = A1 (β)A2 (β)A3 (β) − A23 (β) − A4 (β)A21 (β). Let β0 be a positive root of φ(β) = 0. Therefore, the system (1) enters into Hopf - bifurcation around P ∗ at β = β0 ∈ (0, ∞) if and only if (i) Ai (β0 ) > 0, i = 1, 2, 3, 4, 5, (ii) A1 (β0 )A2 (β0 ) > A3 (β0 ), (iii) φ(β) = 0, which gives two purely imaginary eigenvalues of the variational matrix A. 0 0 0 0 (iv) A21 (A1 A4 − A2 A3 ) − (A1 A2 − 2A3 )(A1 A3 − A1 A3 ) 6= 0, and all other eigenvalues has purely negative real parts. Proof. By using the condition φ(β0 ) = 0 the characteristic equation (10) takes the following form, 

x2 +

A3  2 A1 A4  x + A1 x + = 0. A1 A3

The roots of the above equation are, xi for i = 1, 2, 3, 4. Let the pair of purely imaginary roots at β = β0 are x1 and x2 , that is, x1 = x¯2 , then we have x3 + x4 = −A1 , 11

(11)

ω02 + x3 x4 = A2 , (12) 2 (13) ω0 (x3 + x4 ) = −A3 , 2 ω0 x3 x4 = A4 , (14) q 3 where ω0 = Im x1 (β0 ). From above ω0 = A . Now, if x3 and x4 are complex conjugate, A1 then from the equation (12), it follows that 2Re(x3 ) = −A1 < 0; if they are real roots then by (12) and (13), x3 , x4 both are negative. Still we need to verify the transversality condition, to complete the discussion. Since φ(β) is continuously differentiable function of all its roots, then there exists an open interval β ∈ (β0 − , β0 + ), where x1 and x2 are complex conjugate for β. Suppose their general forms in this neighborhood are x1 (β) = α1 (β) + iα2 (β), x2 (β) = α1 (β) − iα2 (β). Now we shall verify the transversality condition i d h Re(xj (β)) 6= 0, for j = 1, 2. dβ β=β0

Substituting xj (β) = α1 (β) ± iα2 (β) into the characteristic equation (10) and then after differentiation, separating the real and imaginary parts, we have 0

0

M1 (β)α1 (β) − M2 (β)α2 (β) + M3 (β) = 0, 0 0 M2 (β)α1 (β) + M1 (β)α2 (β) + M4 (β) = 0,

(15)

where M1 (β) M2 (β) M3 (β) M4 (β)

= = = =

4α12 − 12α1 α22 + 3A1 (α12 − α22 ) − 2A2 α1 + A3 , 12α1 α22 − 4α13 + 6A1 α1 α2 + 2A2 α2 , 0 0 0 0 0 A1 α13 − 3A1 α1 α22 + A2 (α12 − α22 ) + A3 α1 + A4 , 0 0 0 0 3A (β)α12 α2 − A1 α23 + 2A2 α1 α2 + A3 α2 .

Now, at β = β0 we obtain M1 (β0 ) = −2A3 , M2 (β0 ) = 2 q   0 0 0 0 A A A1 A3 3 A4 − A2 1 3 and M4 (β0 ) = A A − . 3 A1 A1

q

A3 A1

 A2 −

2A3 A1

 , M3 (β0 ) =

0

Solving the system (15) for α1 (β) at β = β0 , we have

i d h M2 (β0 )M4 (β0 ) + M1 (β0 )M3 (β0 ) 0 Re(xj (β)) = α1 (β0 ) = − dβ M12 (β0 ) + M22 (β0 ) β=β0 0 0 0 0 A21 (A1 A4 − A2 A3 ) − (A1 A2 − 2A3 )(A1 A3 − A1 A3 ) = 6= 0, 2A31 A3 + 2(A1 A2 − 2A3 )2 12

0

0

0

0

if A21 (A1 A4 − A2 A3 ) − (A1 A2 − 2A3 )(A1 A3 − A1 A3 ) 6= 0. Thus the tranversality condition holds and the system undergoes Hopf - bifurcation at β = β0 . Hence, the contact rate β crosses its critial/threshold value, β = β0 , then all the individuals start oscillating around the endemic equilibrium point P ∗ . This finishes the proof of the theorem. 4.1. Direction of Hopf-bifurcation Now we shall derive the direction of Hopf bifurcation, stability and period of bifurcating periodic solutions for system (1), by using the normal form theory and center manifold theorem introduced by Hassard et al. [16]. We used the Poincare’s method to put the system of differential equations (1) into the normal form following the procedure outlined by Hassard et al. [16]. For the sake of simplicity, we replace the current variables to the new variables (n1 , n2 , n3 , n4 ) by considering the following transformations n1 = S − S ∗ , n2 = E − E ∗ , n3 = I − I ∗ , n4 = R − R∗ . Then the system (1) can be written in matrix form as X˙ = AX + B.

(16)

Here A represents the linear part and B stands for the nonlinear part of (1). Where    n1 −βI ∗ − δ 0 −βS ∗ 0  n2   (1 − c)βI ∗ −αβI ∗ − k − δ (1 − c)βS ∗ − αβE ∗ + γβR∗ γβI ∗ , A =  X= ∗ ∗ ∗ ∗  n3   cβI αβI + k cβS + αβE − (δ + d + h) 0 ∗ n4 0 0 h − γβR −δ − γβI ∗ 

 −βn1 n3  (1 − c)βn1 n3 − αβn2 n3 + γβn3 n4  . B=   αβn2 n3 + cβn1 n3 −γβn3 n4

At β = β0 and using the equation φ(β) = A1 (β)A2 (β)A3 (β) − A23 (β) − A4 (β)A21 (β) = 0 in to the characteristic equation (10), it becomes A1 A4  A3  2 x + A1 x + = 0. A1 A3 q 3 From the above equation it is clear that x1,2 = ±i A = ±iω0 and the other eigenvalues A1 have negative real parts, say −ωi for i = 1, 2. Now, we want a transformation matrix P in such a way that reduces the matrix A in the form 

x2 +



 0 −ω0 0 0  ω0 0 0 0  , P −1 AP =   0 0 −ω1 0  0 0 0 −ω2 13



 , 

where the nonsingular matrix P can be written as  1 0 1 1  p21 p22 p23 p24 P =  p31 p32 p33 p34 p41 p42 p43 p44



 , 

where

p22 p22 p23 p24 p31 p32 p33 p34

= = = = = = = =

p41 = p42 = p43 = p44 = l1 = l2 =

a33 a11 + ω02 − cβ 2 I ∗ S ∗ , a11 − a33 , l1 − cβ 2 I ∗ S ∗ , l2 − cβ 2 I ∗ S ∗ , −a11 , ω0 , w1 − a11 , w2 − a11 , a43 (a11 + ω0 ) − , a244 + ω02 ω0 (a44 − a11 ) , a244 + ω02 a43 (w1 − a11 ) , a44 − w1 a43 (w2 − a11 ) , a44 − w2 (a33 + w1 )(a11 − w1 ), (a33 + w2 )(a11 − w2 ).

In order to obtain the normal form for the equation (16), we use another transformation X = P Y , where Y

= (y1 , y2 , y3 , y4 )T .

After some algebraic manipulations, the equation (16) leads us to the following form Y˙ = ΩY + R, where Ω = P −1 AP and

with f is given by



R1 (y1 ,  R2 (y1 , R = P −1 f =   R3 (y1 , R4 (y1 , 

f 1 (y1 ,  f 2 (y1 , f =   f 3 (y1 , f 4 (y1 , 14

y2 , y2 , y2 , y2 ,

y2 , y2 , y2 , y2 ,

y3 , y3 , y3 , y3 ,

y3 , y3 , y3 , y3 ,

 y4 ) y4 )  , y4 )  y4 )

 y4 ) y4 )  , y4 )  y4 )

(17)

where f 1 (y1 , y2 , y3 , y4 ) = −β(y1 + y3 + y4 )(p21 y1 + p22 y2 + p23 y3 + p24 y4 ), f 2 (y1 , y2 , y3 , y4 ) = (1 − c)β(y1 + y3 + y4 )(p21 y1 + p22 y2 + p23 y3 + p24 y4 ) − αβ(p21 y1 + p22 y2 + p23 y3 + p24 y4 )(p21 y1 + p32 y2 + p33 y3 + p34 y4 ) + γ(p31 y1 + p32 y2 + p33 y3 + p34 y4 )(p41 y1 + p42 y2 + p43 y3 + p44 y4 ), 3 f (y1 , y2 , y3 , y4 ) = αβ(p21 y1 + p22 y2 + p23 y3 + p24 y4 )(p31 y1 + p32 y2 + p33 y3 + p34 y4 ) + cβ(y1 + y3 + y4 )(p31 y1 + p32 y2 + p33 y3 + p34 y4 ), 4 f (y1 , y2 , y3 , y4 ) = −γβ(p31 y1 + p32 y2 + p33 y3 + p34 y4 )(p41 y1 + p42 y2 + p43 y3 + p44 y4 ). Equation (17) is the normal form of (16) from which the stability and direction of Hopfbifurcation can be calculated. To find the direction of Hopf bifurcating periodic solution, we need to compute the following quantities at β = β ∗ and origin:

s11 = s02 = s20 = S21 = hj11 = hj20 = λ111 = λ120 = λ320 = j S110 = j S101 =

  ∂R2 ∂R2  1 ∂R1 ∂R1 + +i + , 4 ∂y12 ∂y22 ∂y12 ∂y22    ∂R2 ∂R2 ∂R2 ∂R1  1 ∂R1 ∂R1 − −2 − +2 +i , 4 ∂y12 ∂y22 ∂y1 ∂y2 ∂y12 ∂y22 ∂y1 ∂y2    ∂R2 ∂R2 ∂R2 ∂R1  1 ∂R1 ∂R1 − +2 − −2 +i , 4 ∂y12 ∂y22 ∂y1 ∂y2 ∂y12 ∂y22 ∂y1 ∂y2    ∂R2 ∂R1 1 ∂R1 ∂R2 ∂R2 ∂R1 ∂R1 ∂R2  + + 2 + +i + + 2 + , 8 ∂y13 ∂y23 ∂y1 ∂y2 ∂y1 ∂y22 ∂y13 ∂y23 ∂y1 ∂y2 ∂y1 ∂y22   1 ∂ 2 Rj ∂ 2 Rj + , 4 ∂y12 ∂y22   1 ∂ 2 Rj ∂ 2 Rj ∂ 2 Rj − − 2i , 4 ∂y12 ∂y22 ∂y1 ∂y2 h111 h2 h3 h4 , λ211 = 11 , λ311 = 11 , λ411 = 11 , w1 w2 w3 w4 1 2 h20 h20 , λ220 = , w1 + 2iω0 w2 + 2iω0 h320 h420 , λ420 = , w2 + 2iω0 w2 + 2iω0  2 1   ∂ 2 R2 1 ∂ R ∂ 2 R2 ∂ 2 R1  + +i − , 2 ∂y1 ∂yj ∂y2 ∂yj ∂y1 ∂yj ∂y2 ∂yj    ∂ 2 R2 1 ∂ 2 R1 ∂ 2 R2 ∂ 2 R1  − +i + , j = 1, 2, 3, 2 ∂y1 ∂yj ∂y2 ∂yj ∂y1 ∂yj ∂y2 ∂yj

and s21 = S21 +

4 X

j j (2S110 λj11 + S101 λj20 ).

j=1

15

Furthermore, we can compute s11 , s20 , s02 , s21 using the above expressions. Thus, we can obtain the following quantities described by Hassard et al. [16]. i s i h 1 21 s20 s11 − |s11 |2 − |s02 |2 + , 2ω0 3 2 Re{C1 (0)} , = − Re{λ0 (r0 )} = Re{C1 (0)}, 1 0 = − [Im{C1 (0)} + µ2 Im{λ (β ∗ )}]. v

C1 (0) = µ2 β2 T2

(18)

Theorem 4.2. Recall the expressions in (18), the following results hold: (a) The sign of µ2 determines the direction of the Hopf-bifurcation: if µ2 > 0 (µ2 < 0), then it is supercritical (subcritical) Hopf-bifurcation and the bifurcating periodic solution exits for β > β0 (β < β0 ). (b) β2 indicates the stability of periodic solutions. The periodic solutions are orbitally asymptotically stable (unstable) according β2 < 0 (β2 > 0). (c) Finally, T2 gives the period of the bifurcating periodic solutions, the period increases (decreases) as T2 > 0 (T2 < 0). 5. Numerical simulations To investigate the feasibility of our analytical findings regarding stability and bifurcation of the system (1), we performed some numerical simulations using MATLAB by choosing the following set of parameter values: λ = 18, δ = 0.08, β = 0.28, α = 0.23, k = 0.002, h = 9.71, γ = 0.5, d = 0.2, c = 0.13. For the above set of simulated parameter values we calculate the sensitivity key, that is, basic reproduction number R0 ≈ 0.95297 < 1, and we find that the system has a disease free equilibrium point P0 (225, 0, 0, 0) and two endemic equilibrium points, say, P1u (177.445, 39.0856, 0.0765713, 8.20153) and P1s (31.6219, 134.961, 1.74724, 52.3019). Therefore, it can be noticed that the system (1) has multiple endemic equilibrium point for R0 < 1 and it exhibits that the system undergoes a subcritical bifurcation known as backward bifurcation, which suggests that for R0 < 1 may not be sufficient to eradicate the tuberculosis transmission. The corresponding eigenvalues for P0 are given by −1.86836, −0.08, −0.08, − 0.02064, which indicates that the disease free equilibrium point P0 is asymptotically stable. Also, the eigenvalues for P1u are −1.16486, − 0.0860223, − 0.0795685, 0.304639 and P1s are −0.096592, − 0.125408, − 0.0779184 ± 0.261661i. Thus, the system (1) has one unstable and one stable endemic equilibrium point, even though R0 < 1. In order to investigate the backward bifurcation, we verify that both the backward coefficients a = 371.741 and b = 1395 are nonnegative. To validate our analytical findings we perform numerical simulations for the above set of parameter values with fixed β > β ∗ and the figure 2 shows that, there is a branch of endemic equilibrium points bifurcating backwards from the disease free equilibrium point P0 at R0 = 1. From the figure 2, it is clear 16

that the system (1) has multiple endemic equilibrium points for R0 < 1, the dotted line of the figure 2 describes the unstable branch and the bold line indicates the stable branch of the multiple endemic equilibria. If we slightly increase the value of β = β ∗ = 0.34, then the basic reproduction number R0 become greater that one, that is, R0 ≈ 1.15718 > 1 and the system (1) has only one stable endemic equilibrium point. Thus, the system (1) undergoes forward bifurcation which has been shown in the figure 3. In the section 4, we analyzed that our system (1) undergoes Hopf bifurcation with respect to the contact rate β. From the bifurcation analysis, it is clear that at β = 0.28, the endemic equilibrium point P ∗ undergoes Hopf bifurcation. Using the software MATCONT, we indicate the point at which the system (1) undergoes Hopf bifurcation is H = (119.761076, 76.987406, 0.255897, 27.355879) and the threshold value for the contact rate β = 0.2785 ≈ 0.28. The oscillation of the Hopf bifurcation persists at the critical value β = 0.351579 ≈ 0.35, which has been labeled at the point H = (29.644423, 100.005997, 1.499514, 90.101282). This is shown in the figures 4 A, B, C, D. The importance of Hopf bifurcation in the context of tuberculosis transmission is that at the bifurcation/critical point a limit cycle is formed around the endemic equilibrium point, hence resulting in stable periodic solutions. The occurrence of periodic oscillating solutions is relevant in tuberculosis models, as it indicates that the disease levels may oscillate around the endemic equilibrium point P ∗ even in absence of any treatment. The system (1) exhibits Hopf bifurcation, resulting a period - 1 limit cycle occur at the bifurcation point βcr = 0.28, which has been shown in the figure 5(a) and the period - 1 limit cycle persist at the critical point βcr = 0.35. If we further increase the contact rate β = 0.36, the model system is locally asymptotically stable around the nonnegative endemic equilibrium point P ∗ , which also has been shown in the figure 5(b). Next we try to investigate the model behavior of susceptible individuals (S), exposed indiviuals E, infected individuals (I) and recovered individuals (R) with respect to the different values for the contact rate or the rate of infection parameter β. Let us consider different values for the parameter β, for example, β = 0.28 and β = 0.36. For these two values of β and other parameters are specified earlier we plot the solution curves of susceptible individuals, exposed class, infected individuals and recovered class in figure 6. From the time series solution it can be observed that when the value of β tends very close to βcr , all the populations oscillate in same length of magnitude, which means that a limit cycle occur at the threshold value βcr and the system undergoes Hopf bifurcation. When β = 0.36 crosses the threshold value βcr the system become asymptotically stable which has been shown in the figure 7. We also investigate that our system (1) undergoes Hopf bifurcation with respect to the reinfection level α. It is observed from the bifurcation diagram (see figure 8) that at α = 0.88, the endemic equilibrium point P ∗ undergoes Hopf bifurcation. Using the software MATCONT (a MATLAB package for numerical bifurcation analysis of ODEs) [35], we indicate the point at which the system (1) undergoes Hopf bifurcation is H = (15.930211, 39.795434, 3.888622, 155.664169) and the threshold value for the contact rate α = 0.87646 ≈ 0.88. The importance of Hopf bifurcation in the context of tuberculo17

sis transmission is that at the bifurcation point a limit cycle is formed around the endemic equilibrium point, hence resulting in stable periodic solutions. The system (1) exhibits Hopf bifurcation, resulting a period - 1 limit cycle occur at the bifurcation point α = 0.88, which has been shown in the figure 9(a) and the period - 1 limit cycle persist at the critical point αcr = 0.88. If we further increase the contact rate α = 0.88, the model system is locally asymptotically stable around the nonnegative endemic equilibrium point P ∗ , which shown in the figure 9(b). It can also be noticed that our system (1) undergoes Hopf bifurcation at the threshold value α = 0.2743 ≈ 0.27 labeled as H = (138.337008, 65.463090, 0.185619, 20.550237). This is shown in the figure 8, where H is the neutral saddle, which has no meaning in the biological context. Next we plot the time series solutions to investigate the behavior of susceptible individuals (S), exposed class E, infected individuals (I) and recovered class (R) with respect to the different values for the reinfection rate α. Consider two different values for the parameter α, for example, α = 0.88 and α = 0.90. For these two values for α and other parameters are specified in the text we plot the solution curves of susceptible individuals, exposed class, infected individuals and recovered class in figure 10h. From the time series solution it can be observed that when the value of α tends very close to αcr , all the populations oscillate in same length of magnitude, which means that a limit cycle occurs at the threshold value αcr and the system undergoes Hopf bifurcation. When α = 0.90 crosses the threshold value αcr the system become asymptotically stable which has been shown in the figure 11. The figure 12A indicates that for different values of transmission rate β the system (1) become asymptotically stable, which also shows that the infection pick become high for higher values of β. The figure 12B represents the impact of exogenous reinfection rate α for different values of α. The figure 12B shows that the system (1) becomes asymptotically stable, which also indicates that the infection pick become high for greater values of α.

6. Discussion and conclusions In this study we have proposed and analyzed a TB transmission model through a system of nonlinear differential equations, considering the role of exogenous reinfection and endogenous reactivation among the treated individuals. It is investigated that the basic reproduction number R0 plays a vital role in epidemiology for the extinction and persistence of TB transmission and it is known as threshold value above which the disease persist and below which the disease goes extinction. The basic reproduction number is the exchange of stability of the disease free equilibrium and the system undergoes a transcritical bifurcation at R0 = 1. The introduction of exogenous reinfection into our model system (1) admits qualitative behaviors unexpectedly nonidentical from those derived from model without reinfection. In our study we observed that R0 < 1 is not enough to eradicate the disease from the individuals. Moreover, we noticed that two endemic equilibrium points exist for R0 < 1 among which one is asymptotically stable and other is unstable. To verify this situation we plotted a subcritical bifurcation, known as backward bifurcation due to sufficient level of exogenous reinfection rate and contact rate, which has been shown in the figure 2. Consequently, exogenous reinfection rate α and transmission rate β may in fact assist re-establish TB transmission 18

to a new endemic level via the flux of infectious individuals into the system. We also observed that our system (1) undergoes forward bifurcation (see figure 3) at R0 = 1, which means that the infected population size will be approximately proportional to the difference | R0 − 1 |. From the expression of basic reproduction number R0 , it is clear that R0 is an increasing function of contact rate β. We have described in the theorem 3.3, that the system (1) has multiple, unique or non-existence of endemic equilibrium points when the threshold value of R0 is below one, one or greater than one and the transmission rate β passes through another threshold value. Thus, in our model system (1) the contact rate β plays a vital role for the existence of subcritical/backward bifurcation which has been shown in the figures 2, 3. For our model (1), the transmission rate β and the exogenous reinfection rate α plays an important role. The contact rate β can drive a unstable equilibrium to a stable one; that is, there is a critical value βcr , such that for β > βcr the nonnegative endemic equilibrium E ∗ become stable, and it loses stability for the lower threshold value of βcr . It was shown that the system (1) experiences the Hopf - bifurcation as the parameter β crosses a critical value βcr . Further increase in β beyond the bifurcation point leads to simple dynamic behavior and system become locally asymptotically stable. We have shown this scenario numerically in the figures 4. Also, we have shown that the system (1) undergoes Hopf-bifurcation as the exogenous reinfection rate α crosses some threshold value αcr . For the lower threshold value of αcr the system become unstable and if we gradually increase the value of α the system become asymptotically stable. These scenarios shown numerically in the figures 8. The occurrence of Hopf bifurcation in the context of tuberculosis transmission is that at the bifurcation point a limit cycle is formed around the endemic equilibrium point, hence resulting in stable periodic solutions. The importance of periodic solution is relevant in epidemic models, as it indicates that the disease levels may oscillate around the equilibrium points even in the absence of any treatment or control. The main thrust of this article is to study the role of transmission rate β and the rate of exogenous reinfection α in the dynamics of epidemic model. To achieve our goal, we investigate a simple TB transmission model consisting of susceptible, exposed, infectious and recovered class. Our results suggest that for the lower threshold value of the contact rate β and the exogenous reinfection rate α, the infection pick is low and for the higher values of β and α the infection pick is high (see figures 12), which indicates that the effectiveness of β and α is very helpful for the control of TB transmission. Our results also suggest that the parameters β and α can destabilize the system (1) through Hopf-bifurcation (see figures 4 and 8).

References [1] B.R. Bloom, Tuberculosis: Pathogenesis, Protection, and Control, ASM Press, Washington, D.C., 1994. [2] S.M. Arend, J.T. van Dissel, Evidence of endogenous reactivation of tuberculosis after a long period of latency, J. Infect. Dis. 186(6) (2002) 876–877. 19

3.5

Stable Steady State

3 2.5

I

2 1.5 1

Unstable Steady State

0.5 0 0

0.5

1

1.5

2

2.5

3

3.5

R0 Figure 2: The figure shows the backward bifurcation of endemic equilibrium points, I stands for the number of infective individuals. The dotted line represents unstable steady state whereas the solid line represents the stable steady state.

[3] B. Miller, Preventive therapy for tuberculosis, Med. Clin. N. Am. 77(6) (1993) 1263–1275. [4] E. Vynnycky, P. Fine, The long-termdynamics of tuberculosis and other diseases with long serial: The implications of and for changing reproduction numbers, Epidemiol. Infect. 121 (1998) 309–324. [5] K. Styblo, Epidemiology of Tuberculosis: Selected Papers, Royal Netherlands Tuberculosis Association, (1991). [6] T. Lillebaek, A. Dirksen, I. Baess, B. Strunge, V.O. Thomsen, A.B. Andersen, Molecular evidence of endogenous reactivation of Mycobacterium tuberculosis after 33 years of latent infection, J. Infect. Dis. 185(3) (2002) 401–404. [7] J.W. Raleigh, R.H. Wichelhausen, Exogenous reinfection with mycobacterium tuberculosis confirmed by phage typing, Am. Rev. Respir. Dis. 108 (1973) 636–642. [8] E. Nardell, B. Mc Innis, B. Thomas, S. Weidhaas, Exogenous reinfection with tuberculosis in a shelter for the homeless, N. Engl. J. Med. 315 (1986) 1570–1573. [9] I. Sutherland, E. Svandova, Endogenous reactivation and exogenous reinfection: Their relative importance with regard to the development of non-primary tuberculosis, Bull. Int. Un. Tuberc. 46 (1971) 74–114. [10] C. Castillo-Chavez, Z. Feng, Global stability of an age-structure model for TB and its applications to optimal vaccination strategies, Math. Biosci. 151 (1998) 135–154. 20

2

Stable Steady State

I

1.5

1

0.5

0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

R0 Figure 3: The figure shows the forward bifurcation of endemic equilibrium point with respect to β = 0.34 and the bifurcation diagram demonstrates the exchange of stability at R0 = 1. All other parameters are same as in figure 2.

[11] N. Nagumo, ber die lage der integralkurven gewhnlicher differentialgleichungen, Proc. Phys.-Math. Soc. Japan. 24 (1942) 551. [12] G. Birkhoff, G.-C. Rota, Ordinary Differential Equations, Ginn, Boston, 1982. [13] S. Khajanchi, Uniform Persistence and Global Stability for a Brain Tumor and Immune System Interaction, Biophys. Rev. Lett. 14(4) (2017) 187–208. [14] J. Sotomayor, Generic bifurcations of dynamical systems, In: Peixoto, M.M. (Ed.), Dynamical Systems, Academic Press, New York, 1973. [15] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1(2) (2004) 361–404. [16] B.D. Hassard, N. Kazarinoff, Y. Wan, Theory and Application of Hopf Bifurcation, CUP archive, 1981. [17] Z. Feng, C. Castillo-Chavez, A.F. Capurro, A model for tuberculosis with exogeneous reinfection, Theor. Popul. Bio. 57 (2000) 235–247. [18] T.K. Kar, P.K. Mondal, Global dynamics of a tuberculosis epidemic model and influence of backward bifurcation, J. Math. Model. Algor. 11 (2012) 433–459.

21

200

120

A

B H

100 150 80

E(t)

S(t)

LP H 100

H LP

60 40

50

H 20 0 0.2

0.25

0.3

0.35

0.4

0.2

2

1.6

0.35

0.4

H

70

1.2

60

1

50

0.8

40

0.6

30

0.4

20

H

H LP

10

LP 0.25

D

80

1.4

0.2

H

90

R(t)

I(t)

0.3

100

C

1.8

0.2

0.25

0.3

β

0.35

0.4

0.2

0.25

0.3

β

0.35

0.4

Figure 4: The figure shows the points of bifurcation for the system parameter β. The Hopf bifurcation occur at the points β = 0.28 and β = 0.36, which has been labeled by “H”. At β = 0.28 the Hopf bifurcation labeled by “H” is neutral saddle, which has no meaning in the biological context. When the contact rate β crosses the threshold value βcr = 0.36, the system (1) become asymptotically stable. The parameter values as specified in the numerical section.

[19] P.K. Mondal, T.K. Kar, Optimal treatment control and bifurcation analysis of a tuberculosis model with effect of multiple re-infections, Int. J. Dynam. Control. (2015) 1–14. doi.10.1007/s40435-015-0176-z. [20] S. Khajanchi, D. Ghosh, The combined effects of optimal control in cancer remission, Appl. Math. Comput. 271 (2015) 375–388. [21] T.K. Kar, A. Batabyal, Stability analysis and optimal control of an SIR epidemic model with vaccination, Biosystems 104 (2011) 127–135. [22] S. Bowong, J.J. Tewa, Global analysis of a dynamical model for transmission of tuberculosis with a general contact rate, Commun. Nonlinear Sci. Numer. Simulat. 15 (2010) 3621– 3631. [23] D.J. Gerberry, Practical aspects of backward bifurcation in a mathematical model for tuberculoisis, J. Theor. Biol. 388 (2016) 15–36. 22

3 6 2

I

I

4 1

2 0 120 100 100

E

50 80

0

0 120 100

100

S

50

80

E

(a) β = 0.28

60

S

0

(b) β = 0.36

Figure 5: The figures shows the 3D phase diagram for the system (1) with respect to the contact rate β. The figure (a) exhibits the period - 1 limit cycle arising at the Hopf bifurcation point β = 0.28 and the figure (b) represents that the system (1) is locally asymptotically stable when β crosses the threshold value βcr = 0.36. All other parameters are same as specified in the numerical section.

[24] H.F. Huo, M.X. Zou, Modelling effects of treatment at home on tuberculosis transmission dynamics, Appl. Math. Model. 40 (2016) 9474–9484. [25] J.A. Romeyn, Exogenous Dis. 101 (1970) 923–927.

re-infection

in

tuberculosis,

Am.

Rev.

Respir.

[26] P.W. Uys, P.D. van Helden, J.W. Hargrove, Tuberculosis reinfection rate as a proportion of total infection rate correlates with the logarithm of the incidence rate: a mathematical model, J. R. Soc. Interface. 6 (2009) 11–15. [27] A. van Rie, R. Warren, M. Richardson, T.C. Victor, R.P. Gie, D.A. Enarson, N. Beyers, P.D. van Helden, Exogenous reinfection as a cause of recurrent tuberculosis after curative treatment, N. Engl. J. Med. 341 (1999) 1174–1179. [28] M. Martcheva, H.R. Thieme, Progression age enhanced backward bifurcation in an epidemic model with super-infection, J. Math. Biol. 46 (2003) 385–424. [29] Z. Wang, C.T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salath, D. Zhao, Statistical physics of vaccination, Phys. Rep. 664 (2016) 1–113. [30] Z. Wang, Y. Moreno, S. Boccaletti, M. Perc, Vaccination and epidemics in networked populationsAn introduction, Chaos Soliton Fract. 103 (2017) 177– 183. 23

140 150

120

E(t)

S(t)

100 100

50

80 60 40

0

500

1000

1500

0

5

100

4

80

3

R(t)

I(t)

0

2

500

1000

1500

1000

1500

60 40

1 20 0

0

500

1000

1500

0

Time

500

Time

Figure 6: Waveform plot for the system (1) with respect to the contact rate β = 0.28 around the interior equilibrium point P ∗ . The figures depicts the existence of periodic solutions around P ∗ . The parameter values are specified in the numerical section.

[31] D. Ghosh, S. Khajanchi, S. Mangiarotti, F. Denis, S.K. Dana, C. Letellier, How tumor growth can be influenced by delayed interactionsbetween cancer cells and the microenvironment?, BioSystems 158 (2017) 17–30. [32] H.M. Yang, S.M. Raimundo, Assessing the effects of multiple infections and long latency in the dynamics of tuberculosis, Theor. Biol. Med. Model. 7 (2010) 41. [33] R.M. Anderson, R.M. May, Infectious diseases of humans, Vol. 1, Oxford university press Oxford, 1991. [34] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42(4) (2000) 599– 653. [35] S. Khajanchi, S. Banerjee, Quantifying the role of immunotherapeutic drug T11 target structure in progression of malignant gliomas: Mathematical modeling and dynamical perspective, Math. Biosci. 289 (2017) 69–77.

24

110

80

100

60

90

E(t)

S(t)

100

40 20

70

0

200

400

600

60

800

6

110

5

100

4

90

R(t)

I(t)

0

80

3

70

1

60 0

200

400

600

50

800

Time

200

400

600

800

0

200

400

600

800

80

2

0

0

Time

Figure 7: The figures depict local asymptotic stability of the endemic equilibrium for the system (1) with respect to the contact rate β = 0.36. All other parameters are specified in the numerical section.

[36] P.V.D. Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compermental models of disease transmission. Math. Biosci. 180 (2002) 29–48. [37] B.H. Singer, D.E. Kirschner, Influence of backward bifurcation on interpretation of R0 in a model of epidemic tuberculosis with reinfection, Math. Biosci. Eng. 1(1) (2004) 81–93. [38] P.V.D. Driessche, J. Watmough, A simple SIS epemic model with a backward bifurcation, J. Math. Biol. 40 (2000) 525–540.

25

250

100

200

80

A 150

60

E(t)

S(t)

LP

B

LP H

H

H 40

100

20

50

H 0 0.2

0.4

0.6

0.8

0 0.2

1

4.5

0.8

1

H

160

3.5

140

C

3

120

R(t)

2.5 2

D

100 80

1.5

60

1

40

LP

20

H

0.5

0.6

180

H

4

I(t)

0.4

LP

0 0.2

H 0.4

0.6

0.8

0 0.2

1

0.4

0.6

0.8

1

α

α

Figure 8: The figure depict the points of bifurcation for the system parameter α. The Hopf bifurcation occur at the points α = 0.88 and α = 0.27, which has been labeled by “H”. At α = 0.27 the Hopf bifurcation labeled by “H” is neutral saddle, which has no meaning in the biological context. When the reinfection rate α crosses the threshold value αcr = 0.88, the system (1) become asymptotically stable. The parameter values are specified in the numerical section.

26

80

60

60

40

40

I

I

80

20

20

0 100

0 100 200

200

E

50

100 0

50

S

100

E 0

0

0

S

(b) α = 0.90

(a) α = 0.87

Figure 9: The figure depicts the 3D phase diagram of the system (1) with respect to the reinfection rate α. The figure (a) represents the stable limit cycle oscillations for α = 0.87, whereas figure (b) indicates the local asymptotic stability of the system (1) for α = 0.90. All other parameters are specified in the numerical section.

27

80 70

150

100

E(t)

S(t)

60 50 40 30

50

20 0

0

500

1000

10

1500

0

500

1000

1500

500

1000

1500

200

60 50

150

R(t)

I(t)

40 30

100

20 50

10 0

0

500

1000

1500

Time

0

Time

Figure 10: The plot represents the existence of periodic solutions around the endemic equilibrium point P ∗ for the model system (1), with respect to the reinfection rate α = 0.88, where the other parameter values are specified in the numerical system.

28

70 150

60

E(t)

S(t)

50 100

40 30

50

20 0

0

200

400

600

10

800

0

200

0

200

400

600

800

400

600

800

60 50

150

R(t)

I(t)

40 30

100

20 50

10 0

0

200

400

600

800

Time

Time

Figure 11: The figure depicts local asymptotic stability of (1) around the endemic equilibrium point P ∗ , with respect to the reinfection rate α = 0.90, where the other parameters are specified in the numerical section.

29

35 30 25

I(t)

30

β = 0.18 β = 0.20 β = 0.23 β = 0.25 β = 0.27 β = 0.30

A

20

α = 0.2 α = 0.3 α = 0.4 α = 0.5 α = 0.6 α = 0.7 α = 0.8

B 25

20

15 15 10 10 5

5 0

0

0.5

1

1.5

2

2.5

3

0

3.5

Time

0

1

2

3

4

5

6

Time

Figure 12: The figure A) represents the effect of contact rate β and figure B) represents the effect of exogenous reinfection rate α on the infected population where all other parameters are specified in the numerical section.

30