Mathematics in Applied Sciences and Engineering https://doi.org/10.5206/mase/10221 Volume 1, Number 1, March 2020, pp.39-49 http://www.uwo.ca/lib/mase ANALYSIS OF SOLUTIONS AND DISEASE PROGRESSIONS FOR A WITHIN-HOST TUBERCULOSIS MODEL WENJING ZHANG, FEDERICO FRASCOLI, AND JANE M HEFFERNAN Abstract. Mycobacterium tuberculosis infection can lead to different disease outcomes, we analyze a within-host tuberculosis infection model considering interactions among macrophages, T lymphocytes, and tuberculosis bacteria to understand the dynamics of disease progression. Four coexisting equilibria that reflect TB disease dynamics are present: clearance, latency, and primary disease, with low and high pathogen loads. We also derive the conditions for backward and forward bifurcations and for global stable disease free equilibrium, which affect how the disease progresses. Numerical bifurcation analysis and simulations elucidate the dynamics of fast and slow disease progression. 1. Introduction Mycobacterium tuberculosis (Mtb) is a bacterium that causes an ancient and deadly infectious disease in humans, called tuberculosis (TB) [9]. Currently, TB affects approximately one third of the world’s population [10, 6]. In 2018, the World Health Organization (WHO) estimated approximately 10 million infections globally, and 1.2 million deaths among HIV-negative people [12]. It has also been found that TB susceptibility and disease are increased in HIV-AIDS infected individuals, resulting in higher mortality rates [8, 1, 14, 16]. The pathological outcomes of TB infection include clearance, latent infection, and primary disease with fast or slow progression [13]. After initial infection, 5−10% of infected subjects can clear the disease. Of the remaining individuals, 5 − 10% will progress to primary disease, and the rest will remain in a latently infected state with no clinical symptoms, with the possibility of re-activation to primary disease later in their life. A large number of mechanisms have been proposed to explain TB disease progression considering individual factors, including bacterial and immune response mechanisms. However, the most influential factors for TB outcomes are not currently known. Motivated by this, we analyze a TB host-pathogen model first proposed in Ref. [3]. The model incorporates known mechanisms of host-pathogen interaction in TB dynamics, and includes all realistic disease outcomes. Analysis is performed to determine the driving factors behind disease progression and outcome, especially fast or slow progression to primary disease. The paper is organized as follows. In Section 2, we introduce the established tuberculosis progression model. In Section 3, model dynamics are shown through the proofs of the well-posedness of solutions, the existence of equilibrium solutions, and analyses of the disease free equilibrium. The basic reproduction number R0 and the vector field on the center manifold for the disease free equilibrium when R0 = 1 are derived analytically. The conditions for the occurrence of the backward and forward bifurcations Received by the editors 13 February 2020; revised 4 March 2020; accepted 5 March 2020; published online 7 March 2020. 2000 Mathematics Subject Classification. Primary 37G15; Secondary 92-08. Key words and phrases. Tuberculosis, Disease progression, Stability, Bifurcation. W. Zhang was supported in part by Texas Tech University New Faculty Startup Fund; J. M. Heffernan was funded by NSERC and supported by a York Research Chair program. 39 40 W. ZHANG, F. FRASCOLI, AND J. M. HEFFERNAN are also derived. In Section 4, numerical continuations are carried out for the infected equilibrium to confirm the existence of a backward bifurcation. The corresponding numerical simulations show the fast and slow disease progressions to latency and primary diseases. Finally, conclusions are drawn in Section 5. 2. Model The 4-dimensional model (2.1) includes the MTb ideal target cell population, macrophages (their uninfected Mu and infected Mi populations). It also includes the Mtb bacterial population B, and a population of CD4 T cells, which aid in TB clearance. The model is as follows: dMu dt = sM −µMMu −βMuB dMi dt = βMuB − bMi −γMi T/Mi T/Mi + c dB dt = δB ( 1 − B K ) + Mi ( N1b + N2γ T/Mi T/Mi + c ) −MuB(η + N3β) dT dt = sT + cMMiT eMT + 1 + cBBT eBT + 1 −µTT. (2.1) Briefly, uninfected macrophages Mu enter the system with constant rate sM , and can die naturally (µM ), or be infected by the pathogen B (βMuB). It is assumed that infected macrophages can release new bacteria into the system in two different ways: (1) through cell death and bursting b, producing N1 new bacteria, and (2) through cytotoxic T-lymphocyte killing (represented by the ratio T/Mi) with rate γ and saturating factor c, which releases N2 new bacteria into the system. It is assumed that the bacteria population can divide (δB(1 − B/K) and that bacteria can be lost due to interaction with macrophages. This occurs through immune system neutralization ηBMu or macrophage infection βBMu involving, on average, N3 individual bacteria. Finally, it is assumed that T-cells are produced at a constant rate sT by the thymus, can be stimulated to proliferate through interactions with the infected macrophage cMMiT/(eMT + 1) and bacteria cBBT/(eBT + 1), and can die naturally, with rate µT . Infection is initiated with an initial pathogen load. We refer the reader to Du et al. [3] for more detail on the biology and model assumptions. Parameters and their values are listed in Table 1. In previous work, Du et al. [3] found four biologically realistic equilibria and determined the basic reproduction number. Note that, in the original contribution, there is no mention of the driving factors behind the different outcomes of disease (namely, clearance, latency, and primary disease with fast or slow progression) and only an asymptotic version of the model that neglects the effects of the CD4 T-cell population is used/analyzed. In the following, we expand and elaborate on the four disease outcomes and other interesting aspects of the model using the full model system Eq. 2.1. 3. Model Dynamics 3.1. Well-posedness of solutions. Let D = { (Mu, Mi, B, T) ∈ R4+ : Mu + Mi ≤ Mmax, B ≤ Bmax, T ≤ Tmax } , where Mmax = sM min{µM, b} , Tmax = 1 µT ( sT + cM eM Mmax + cB eB Bmax ) , Bmax = K 2 + √ (4KδMmax(N1b + N2γ) + K2δ2 2δ . (3.1) ANALYSIS OF SOLUTIONS AND DISEASE PROGRESSIONS FOR A WITHIN-HOST TUBERCULOSIS MODEL 41 Proposition 3.1. Under the flow of (2.1), there exists a positive invariant set D that attracts all solutions in R4+ as time moves forward. Proof. The smoothness of the right hand side of model (2.1) guarantees the local existence and unique- ness of the solution of the initial value problem of model (2.1). The trajectories starting from positive initial values never cross the boundary of R4+, since dMu dt |Mu=0 = sM > 0, dMi dt |Mi=0 = β Mu B ≥ 0, and dB dt |B=0 = Mi ( N1b + N2γ T T + cMi ) ≥ 0, dT dt |T=0 = sT > 0. Next, we show that positive solutions are bounded. Due to the positiveness, we have d dt (Mu + Mi) < sM −µMMu − bMi ⇒ lim t→+∞ sup(Mu + Mi)(t) = sM min{µM, b} := Mmax. Moreover, dB dt < δB ( 1 − B K ) + Mi (N1b + N2γ) −MuB(η + N3β), TT+cMi ∈ (0, 1) < − δ K B2 + δB + Mmax (N1b + N2γ) ⇒ B(t) = K/2 + tanh [√ (4KδMmax(N1b + N2γ) + K2δ2)(C0 + t)/(2K) ] × √ (4KδMmax(N1b + N2γ) + K2δ2/(2δ), where C0 is determined by initial condition and C0 + t > 0 for sufficiently large t. We have B(t) = K 2 + √ (4KδMmax(N1b + N2γ) + K2δ2 2δ := Bmax. Then, the last equation in (2.1) satisfies dT dt < sT + cMMmaxT eMT + 1 + cBBmaxT eBT + 1 −µTT < sT + cM eM Mmax + cB eB Bmax −µTT. It hence follows that T(t) < 1 µT ( sT + cM eM Mmax + cB eB Bmax ) := Tmax, and the proposition is proven. � 3.2. Equilibrium Solutions. Denote model (2.1) as M′u = f1, M ′ i = f2, B ′ = f3, T ′ = f4. The corresponding steady states are derived as follows: f1 = 0;⇒ M̄u(B) = sM βB + µM . (3.2) Case 1: If (b + γ)Mi −βMuB 6= 0 or βsMB − (b + γ)(βB + µM )Mi 6= 0, we have f2 = 0 ⇒ T̄(M̄u) = cMi [ γMi (b + γ)Mi −βM̄uB − 1 ] (3.2) ==⇒ T̄(B) = [βsMB − (βB + µM )bMi]cMi βsMB − (b + γ)(βB + µM )Mi , T̄(B) > 0 if βsMB < (βB + µM )bMi or βsMB > (b + γ)(βB + µM )Mi. (3.3) 42 W. ZHANG, F. FRASCOLI, AND J. M. HEFFERNAN Considering the preceding results (3.2) and (3.3), we obtain f3 = cγM 2 i f3af3b = 0, where f3a = (b + γ)(βB + µM )Mi −βsMB f3b = Kb (βB + µM ) Mi + ((βB + µM )δ + sM [(N2 −N3)β −η])KB −B2δ(βB + µM ). The existence of T̄ in (3.3) implies f3a 6= 0. Further, Mi = 0 induces that f4(Mu = M̄u(B), Mi = 0, B, T̄(B) = 0) = sT 6= 0. This indicates that the equilibrium does not exist. Therefore f3 = 0 only implies f3b = 0 followed by M̄i(B) = ( Bδ K − δ + sM (N2 −N3)β −sMη βB + µM ) B b(N1 −N2) , M̄i(B) > 0 if Bδ K > δ + sM (N2 −N3)β −sMη βB + µM and N1 > N2. (3.4) The B in (3.4) satisfies f4(M̄u(B), M̄i(B), B, T̄(B)) = 0, the following is true: F(B) = −eBeMµT T̄3(B) + [ (cMM̄i(B) + eMsT −µT )eB + eM (cBB −µT ) ) T̄2(B) + [ cBB + cMM̄i(B) + eBsT + eMsT −µT ] T̄(B) + sT = 0. (3.5) Then, we find the infected equilibrium E∗ = (M̄u(B), M̄i(B), B, T̄(B)). We note that there could be more than one solution, and up to three feasible infected equilibria. Case 2: If βsMB − (b + γ)(βB + µM )Mi = 0, we have f2 = 0 ⇒ M̄i0 = βsMB (b + γ)(βB + µM ) . (3.6) Then substituting M̄u(B) in (3.2) and M̄i0 in (3.6) into f3(M̄u(B), M̄i0) = 0, yields B̄0 = 0. (3.7) This is followed by f4(M̄u(B), M̄i0, B̄0) = 0, which yields T̄0 = sT µT . (3.8) We thus find the disease free equilibrium (DFE) E0 = (M̄u(B̄0), M̄i0(B̄0), B̄0, T̄0), where M̄u(B̄0) = sM/µM and M̄i0(B̄0) = 0. 3.3. Analysis of the disease free equilibrium. 3.3.1. Calculation of the basic reproduction number. Following the next-generation matrix approach in Ref. [11], the basic reproduction number R0 is the spectral radius of FV −1, where FV −1 =   0 βsMµM N1b + N2γ δ  [b + γ 0 0 sM µM (N3β + η) ]−1 =   0 β N3β + η N1b + N2γ b + γ µMδ sM (N3β + η)   , and R0 = ρ(FV −1) = δµM 2sM (N3β + η) + 1 2 [ δ2µ2M s2M (N3β + η) 2 + 4β(N1b + N2γ) (N3β + η)(b + γ) ]1/2 . (3.9) ANALYSIS OF SOLUTIONS AND DISEASE PROGRESSIONS FOR A WITHIN-HOST TUBERCULOSIS MODEL 43 The Jacobian matrix of model (2.1) at the disease free equilibrium is: J0 =   −µM 0 − βsM µM 0 0 −b−γ βsM µM 0 0 N1b + N2γ δ − sM µM (N3β + η) 0 0 cMsT eMsT + µT cBsT eBsT + µT −µT   , (3.10) and gives the following characteristic equation (z + µT ) (z + µM ) ( z2 + Pz + Q ) = 0, (3.11) where P = b + γ −δ + sM µM (N3β + η) , Q = [(−N2 + N3)γβ − b(N1 −N3)β + η(b + γ)] sM µM − δ(b + γ). Equation (3.11) admits at least two negative roots, z = −µT and z = −µM . The third root, z = δ − b−γ − sM µM (N3β + η), is negative if b + γ + sM µM (N3β + η) > δ. The last root is zero, if [(−N2 + N3)γβ − b(N1 −N3)β + η(b + γ)] sM µM − δ(b + γ) = 0, (3.12) which is equivalent to R0 = 1. Theorem 3.1. Under the condition b+γ + sM µM (N3β + η) > δ, the disease free equilibrium E0 is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. 3.3.2. Existence of a backward bifurcation. Following Theorem 4.1 in Ref. [2], we first shift the disease free equilibrium to the origin by letting x1 = Mu − sMµM , x2 = Mi − 0, x3 = B − 0, x4 = T − sT µT , and φ = β −βT . Here R0(βT ) = 1 and βT = (−δµM + ηsM )(b + γ) sMγ(N2 −N3)γ + sMb(N1 −N3) . Then we compute the approximated center manifold for the system near the origin with one simple zero eigenvalue at R0 = 1, and three negative eigenvalues. We choose a right eigenvector associated with the simple zero eigenvalue, w, and the left eigenvector, v, satisfying vw = 1 as follows: w = 1 n   (δµM −ηsM )(b + γ) µ2Mw̃ δµM −ηsM µMw̃ 1 sT {[(w̃cB − δcM )µT + (eMw̃cB −δcMeB)sT ] µM + ηcMsM (eBsT + µT )} (eBsT + µT )(eMsT + µT )µMµTw̃   , v = [ 0, N1b + N2γ b + γ , 1, 0 ] , 44 W. ZHANG, F. FRASCOLI, AND J. M. HEFFERNAN where n = (N2 −N3)µMγ + [(N1 + N2 − 2N3)b−N2δ]µMγ + N2sMηγ µM (b + γ)w̃ + (N1 −N3)b2µM −N1δbµM + N1sMηb µM (b + γ)w̃ and w̃ = (N2 −N3)γ + b(N1 −N3). Further, the flow of the center manifold y(t) truncated at the quadratic term is written as ẏ = Ay2 + Bφy, (3.13) A = v 2 [ w′ ( ∂f1 ∂xi∂xj )∣∣∣ E0 w, w′ ( ∂f2 ∂xi∂xj )∣∣∣ E0 w, w′ ( ∂f3 ∂xi∂xj )∣∣∣ E0 w, w′ ( ∂f4 ∂xi∂xj )∣∣∣ E0 w ]′ = An Ad , An = ((Ã− [µTcsM (N1 −N2)]bγ)Kδ −sM [(N2 −N3)γ + b(N1 −N3)]2(b + γ)sT )µ2Mδ, −sMδKηµMÃ + KµTγbcs3Mη 2(N1 −N2), Ã = −sT (N2 −N3)γ3 − bsT (N1 + 2N2 − 3N3)γ2 − b3sT (N1 −N3) − (2N1 + N2 − 3N3)sTb2γ + 2(N1 −N2)µTcsMbγ, Ad = sMµ 2 M [(N2 −N3)γ + b(N1 −N3)] 2(b + γ)sTK, and B = v ( ∂fi ∂xi∂β ) E0 w = [(N2 −N3)γ + b(N1 −N3)]sM (b + γ)µM , where i, j = 1 . . . 4. The non-zero terms in ( ∂fk ∂xi∂xj )∣∣∣ E0 , where i, j, k = 1 . . . 4, are ∂f1 ∂x1∂x3 ∣∣∣ E0 = ∂f1 ∂x3∂x1 ∣∣∣ E0 = − ∂f2 ∂x1∂x3 ∣∣∣ E0 = − ∂f2 ∂x3∂x1 ∣∣∣ E0 = (δµM −ηsM )(b + γ) [(N2 −N3)γ + b(N1 −N3)]sM , ∂f3 ∂x1∂x3 ∣∣∣ E0 = ∂f3 ∂x3∂x1 ∣∣∣ E0 = N3 ∂f1 ∂x1∂x3 ∣∣∣ E0 −η, ∂f2 ∂x2∂x2 ∣∣∣ E0 = 2γµTc sT , ∂f3 ∂x2∂x2 ∣∣∣ E0 = −N2 2γµTc sT , ∂f3 ∂x3∂x3 ∣∣∣ E0 = −2 δ K , ∂f4 ∂x2∂x4 ∣∣∣ E0 = ∂f4 ∂x4∂x2 ∣∣∣ E0 = cMµ 2 T (eMsT + µT )2 , ∂f4 ∂x3∂x4 ∣∣∣ E0 = ∂f4 ∂x4∂x3 ∣∣∣ E0 = cBµ 2 T (eBsT + µT )2 . ANALYSIS OF SOLUTIONS AND DISEASE PROGRESSIONS FOR A WITHIN-HOST TUBERCULOSIS MODEL 45 Theorem 3.2. Under the condition B > 0, we have Ad > 0. Then the model (2.1) at the disease free equilibrium E0, when R0 = 1 undergoes (1) a backward bifurcation if An > 0 and (2) a forward bifurcation if An < 0. Furthermore, An(c = 0, β = βT , b = bb) = 0, where bb = γ [(N2 −N3)sM + Kδ]µM −ηKsM [(N3 −N1)sM −Kδ]µM + ηKsM . (3.14) 3.3.3. Global stability analysis for the disease free equilibrium E0. Proposition 3.1 shows that state variables Mu, Mi, B, and T are bounded for sufficiently large time. That is, there exists a time T > 0 such that Mu < Mmax, Mi < Mmax, B < Bmax, and T < Tmax. Applying the “fluctuation lemma” [5], there exists time sequences τn →∞ and σn → +∞ such that M∞i := lim supt→∞Mi(t) = limn→+∞Mi(τn) and limn→+∞ dMi(τn) dt = 0, B∞ := lim supt→∞B(t) = limn→+∞B(σn) and limn→+∞ dB(σn) dt = 0. (3.15) The preceding equations are followed by βMu(τn)B(τn) − bM∞i −γM ∞ i T(τn) T(τn) + cM ∞ i = 0 =⇒ βMu(τn)B(τn) = ( b + γ T(τn) T(τn) + cM ∞ i ) M∞i ≤ βMmaxB ∞ −→ M∞i ≤ βMmaxB ∞ b + γ T(τn) T(τn) + cM ∞ i < β b MmaxB ∞, (3.16) and δB∞ ( 1 − B ∞ K ) + Mi(σn) ( N1b + N2γ T(σn) T(σn)+cMi ) −Mu(σn)B∞(η + N3β) = 0 =⇒ −δB∞ ( 1 − B ∞ K ) + Mu(σn)B ∞(η + N3β) = Mi(σn) ( N1b + N2γ T(σn) T(σn)+cMi ) ≤ Mi(σn) (N1b + N2γ) =⇒ [Mu(σn)(η + N3β) − δ] B∞ ≤ Mi(σn) (N1b + N2γ) =⇒ B∞ ≤ Mi(σn) N1b + N2γ Mu(σn)(η + N3β) − δ ≤ N1b + N2γ Mu(σn)(η + N3β) −δ M∞i , (3.17) where T(τn) T(τn)+cM ∞ i ≤ 1. Subsequently, M∞i ≤ β b MmaxB ∞ ≤ β b Mmax N1b + N2γ Mu(σn)(η + N3β) − δ M∞i . (3.18) If β b Mmax N1b + N2γ Mu(σn)(η + N3β) −δ ≤ 1, or equivalently Mu(σn) ≥ 1 η + N3β [ βsM bµM (N1b + N2γ) + δ ] := Mmaxu , (3.19) then M∞i = 0, implying B ∞ = 0, and the disease free equilibrium E0 is globally stable. Theorem 3.3. If b + γ + sM µM (N3β + η) > δ and R0 < 1, the uninfected macrophage population Mu should satisfy Mu ≥ Mmaxu to completely eliminate TB infection. 46 W. ZHANG, F. FRASCOLI, AND J. M. HEFFERNAN Symbol Description (Unites) Value (Range) Source sM recruitment rate of Mu (1/ml day) 5000 (0.33, 33) [13] [7] [4] sT recruitment rate of T (1/ml day) 6.6 (3300, 7000) [13] [7] [4] µM death rate of Mu (1/day) 0.01 (0.01, 0.011) [13] [7] [4] b loss rate of Mi (1/day) 0.11 (0.05, 0.5) [13] [7] [4] µT death rate of T (1/day) 0.33 (0.05, 0.33) [13] [7] [4] β infection rate by B (1/day) 2 × 10−7 (10−8, 10−5) [13] [7] [4] η bacteria killing rate by Mu rate (1/ml day) 1.25 × 10−8 (1.25 × 10−9, 1.25 × 10−7) [13] [7] [4] γ cell-mediated immunity rate (1/day) 0.5 (0.1, 2) [13] [7] [4] δ proliferation rate of B (1/day) 5 × 10−4 (0, 0.26) [13] [7] [4] cM expansion rate of T induce by Mi (1/day) 10 −3 (10−8, 1) Estimated cB expansion rate of T induce by B (1/day) 5 × 10−3 (10−8, 1) Estimated eM saturating factor of T expansion related to Mi 10 −4 (10−6, 10−2) Estimated eB saturating factor of T expansion related to B 10 −4 (10−6, 10−2) Estimated c half-saturation ratio for Mi lysis (T/Mi) 3 (0.3, 30) Estimated K carrying capacity of B (1/ml) 10−8 (106, 1010) Estimated N1 max MOI of Mi (B/Mi) 50 (50, 100) [13] [7] [4] N2 max No. of B released by apoptosis (T/Mi) 20 (20, 30) [13] [7] [4] N3 N3 = N1/2 (B/Mi) 25 (25, 50) [13] [7] [4] Mu uninfected macrophages Mi infected macrophages B extra and intra-cellular bacteria T CD4 T-cells Table 1. Parameter Symbol, Descriptions, Values, and Sources [3] The occurrence of a backward bifurcation destabilizes the globally stable disease free equilibrium E0 under the condition R0 < 1 and an extra condition to regain stability is needed, i.e. b + γ + sM µM (N3β + η) > δ, as shown in Theorem 3.3. In the next section, we verify the existence of a backward bifurcation computationally and investigate the associated dynamical behaviors by numerical simula- tions. 4. Bifurcation analysis and numerical simulations Consider the n-dimensional nonlinear system with m parameter values dx dt = f(x,p), x ∈ Rn, p ∈ Rm, f : Rn+m → Rn. (4.1) The equilibrium solutions xe = xe(p) are derived from the equilibrium condition f(xe(p),p) = 0, x ∈ Rn, p ∈ Rm. (4.2) The local stability of the equilibrium points xe(p) is determined by the eigenvalues of the Jacobian J(p) = [∂fi(xe(p),p)/∂xj], which are the roots of the corresponding characteristic polynomial equation Pn(λ) = det[λI −J(p)] = λ + a1(p)λn−1 + a2(p)λn−2 + · · · + an−1(p)λ + an(p). (4.3) The necessary and sufficient conditions for zero-eigenvalue bifurcation (zero-singularity) are given in Ref. [15]. ANALYSIS OF SOLUTIONS AND DISEASE PROGRESSIONS FOR A WITHIN-HOST TUBERCULOSIS MODEL 47 (A) (B) B b bT1 bT4 b=0.035 b=0.0298 Time (days) 0 200 400 600 800 1000 1200 0 2 4 6 8 10 12 14 16 B s iz e ( m m -3 ) b= 0.035, IC h b= 0.035, IC l b = 0.0298, ICh b = 0.0298, IC l 1 3 5 7 9 B b 4 1 3 5 7 9 B b bT3 bT2 0 0.5 1 1.5 2 2.5 3 Time (days) 10 4 0 2 4 6 8 10 12 B s iz e ( m m -3 ) 10 7 b=0.31 b=0.288 b=0.288 b=0.4 b=0.288 Figure 1. Bifurcation diagrams of model (2.1) with B vs b and simulations. E0 and E1 are in green and red curves. Four zero-eigenvalue bifurcation are denoted as black points as bT1, bT2, bT3, and bT4. Simulations are carried out for five fixed b values as b = 0.0298, b = 0.035, b = 0.288, b = 0.31, and b = 0.4. The first three b values shows bistability. The last three b values show different progression speed for the bacterial population B. Theorem 4.1. The necessary and sufficient conditions for system (4.1) to have a k-zero singularity at a fixed point (equilibrium), x = xe(p), of the system are given by an(p) = an−1(p) = · · · = an+1−k(p) = 0, (4.4) which ai(p)’s are the coefficients of the characteristic polynomial (4.3). Further, if the remaining co- efficients a1, a2, . . . an−k still obey the Hurwitz conditions for order n − k, then all the remaining eigenvalues of the Jacobian have negative real parts. Based on the results of the uncertainty and sensitivity analysis in Ref. [3], the model is significantly affected by the change of macrophage loss rate b, the infection rate β, cell-mediated immunity rate γ, and bacterial killing rate η. We thus choose the macrophage loss rate b as a bifurcation parameter to verify the analytical result for the backward bifurcation discussed in Theorem 3.1. The other parameter values are fixed and shown in Table 1. Using Theorem 4.1, we numerically find four zero-eigenvalue bifurcation critical points at bT1 = 0.0295, with ET1 = (497833, 122, 8.7, 40), bT2 = 0.2993, with ET2 = (21089, 2664, 45417, 6952168), bT3 = 0.1363, with ET3 = (20, 3055, 49998972, 7575684467), and bT4 = 0.3000036, which yields ET4 = (500000, 0, 0, 20). The summarized bifurcation, equilibria, and their stability are shown in Figure 1. Two values, i.e. b = 0.0298, b = 0.035, are chosen near bT1 and time series show that bistability occur for both b values, with solutions landing onto different equilibria depending on their initial conditions (see panel (A)). This is of interest because it means that the disease can die out or persist to latency depending on the initial infection status. Interestingly, the progression to latency shows different dynamics and lasts different periods of time for the chosen b values. The red and yellow curves take different time to stabilize at their latency levels. 48 W. ZHANG, F. FRASCOLI, AND J. M. HEFFERNAN We then take three b values, i.e. b = 0.288, b = 0.31, and b = 0.4, close to bT2 (see panel (B)). Again, bistability occurs when b = 0.288 on the left of bT2. There is an obvious difference in the speed of disease progression for the three different b values, as shown by the curves in the inset of Figure 1(B). These examples of fast and slow disease progression dynamics seem to confirm the numerical findings in Ref. [3] and are the object of current investigation. 5. Conclusion In this paper, we analyze a four-dimensional within-host model (2.1) for tuberculosis infection, which has been previously proposed and studied numerically in Ref. [3]. We carry out analyses for the well-posedness and boundedness for solutions, existence of the disease free and infected equilibriums and local and global stability analysis. A bifurcation analysis for the disease free equilibrium is also conducted, and a numerical continuation for the infected equilibrium shows when a backward bifurcation occurs. Numerical simulations finally show how fast and slow disease progressions take place close to the bifurcation, with examples of bistability behaviour. This is important because different initial infections can lead to different disease progressions, with considerable differences among latency times. An in-depth analysis of the bifurcation scenario of this model is currently under progress, with the aim of characterising the different, possible behaviours towards infection that TB shows. References [1] M. W. Borgdorff and D. Van Soolingen, The re-emergence of tuberculosis: what have we learnt from molecular epidemiology? Clin. Microbiol. Infect., 9(2013), 889-901. [2] C. Castillo-Chavez and B. Song, Dynamical models of tuberculosis and their applications. Math. Biosci. Eng., 1(2004), 361-404. [3] Y. Du, J. Wu and J. M. Heffernan, A simple in-host model for Mycobacterium tuberculosis that captures all infection outcomes. Math. Popul. Stud., 24(2017), 37-63. [4] D. Gammack, S. Ganguli, S. Marino, J. Segovia-Juarez and D. E. Kirschner, Understanding the immune response in tuberculosis using different mathematical models and biological scales. Multiscale Model Sim, 3(2005), 312-345. [5] W. M. Hirsch, H. Hanisch and J. P. Gabriel, Differential equation models of some parasitic infections: methods for the study of asymptotic behavior. Commun. Pur. Appl. Math., 38(1985), 733-753. [6] P. L. Lin and J. L. Flynn, Understanding latent tuberculosis: a moving target. J. Immunol., 185(2010),15-22. [7] S. Marino and D. E. Kirschner, The human immune response to Mycobacterium tuberculosis in lung and lymph node. J. Theor. Biol., 227(2004), 463-486. [8] J. DH Porter and K. PWJ McAdam, The re-emergence of tuberculosis. Annu. Rev. Publ. Health, 15(1994), 303-323. [9] A. Sakula, Centenary of the discovery of the tubercle bacillus. The Lancet, 319(8274)(1982), 750. [10] D. Sud, C. Bigbee, J. L. Flynn and D. E. Kirschner, Contribution of CD8+ T cells to control of Mycobacterium tuberculosis infection. J. Immunol., 176(2006), 4296-4314. [11] P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compart- mental models of disease transmission. Math. Biosci., 180(2002), 29-48. [12] WHO, Golbal tuberculosis report 2019. WHO, 2019. [13] J. E. Wigginton and D. Kirschner, A model to predict cell-mediated immune regulatory mechanisms during human infection with Mycobacterium tuberculosis. J. Immunol., 166 (2001), 951-1967. [14] J. Yang, T. Kuniya, F. Xu and Y. Chen, Evaluation of the tuberculosis transmission of drug-resistant strains in mainland china. J. Biol. Syst., 26(2018), 533-552. [15] P. Yu, Closed-form conditions of bifurcation points for general differential equations. Int. J. Bifurcat. Chaos, 15(2005), 1467-1483. [16] K. Zaman, Tuberculosis: a global health problem. J. Health Popul. Nutr., 28 (2) (2010), 111. ANALYSIS OF SOLUTIONS AND DISEASE PROGRESSIONS FOR A WITHIN-HOST TUBERCULOSIS MODEL 49 Corresponding author, Department of Mathematics and Statistics, Texas Tech University, Broadway and Boston, Lubbock, TX 79409-1042. E-mail address: wenjing.zhang@ttu.edu Department of Mathematics, Faculty of Science, Engineering and Technology, Swinburne University of Technology, John St, 3122, Hawthorn, VIC, Australia. E-mail address: ffrascoli@swin.edu.au Department of Mathematics and Statistics, Centre for Disease Modelling, York University. 4700 Keele St, Toronto, ON, Canada, M3J 1P3. E-mail address: jmheffer@yorku.ca