An investigation of tuberculosis progression revealing the role of macrophages apoptosis via sensitivity and bifurcation analysis.

Mycobacterium tuberculosis infection features various disease outcomes: clearance, latency, active disease, and latent tuberculosis infection (LTBI) reactivation. Identifying the decisive factors for disease outcomes and progression is crucial to elucidate the macrophages-tuberculosis interaction and provide insights into therapeutic strategies. To achieve this goal, we first model the disease progression as a dynamical shift among different disease outcomes, which are characterized by various steady states of bacterial concentration. The causal mechanisms of steady-state transitions can be the occurrence of transcritical and saddle-node bifurcations, which are induced by slowly changing parameters. Transcritical bifurcation, occurring when the basic reproduction number equals to one, determines whether the infection clears or spreads. Saddle-node bifurcation is the key mechanism to create and destroy steady states. Based on these two steady-state transition mechanisms, we carry out two sample-based sensitivity analyses on transcritical bifurcation conditions and saddle-node bifurcation conditions. The sensitivity analysis results suggest that the macrophage apoptosis rate is the most significant factor affecting the transition in disease outcomes. This result agrees with the discovery that the programmed cell death (apoptosis) plays a unique role in the complex microorganism-host interplay. Sensitivity analysis narrows down the parameters of interest, but cannot answer how these parameters influence the model outcomes. To do this, we employ bifurcation analysis and numerical simulation to unfold various disease outcomes induced by the variation of macrophage apoptosis rate. Our findings support the hypothesis that the regulation mechanism of macrophage apoptosis affects the host immunity against tuberculosis infection and tuberculosis virulence. Moreover, our mathematical results suggest that new treatments and/or vaccines that regulate macrophage apoptosis in combination with weakening bacillary viability and/or promoting adaptive immunity could have therapeutic value.

Introduction

Tuberculosis (TB), a disease caused by Mycobacterium tuberculosis (Mtb), remains the deadliest infectious disease in the world (Sud et al. 2006; Lin and Flynn 2010). According to a CDC (Centers for Disease Control and Prevention) report (WHO 2019), in the United States, the reactivation of latent Tuberculosis infection was a major cause for TB incidence in 2019. A detailed understanding of the most significant driving factors for TB disease progression is crucial for the improvement of the treatment strategies to help host fight against the invader pathogen.

To serve this motivation, we adopt an established in-host tuberculosis infection model that generates various disease outcomes. Since transcritical bifurcations determine the fate of infection and saddle-node bifurcations create multiple infected equilibrium solutions, we carry out a parameter sensitivity analysis on the conditions for transcritical and saddle-node bifurcations. We then further investigate model behavior through bifurcation analysis to examine how the macrophage apoptosis rate and other statistically significant parameters impact disease progression and examine how the mathematical results relate to the complexities of the disease dynamics. Finally, we utilize a 2-dimensional bifurcation analysis to explore how combination therapies can enhance macrophage apoptosis modulation therapy to inhibit bacillary viability and promote adaptive immunity.

Before we present our methodology and results, we provide a brief introduction of tuberculosis, including the various disease outcomes, in-host Mtb life cycles, and therapeutic strategies, in Sect. 1.1 and discuss the previous in-host Mtb mathematical modeling literature in Sect. 1.2.

TB epidemiology and immunology background

Tuberculosis (TB), a disease caused by pathogenic bacterium species Mycobacterium tuberculosis (Mtb), has been spreading among humans for thousands of years (Donoghue et al. 2004). Despite the availability of effective TB treatments, TB continues claiming 1.5 million lives per year by infecting one quarter of the global population (Sud et al. 2006; Lin and Flynn 2010). A characteristic feature for Tuberculosis infection is that infected individuals experience various disease outcomes, including clearance, latent infection (LTBI), and active disease. Approximately 10% of infected individuals clear the infection and another 10% progress to primary disease, while the majority of infected individuals become latently infected with no symptoms. Individuals with LTBI carry the risk of developing TB later in life. Slow TB progression through LTBI reactivation contributes approximately of cases globally and threatens a TB re-emergence (Cooper 2009; Rich et al. 1951). A better understanding on TB progression, especially LTBI reactivation, can help to improve the efficacy of therapeutic strategies, and enhance both the innate and adaptive immune responses.

The life cycle of tuberculosis begins when aerosol droplets containing Mtb are inhaled by an uninfected individual. After inhalation, Mtb usually reach the alveoli of the lung, and are preferentially ingested by alveolar macrophages (Warner and Mizrahi 2007), then launch proliferation (Canetti 1955). This makes alveolar macrophages the main target cells for Mtb invasion (Gammack et al. 2005; Gideon and Flynn 2011). Even though activated macrophages can initialize phagocytosis, Mtb can also set up various strategies against macrophage killing and replicate within the host macrophages. Moreover, cytokine environment is also modified to favor the intracellular Mtb survival (BoseDasgupta and Pieters 2014; Orme et al. 2015). Innate immune responses are the crucial factor for the TB outcome (Orme et al. 2015). While the clearance of the phagocytized bacteria is also determined by the adaptive immunity, which is activated by innate immune signals (Torrado and Cooper 2013). These infected macrophages release chemotactic signals to recruit and activate migrating lymphocytes to the site of infection (Gammack et al. 2005; Gideon and Flynn 2011). The progressive immune responses lead to an aggregation of lymphocytes around mycobacteria-harboring macrophages, which is called a granuloma. Granulomas are able to isolate infected macrophages from the rest of the lung to limit the bacillary growth and spread of Mtb from intracellular Mtb replication, and so serve as an immune micro-environment to enhance the T cells and macrophages interactions (Russell et al. 2010; Gideon and Flynn 2011). On the other hand, the granuloma concealed Mtb switch to a non-replicating state, which is a typical feature to enable the bacillus persistence during the long TB latency period (Barry et al. 2009). Mtb has evolved new strategy to escape granulomas: leak into the lung and expectorate bacilli (Ehlers and Schaible 2013). The granuloma then serves as a battle ground between mycobacteria and host immune cells. The host-pathogen dynamics decide granuloma to promote or inhibit the infection and then determine the disease outcomes (Ulrichs and Kaufmann 2006; Ehlers and Schaible 2013). An imbalanced interplay between Mtb and immune cells result in the failure of controlling infection, then lead to the LTBI reactivation (O’Garra et al. 2013).

Today, TB disease can be cured with antibiotics. For the newly diagnosed pulmonary TB, four first-line anti-TB drugs isoniazid (INH), rifampin (RIF), ethambutol (EMB), and pyrazinamide (PZA) form the essence of treatment regimens (Organization 2017). The major problem of the current TB drug regimen is its duration, which requires 6 to 9 months. Incorrect TB drug treatment (the wrong length of duration or dose) could lead the remaining alive Mtb resistant to those drugs. These patients are more likely to develop drug-resistant TB, which is more difficult to treat. The goal of shortening the long treatment duration can advance treatment effectiveness by improving patients’ compliance and reducing the chance of the drug resistance and disease relapse. Mathematical models may be useful tools for, first, elucidating in-host TB progression mechanisms, and then also for suggesting crucial factors for advancing therapeutic strategies.

Mathematical modeling background

Numerous TB in-host modeling studies have been done to identify the driving factor for different TB infection outcomes. The TB host-immune modeling began with models considering HIV-1 and Mtb coinfection by Kirschner (1999) and her collaborators Bauer et al. (2008). HIV infection weakens the immune system, and make patients vulnerable to opportunistic infections, such as TB. Untreated latent TB and HIV coinfection significantly increases the possibility of developing TB disease. The study of TB infection is carried out through ODE models to capture the cellular and cytokine network (Antia et al. 1996; Abdelrazec et al. 2016; Wigginton and Kirschner 2001; Magombedze et al. 2006) and CD8+ T cells for host infection control (Sud et al. 2006; Magombedze et al. 2006), PDE models have been developed for cytokine and antibiotics molecules diffusion (Lauffenburger and Linderman 1993), agent-based models for spatial and temporal features of granuloma (Segovia-Juarez et al. 2004; Ray et al. 2009) and multiple scale models for investigating the granulomas dynamics in a smaller scale (such as a lung) and a larger scale (such as a host body) (Gong et al. 2015; Marino and Kirschner 2016). However, little work has been done to investigate the determining factors for LTBI reactivation. One of the challenges is the computational cost due to the large number of considered compartments and parameters, combined with complex nonlinear terms. In this paper, we only consider the key factors involved in the in-host Mtb dynamics, which include the infection interaction and innate immune response between macrophages and bacteria, the adaptive immune response initiated by T cells, and bacterial growth and elimination processes.

A four-dimensional within-host TB model was developed in previous modeling work by Du et al. (2017). This 4-dimensional model, presented in (1), incorporates the uninfected macrophages (the Mtb target cells), infected macrophages, Mtb bacteria and CD4 T cells. The model numerically captures various TB infection outcomes, which are clearance, latency, and active disease. In an analytical study of the same model (Zhang et al. 2020), we showed the existence of four coexisting equilibria that reflect the different disease outcomes. We also determined conditions for forward and backward bifurcations (Zhang et al. 2020). However, neither (Du et al. 2017) or (Zhang et al. 2020) has analyzed the driving factors behind the different outcomes of disease.

Overview

We provide a brief immunology and mathematical modeling background of tuberculosis infection in Sect. 1. In Sect. 2, we introduce an established in-host tuberculosis model, then identify parameters statistically significantly influencing the early stage TB infection through a sample-based sensitivity analysis on the basic reproduction number. Using identified parameters from the first sensitivity analysis, we employ a second sensitivity analysis on the single-zero eigenvalue bifurcation condition and identify that macrophage apoptosis rate is the statistically significantly parameter for disease progression, especially LTBI reactivation. In Sect. 3, using bifurcation analysis and numerical simulations, various disease outcomes are fully unfolded under the variation of the macrophage apoptosis rate. In Sect. 4, we describe how the mathematical results from the previous section reveal the complicated roles of macrophage apoptosis in the complex Mtb-macrophage dynamics. In Sect. 5, we further investigate the disease outcomes under the influence of macrophage apoptosis and one of the other four parameters to give insight into the development of new therapeutic strategies. Our findings suggest that macrophage apoptosis modulation therapy can be enhanced with combination therapies, which inhibit bacillary viability and promote adaptive immunity.

Investigation of the driving factors for disease progression through sensitivity analysis

Tuberculosis progression is characterized by reactivation of LTBI, which occurs when a weakened immune system is no longer capable of controlling the latent stage Mtb. The awoken bacteria become active, defeat the immune defense, and lead to the infected individual to active TB disease. In other words, the LTBI reactivation can be described as the Mtb bacterial level shifted from a lower infection steady state to a higher infection steady state. Our target in this section is to identify the important factors determining the shifting of Mtb bacterial level by evaluating the influence/sensitivity level of each involving mechanism applied on the shift.

Model

The TB model with various disease outcomes (Du et al. 2017) is written as follows:Here, uninfected macrophages are generated and die naturally at constant rates and . Moreover, they become infected through the engulfment of Mtb at a rate of . The internalized Mtb reproduction can kill infected macrophages at a rate b and produce intracellular Mtb bacteria per cell, on average. Infected macrophages can also be eliminated through T cell-mediated cytotoxicity. Both CD4 and CD8 T cells account for the cytotoxic process. The rate of the cytotoxic effect, which T lymphocytes (E) act on Mtb-stimulated monocytes (T), depends on the E:T ratio, (Lewinsohn et al. 1998; Wigginton and Kirschner 2001). This ratio has a maximum rate and saturating factor c. During the cytotoxic process, infected macrophages die and release on average intracellular Mtb bacteria to the extracellular compartment B. Extracellular bacteria can divide, modeled as a logistic term . The infection in macrophages can lead to a loss of extracellular Mtb. This is modeled using the term . On the other hand, uninfected macrophages phogocytize Mtb bacteria, then initiate phagosome maturation, which causes microbicidal and degrading effects on phogocytized pathogens. The Mtb killing rate by this process is denoted as (Kinchen and Ravichandran 2008; Upadhyay et al. 2018). Host adaptive immune responses are signaled by CD4 T cells, which are assumed to have constant production and death rates and . CD4 T cells can proliferate in the presence of Mtb bacteria () or/and infected macrophages (). More detailed model assumptions are explained in Du et al. (2017) and parameters values are stored in Table 1.

Table 1

Parameter descriptions and values (Du et al. 2017; Zhang et al. 2020)

Symbol	Definition	Unites	Mean value (Range)	References
s_M	Uninfected macrophages’ recruitment rate	ml^-1 day^-1	5000(3300,7000)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
s_T	T-cell recruitment rate	ml^-1 day^-1	6.6(0.33,33)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
μ_M	Uninfected macrophages’ death rate	day^-1	0.01(0.01,0.011)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
b	Infected macrophages’ clearance rate	day^-1	0.11(0.05,0.5)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
μ_T	T-cell death rate	day^-1	0.33(0.05,0.33)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
β	Mtb infection rate	day^-1	5×10^-6 (10^-8,10^-5)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
η	Mtb death rate by uninfected macrophages’ killing	ml^-1 day^-1	1.25×10^-9 (1.25×10^-9,1.25×10^-7)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
γ	T-cell mediated immunity rate	day^-1	1.5(0.1,2)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
δ	Mtb proliferation rate	day^-1	5×10^-4 (0,0.26)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
c_M	T-cell proliferation rate induced by infected macrophages	day^-1	10^-3(10^-8,1)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
c_B	T-cell proliferation rate induced by Mtb	day^-1	5×10^-3 (10^-8,1)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
e_M	Saturating factor for T-cell proliferation rate induced by infected macrophages	−	10^-4(10^-6,10^-2)	Estimated
e_B	Saturating factor for T-cell proliferation rate induced by Mtb	−	10^-4(10^-6,10^-2)	Estimated
c	T cellsuninfected macrophages half-saturation ratio for uninfected macrophages lysis	TM_i	3(0.3,30)	Estimated
K	Mtb carrying capacity	ml^-1	10^8(10^6,10^10)	Estimated
N_1	Average number of Mtb released by one macrophage	BM_i	50(50,100)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
N_2	Maximum number of Mtb released by one macrophage lysis	TM_i	20(20,30)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)
N_3	Average number of Mtb carried by one macrophage	BM_i	25(25,50)	Wigginton and Kirschner (2001), Marino and Kirschner (2004), Gammack et al. (2005)

Basic model properties

A detailed analytical investigation on the property of model (1) is carried out in Zhang et al. (2020). For convenience, we change the notation asThe solutions of model (1) are well-posed and bounded and include a disease free equilibrium () and infected equilibriums () as follows (Zhang et al. 2020)Furthermore, equilibrium solutions and intersect and exchange local stability at a transcritical bifurcation asMoreover, at the transcritical bifurcation the basic reproduction number equals to unit (Du et al. 2017), i.e. , where

Sensitivity analysis of for the early stage of TB infection

After the exposure of TB disease, the critical value to determine the fate of the infection is the basic reproduction number . Starting from one infectious pathogen, the infection will be able to spread to the whole body if , but will be eliminated if . Because is not a constant number for any given disease, in the case of in-host TB dynamics, is affected by in-host immunological factors and Mtb behaviors. We first investigate in (5) for all model parameters in the experimental and/or estimated ranges from Table 1. We then choose in (5) as the first sensitivity/output function of the model parameters (or input factors) to identify which parameters significantly affect the uncertainty on disease spread and elimination in the early stage of the TB infection. The analytical formulation (5) indicates that is linear on the parameters , , , and , along with the ratio . also depends on the term .

To show the uncertainty on due to the uncertainties in model parameters, we perform uncertainty analysis. Because a priori knowledge of the Mtb-host relationship between the model parameter and is unavailable, we choose the most popular sampling-based approach, i.e. Latin Hypercube sampling (LHS) (McKay et al. 1979). LHS belongs to the Monte Carlo class of sampling procedures for the propagation of uncertainty. LHS performs a full stratification over each sampled-parameter range to allow every parameter get equally sampled, and returns unbiased estimates for the sample mean and distribution function over each sampled-parameter range. Moreover, when necessary, a log scale can be applied on the parameter range for LHS to prevent under-sampling (McKay et al. 1979; Marino et al. 2008). Furthermore, LHS nearly eliminates confounding in estimating parameter effects by ensuring that exactly one parameter is sampled from each stratum in a given sample. That is, if there are 3 parameters, then each parameter range is partitioned into 3 strata and one parameter is sampled from the lowest stratum, another from the middle stratum, and the last from the highest stratum. This is repeated many times to ensure that all parameters are sampled from each of their sampling strata. This results in the main effects of each parameter being confounded only with high-order interactions of other parameters, which are classically considered negligible in experimental design. This also makes LHS a computationally inexpensive approach, especially compared to full factorial sampling designs, and allows LHS to perform robustly for even a relatively small sample size, i.e. sample size to 200, which is used in Iman and Helton (1988); Kleijnen and Helton (1999); Helton and Davis (2003). As such, this sampling design facilitates a measurement of parameter’s corresponding impact on the value of the sensitivity functions in (5) and in (12).

Due to the multi-variable property of sensitivity functions, standard correlation coefficients (CCs) are not sufficient for measuring the influence of parameter values on sensitivity functions because CCs are restricted to quantify relationships between just two variables and, as such, cannot account for the impact of the other parameters on the variation in sensitivity functions. Instead, we can calculate partial correlation coefficients (PCCs), which quantify the relationship between each parameter on the portion of the variation in sensitivity functions that is unexplained by the variation in the values of the other parameters and thusly, provide a more appropriate measure of influence. However, PCCs only reflect the linear relationships between the output of sensitivity functions and the input variables/parameters. To capture the nonlinear relationships between sensitivity functions and the model parameters, we can instead apply the rank transformation to all quantities (Conover and Iman 1981; Saltelli and Sobol 1995). Calculating the correlations between these provided us with the Spearman correlation coefficients, which measure the strength and direction of monotonic relationships. PCCs for the ranked quantities yields partial rank correlation coefficients (PRCC), which provide a robust sensitivity measure for quantifying potentially nonlinear, but monotonic, relationships between the output and each of the model parameters in the presence of the other parameters (Marino et al. 2008; Saltelli and Annoni 2011).

Using the MATLAB LHS and PRCC code provided in Marino et al. (2008), we run an LHS uncertainty analysis on a sample of 1000 parameter sets, and set the baseline value/mean, lower bound, and upper bound for each parameter in its range from Table 1. Scatter plots showing the relationships between residuals of the output and input parameters, along with their PRCCs and associated p-values, are shown in Figure 1 in Appendix (see supplementary material). From these figures, it is readily apparent that all of the pairs with near-zero PRCCs display no evidence of non-monotonic relationships, so PRCCs appropriately describe the relationships. A summary of the PRCC values is shown as a bar graph in Fig. 1. Following traditional use, p-values below 0.05 designate pairs of parameters whose partial rank correlations are statistically significantly different than 0. The negative or positive sign of PRCC indicate the negative or positive correlation between and the specific parameter.

Fig. 1

PRCC values for under the variations of all model parameters. n.s. denotes statistically not significant, that is p value

The PRCC between and each parameters in Fig. 1a indicate that the expected number of bacteria, , generated by one single bacterium after exposure to Mtb hasThere is an identifiablility issue with the parameters , and , because they only appear in (5) together as the term . While each of these parameters should be expected to have a non-zero rank correlation with ,since partial correlation accounts for the effects of the other parameters, only one of these three can have a high PRCC value. In this case, the effect appeared present for rather than the other two. Since the uninfected macrophages level in a healthy host is relatively constant, the Mtb proliferation rate is chosen to represent the significantly positive influence on . Also note that the other non-influential parameters are , , c, K, , , , , which do not show in the analytical formula of in (5).

a significantly positive relation with , , and ;

a significantly negative relation with because the infection process involves the loss of extracellular Mtb by macrophages phagocytosis;

a significantly negative relation with ;

a significant relation with the whole term , indicating that more intracellular Mtb are released from apoptosis than necrosis (Lee et al. 2009);

a positive relation with the whole term .

Sensitivity analysis of saddle-node bifurcations for LTBI reactivation

LTBI reactivation is characterized by an in-host Mtb outbreak. The Mtb level shifts from a low-level latent infection stage to a high-level active disease stage by the slowly weakening immune system, which may be caused by AIDS, chemotherapy, aging, and diabetes. In our mathematical model, the LTBI reactivation can be captured by the bacteria level traveling from a low steady state to a high steady state due to slowly changing parameter values. The prerequisite condition is the existence of co-existing steady states. The fundamental mechanism for the generation and destruction of steady states is the occurrence of saddle-node bifurcation. Therefore, the parameter influence on the shifting of the infection steady states is replaced by the parameter influence on the occurrence of saddle-node bifurcation. The sensitivity analysis of the LTBI reactivation is taken over by the sensitivity analysis of the saddle-node bifurcation conditions.

In general, we consider an n-dimensional nonlinear systemwith m parameter values and equilibrium solutions derived from the equilibrium conditionThe local stability of the equilibrium points is determined by the eigenvalues of the Jacobian , which are the roots of the corresponding characteristic polynomial equationThe necessary and sufficient conditions for zero singularities are given in Yu (2005).

Theorem 1

Yu (2005) The necessary and sufficient conditions for system (6) to have a k-zero singularity at a fixed point (equilibrium), , of the system are given bywhich ’s are the coefficients of the characteristic polynomial (8). Further, if the remaining coefficients still obey the Hurwitz conditions for order , then all the remaining eigenvalues of the Jacobian have negative real parts.

For the current case, we only need to study single zero bifurcations, that is . The corresponding necessary condition becomes , since the dimension of model (1) is 4. Single zero bifurcations include saddle-node/turning bifurcation, which creates or destroys steady states, transcritical bifurcation (at ), which serves as the intersection between infection-free and infected steady states (Zhang et al. 2020), and pitchfork bifurcation, which is not common in within-host infectious disease model, due to the symmetric property.

The sensitivity function on LTBI reaction for model (1) is represented by the sensitivity function of the necessary condition of single zero bifurcation (denoted by ) in terms of the parameters (denoted by ) significantly affecting the transcritical bifurcation condition (). The sensitivity analysis in Sect. 2.3 identified the significant parameters affecting to beWe fix the parameters, which are not significantly affecting or the fate of disease as the corresponding mean values in Table 1. There parameters are denoted asDue to the uncertainties of the input parameters , we run Latin Hypercube sampling for 10, 000 sampled parameter set. The sensitivity function/output is set aswhere denotes the infected steady state and satisfies in (3). Moreover, the non-influential parameters in terms of disease outcomes are fixed. Scatter plots showing the relationships between the disease-outcome influential parameters and the sensitivity function in (12) are shown in Fig. 2. As in the preceding sensitivity analysis, these plots show no evidence of non-monotonicity, which again suggests that PRCC scores are appropriate to use. The PRCC scores for these are shown in bar graph form in Fig. 2. The sensitivity analysis identifies thatTo further investigate how the identified significant parameter affects disease progression and what disease outcome will show up, we set b as a bifurcation parameter of interest for a one-parameter bifurcation analysis in the next section. More two-dimensional bifurcation diagrams are provided in the later section on the identified influential parameters to provide insight on TB treatment strategies.

Fig. 2

PRCC values for single zero-eigenvalue bifurcations in (12) on the influential parameters for disease outcomes

b has the largest magnitude PRCC value with the sensitivity function ;

is non-influential to the LTBI reactivation;

the other parameters in are statistically significant but have a moderate absolute PRCC value with the sensitivity function/output .

Bifurcation analysis of the macrophages programmed cell death rate

Investigation of the single-zero eigenvalue bifurcation

The single-zero eigenvalue bifurcations are derived according to Theorem 1. Leaving the most influential parameter b as the bifurcation parameter, fixing the other influential parameters in and non-influential parameters in as the corresponding mean value in Table 1, we rewrite model (1) using the new notation in (2) asSuppose the steady state/equilibrium is given in the formwhich could represent multiple co-existing equilibrium solutions. Then the Jacobian matrix of (13) at the infected equilibrium in (3) iswhere satisfies the equilibrium conditions in (3). We denote the fourth-degree characteristic polynomial in terms of L at the infected equilibrium in (3) aswhere satisfies the equilibrium condition in (3) . Evaluating this with the parameter values in Table 1, the formula of and are shown in Appendix (see supplementary material) Equations (1),(3),(4),(5), and (6). We assume We obtain four single-zero bifurcation critical points: Single-zero bifurcation include saddle-node/turning bifurcation, transcritical bifurcation, and pitchfork bifurcation, through the pitchfork bifurcation is not common in infectious disease models due to the symmetric pattern. To identify the type of the single-zero bifurcation and the local behavior of the corresponding vector field, We apply Theorem 3.4.1 in Guckenheimer and Holmes (1990). That is, for the system (13), there are three hypotheses SN1, where w and v are the left and right eigenvectors corresponding to the 0 eigenvalue of the matrix . Moreover the transversality and degeneracy conditions are defined respectively asIf system (13) satisfies three hypotheses SN1, SN2, and SN3, there are no equilibriums that exist near when () and two equilibria near when (). It implies a saddle-node bifurcation. If we replace SN2 by SN2’, then a transcritical bifurcation occurs.

the Jacobian matrix has a simple zero eigenvalue, at , if (1) , (2) , (3) , and (4) or . Then and is called critical point (Yu 2005).

day, cells/ml.

day, cells/ml.

day, cells/ml.

day, cells/ml.

,

,

For the critical point , the left and right eigenvectors corresponding to the simple zero eigenvalue of are chosen respectively aswhere . The transversality and degeneracy conditions are satisfied sinceThe signs of the transversality and degeneracy conditions indicate that no equilibrium point near when , but two equilibria are near when . Choosing , we test the stability of the two bifurcation equilibria: (associated with four negative real-part eigenvalues , , , and ) and (associated with one positive and three negative eigenvalues , , , and 0.0059). Given these results, a saddle-node bifurcation occurs at shown in Fig. 3a,b as LP. The bifurcating upper branch equilibrium is stable and the lower branch is unstable.

Fig. 3

Bifurcation diagram of Mtb concentration B ( ) vs macrophage apoptosis rate b. Infected () equilibriums are plotted according to Equation 1 in Appendix (see supplementary material) as in blue lines (see (a) and (b)). The disease-free () is in red. The solid and broken line styles distinguish stable and unstable equilibriums. Oscillating bacterial load are due to two supercritical Hopf bifurcations. Two oscillation regions are separated because periodic solutions merge with the lowest branch of (see (b)) and generate blow-up periods (see (c)) (color figure online)

For the critical point , the left and right eigenvectors corresponding to the simple zero eigenvalue of are chosen respectively aswhere . The transversality and degeneracy conditions are satisfied sinceThe other three negative eigenvalues associated with are , , and . Since both the transversality and degeneracy conditions are positive, there is no equilibrium near when , but there are two equilibria near when . To test the stability of the bifurcating equilibria, we choose and obtain an unstable upper equilibrium (associated with four eigenvalues , , , and 0.00048) and a stable lower equilibrium (associated with four eigenvalues , , , and ). There also exist an unstable positive equilibrium (associated with four eigenvalues , , , and 0.0094) and a stable positive equilibrium (associated with four eigenvalues , , , and ). The dynamics are summarized in Fig. 3a near the point LP.

For the critical point , The left and right eigenvectors corresponding to the simple zero eigenvalue of are chosen respectively aswhere . The transversality and degeneracy conditions are satisfied sinceThe rest of the three negative eigenvalues associated with are , , and ). The signs of the transversality and degeneracy conditions indicate no equilibrium near when , and two equilibria near when . Choosing we test the stability of the bifurcating branches and observe that there are two equilibria near , a stable (associated with four eigenvalues , , , and ) and an unstable 4396306898) (associated with four eigenvalues , , , and 0.0002). The other two unstable positive equilibria located further away from are (associated with four eigenvalues , , , and ) and (associated with four eigenvalues , , , and 0.0709). The dynamics are summarized in Fig. 3a near the point LP.

At , the disease free equilibrium and infected equilibrium intersect and exchange stability. The qualitative behavior of the model (1) at is determined by its corresponding center manifold governed by . The formula of the coefficients and are derived in our previous work (Zhang et al. 2020). Using the parameter values in Table 1, we have and . Therefore, the model (1) undergoes a transcritical bifurcation at . Furthermore it shows a backward bifurcation due to the positiveness of the coefficients and in the quadratic terms. The dynamics are summarized in Fig. 3b near the point BP.

Note that the magnitudes of transversality and degeneracy conditions depend on the normalization of the corresponding left and right eigenvalues and (or and and and ). While the signs of transversality and degeneracy conditions (which are related to the opening direction of the saddle node bifurcation) are invariant under the scaling of the left and right eigenvectors. Therefore, we choose the left and right eigenvectors such that their dot product equals to one. Furthermore, the signs of transversality and degeneracy conditions for , , and show that the opening direction of the saddle node bifurcation LP is different with LP and LP (see Fig. 3).

Investigation of the Hopf bifurcations

The necessary and sufficient conditions (Wm 1994; Yu 2005; Yu and Wang 2019) for Hopf bifurcation require Here, , are Hurwitz arguments and the function is the equilibrium condition. The above four criteria yield two Hopf bifurcation points: We further study the post-critical behavior of the Hopf bifurcations to determine the stability of the bifurcating limit cycles. The corresponding Hopf bifurcation normal forms up to third order terms are derived through center manifold and normal form reductions. Due to the high non-linearity and stiffness of the model equations, we employ an established Maple program (Yu 1998) for the model analysis.

,

,

(positive determinant of the corresponding Jacobian matrix),

.

day, cells/ml, with , , and , and . The corresponding eigenvalues are with .

day, cells/ml, with , and , . The corresponding eigenvalues are .

The Hopf bifurcation at

We first apply a parameter-dependent shift of coordinate on the equilibrium condition in Eq. 3, that isSolving the preceding equation obtainsTaking only the linear approximation of , we assume that , where , and are small perturbation parameters. It follows that . We introduce the transformationinto (1) and obtainwhere for are stored in Appendix (see supplementary material) (7)–(10). The corresponding Jacobian matrix of equations evaluated at the equilibrium , takes the Jordan canonical form aswhere . The corresponding Hopf bifurcation normal form up to third order term iswhere and are derived from linear analysis as followsThe first focus value and the coefficient are derived from the Maple program (Yu 1998) taking as input. We obtainThe trivial equilibrium indicates the equilibrium for model (1). While the nontrivial equilibrium exists if and represents the amplitude of the bifurcating limit cycle. SinceThe trivial equilibrium is stable when . When , becomes unstable, then a stable limit cycle merged with an amplitude of .

Employing the similar approach as shown in the previous subsection, we first take the linear approximation of solved from . That is, with ( and are small perturbations). Then we introduce the following transformationinto (1). It yields , for , where are stored in (11)-(14) in the Appendix (see supplementary material). The expression of functions are omitted for brevity. The corresponding Jacobian matrix of equations evaluated at the equilibrium () is in the Jordan canonical form , where . The corresponding Hopf bifurcation normal form up to third order term is written as and , where , , , and . The trivial equilibrium indicating is stable if and unstable if , while the non-trivial equilibrium exists and is stable if . A stable limit cycle is generated from this Hopf bifurcation.

Demonstrating roles of macrophage programmed cell death through numerical simulations

The bifurcation diagrams in Fig. 3 demonstrate the diverse Mtb levels B with regard to the infected macrophage programmed cell death (apoptosis) rate b. Disease-free equilibrium (red line in Fig. 3b) and infected equilibrium (blue curve in Fig. 3a, b) intersect and exchange stability at a transcritical bifurcation (marked as point BP at ). In the first branch from the bottom of the infected equilibrium, is unstable and undergoes two neutral saddle points (marked as NS), then bifurcates to a stable equilibrium solution in the second branch from the bottom at a saddle-node bifurcation (marked as LP at ). The stable loses its stability at a supercritical Hopf bifurcation at (marked as H). The stable bifurcating limit cycles disappear when the periodic solution merges with the lowest branch of (see Fig. 3b) and induces a blow up in the period (see Fig. 3c top). keeps being unstable until reaching to a second supercritical Hopf bifurcation at (marked as H). The corresponding bifurcating stable limit cycles also vanish when the cycles merge with the lowest branch of (see Fig. 3b), and induce a blow up in the period (see Fig. 3c bottom). This leaves a window in the branch at the second from the bottom of (approximately day), where solutions converge to the stable disease-free equilibrium (red line). Bistability is possible if day. The stable branch at the second from the bottom of after H bifurcates an unstable branch (at the second from the top) at another saddle node bifurcation (marked as LP) at . Then gains back stability at the third saddle node bifurcation (marked as LP) at . It yields the top stable branch for the infected equilibrium . Overall, the model (1) admits various disease outcomes (clearance, latent infection (LTBI), and active infection Cadena et al. 2017) through the occurrence of three saddle node bifurcations (LP at , LP at , LP at and BP at ). Interestingly, in the second branch from the bottom, where Mtb level B is in latent stage (), the Mtb level first stabilizes at a low level (), then oscillates , then vanished ( day), and oscillate again ( day), then gain back stability (). These surprisingly complex dynamical behaviors motivate us to investigate the macrophage-Mtb interaction model 1 regarding the macrophage apoptosis against Mtb.

Macrophage programmed cell death in TB infection

Before interpreting the complex macrophage-Mtb interplay shown in the bifurcation diagrams in Fig. 3 and the simulations in Figs. 4 and 5, it may be helpful to review macrophage apoptosis as a defense strategy against Mtb infection. Host immune responses to tuberculosis infection start from the inhalation of air droplets containing Mtb pathogens. These bacilli first manage to enter into resident alveolar macrophages through the phagocytosis process. A phagolysosome is then formed by the fusion of lysosomes with phagosomes containing engulfed Mtb pathogens. In the phagolysosome, internalized bacilli are degraded by reactive oxygen species and the acidification. However, Mtb pathogens evolve to survive from this cell anti-bacterial mechanism and to reproduce inside the macrophage. At this point, infected macrophages prohibit the extracellular Mtb pathogens, but facilitate intracellular bacterial replication at least until the activation of the adaptive immune responses (Sturgill-Koszycki et al. 1994), and provide an intracellular environment in favor of bacterial latency (Kornfeld et al. 1999). To defend against intracellular parasitic Mtb pathogens, a common strategy is to initiate host cell apoptosis. It destroys the protected intracellular environment, which promote intracellular bacterial reproduction and latency (Lee et al. 2009). However, to confront the elimination through host cell apoptosis, Mtb pathogens (especially virulent Mtb strains) evolve to interrupt tumor necrosis factor- (TNF-) signaling (Spira et al. 2003) and promote anti-apoptotic Mcl-1 (Sly et al. 2003). These signals suppress the infected macrophage apoptosis. On the other hand, after maintaining the suitable intracellular environment, virulent strain Mtb pathogens manage to replicate, but still need a mean to spread the newly replicated pathogens. Host cell death is then triggered by heavy load of intracellular bacteria to spread the infection to other host cells. This heavily infected macrophage apoptosis process progresses rapidly and ends up in macrophage necrosis. It is found that the intracellular Mtb pathogens released from this process are not killed by host immune cells and able to replicate as extracellular bacteria (Lee et al. 2006). The macrophage necrosis caused by heavy load of intracellular Mtb benefits the TB progression to active TB disease.

Fig. 4

Clearance and development of infection with different speeds

Fig. 5

Simulated Mtb concentration as the result of the interplay between the host defense mechanism (macrophage apoptosis) and virulent Mtb strains counteract strategies (inhibition of macrophage apoptosis and trigger macrophage necrosis)

TB outcomes with a small bacterial load

Pathogen clearance is possible if , due to the local stability of the disease free equilibrium . To represent the low bacterial level due to low inhaled bacterial load or antibiotic drug treatments, we choose a low initial bacterial load cells/ml for the simulations in Fig. 4. The simulated Mtb infection clearance are shown in Fig. 4 panel (a) for six different macrophage apoptosis rates satisfying day. Simulations in Fig. 4a also demonstrate a positive relation between infection clearance time (and the Mtb peak magnitude) and infected macrophage apoptosis rate b. It is probably because the release of intracellular Mtb pathogens during the Mtb-containing macrophage apoptosis process causes the increase of the total bacterial load, which takes the immune system longer time to eliminate.

Increasing the Mtb-containing macrophage apoptosis rate passing day indicates the occurrence of heavily infected macrophage necrosis. In this scenario, the host immune system loses control of the growth of the invading pathogen. Either a small amount of inhaled Mtb pathogens or a small load of latent Mtb bacteria will eventually progress to active TB disease (see Fig. 4b). A positive relation between the speed (and magnitude) and infected macrophages apoptosis rate b is also observed in Fig. 4b. This result supports the experimental finding that the release of Mtb pathogens through necrosis induces the high Mtb bacterial burden (Lee et al. 2011).

Overall, simulation results confirm the experimental findings that the infected macrophage undergoing apoptosis promotes the clearance of Mtb infection (Keane et al. 1997, 2000; Lee et al. 2009) even though the released Mtb pathogens increase the total bacterial burden and result in longer clearance time. Moreover, the apoptosis for heavily infected macrophages culminates in necrosis, which releases more intracellular Mtb pathogens (Lee et al. 2009). The uncontrolled infection prolongs clearance time and facilitates infection progression to active TB disease.

The tug of war between macrophage apoptosis and TB infection in LTBI state

Latent tuberculosis infection (LTBI) describes an Mtb-infected individual who does not develop active disease and is not contagious. In this stage, the Mtb-infected macrophages provide a suitable environment for the latency and persistence of intracellular bacteria. However, diseases impairing the immune system (such as AIDS), treatments causing immune system suppression (such as chemotherapy for cancer), or advanced age can trigger LTBI reactivation. The LTBI status is demonstrated as low levels of Mtb pathogen load in Fig. 3 when the macrophage apoptosis rate is between LP and LP. Low levels of Mtb load can be both persistent (b between LP and or between and LP) and oscillatory (generated from two Hopf bifurcations and ). Interestingly, the oscillations won’t persist in the whole parameter interval of (, ), but disappear when oscillations starting from and starting from merge with the lower branch of the infected equilibrium (see Fig. 3b). It follows blowing-up oscillation periods (see Fig. 3c), where the coexisting unstable branches of infected equilibrium repel nearby solutions to the locally asymptotically stable disease free equilibrium . A macrophage apoptosis rate window (approximately day) between two Hopf bifurcations and are generated, where oscillations disappear. Moreover, the risk to progress to active disease is possible if the macrophage apoptosis rate b is between LP and LP, and if large load of Mtb pathogens is released due to the weakened immune system.

The role of macrophage apoptosis on the Mtb pathogen dynamics can be demonstrated by time-history simulations (see Fig. 5) based on the analytical predictions in bifurcation diagrams in Fig. 3. The clearance of internalized Mtb pathogens after phagocytosis by the acidification is demonstrated (see Fig. 5a) by a low level of macrophage apoptosis rate b near LP ( and day) and a low initial Mtb concentration in LTBI state. To confront the parasite pathogen after the infection establishment, Mtb-containing macrophages adopt a defense strategy to activate the host cell programmed death (macrophage apoptosis). This host defense mechanism is demonstrated by increasing the macrophage apoptosis rate b to the window (approximately day). Disease clearance is shown in both LTBI case (with low initial Mtb load) for day and day) (see Fig. 5b) and high initial Mtb load case for day (see Fig. 5c). Infected macrophage apoptosis effectively inhibits intracellular Mtb reproduction. To counteract this survival disadvantage, Mtb pathogens evolve virulent strains, which are capable of suppressing Mtb-containing macrophage apoptosis and continuing intracellular replication. The reduced macrophage apoptosis rate by virulent Mtb strains is between day (LP) and 0.163 day. The low level persistent and oscillatory Mtb levels are shown in Fig. 5a for day and 0.043 day). Virulent Mtb strains manage to suppress macrophage apoptosis for the intracellular replication, but also provoke host macrophage necrosis when the intracellular pathogen load passes a threshold or the host cell resources are depleted. The rapid macrophage apoptosis progresses to macrophage necrosis, which induces a large amount of intracellular pathogen release and benefits the infection. Figure 5b and c show that the low-level oscillatory Mtb level for day and day ( day) and an active disease case for day ( day) if the patient is exposed with high level Mtb pathogen. Even undergoing host cell necrosis, the immune system can still manage to control the growth of Mtb pathogens (shown in Fig. 5 for day and day), but fails if the infected macrophage loss rate passes through the threshold (BP) (shown in Fig. 5 for day).

Our analytical and numerical results confirm the mechanism that the host immune system employs programmed cell death (apoptosis) as a strategy to defend against intracellular Mtb parasitism (Lee et al. 2009). However, Mtb pathogens counteract macrophage apoptosis by evolving virulent Mtb strains. Virulent Mtb strains can both inhibit macrophage apoptosis to facilitate intracellular Mtb replication and trigger heavily-infected macrophage necrosis to release intracellular pathogens (Zychlinsky 1993; McCormick 2008; Lee et al. 2009). This complex parasitic relationship determines the various TB disease outcomes. Moreover, a comprehensive understanding of the interplay between macrophage apoptosis and Mtb pathogens is fundamental for the TB treatment design.

New treatments suggested by the new insights of macrophages programmed cell death

The current standard treatment of TB involves four first-line antibiotics (WHO 2019). Effective treatment for curing TB, however, can be hampered by the development of drug resistance, which can occur when treatment is not properly completed or adhered to. Host directed therapy (HDT), which includes small molecule therapies that can modulate the host immune response, can be used to augment treatment regimens to better control TB infection of both antibiotic sensitive and resistant TB strains (Kolloli and Subbian 2017). Antibiotic and HDT treatment can increase bacterial elimination (), infected macrophage loss (b, ), T-cell generation (, , ), and decrease macrophage infection (), and intracellular and extracellular bacterial proliferation (, , ). Sensitivity analysis has shown that , , , , and do not significantly affect the reproduction number (Fig. 1) or LTBI reactivation (Fig. 2). These parameters are therefore ignored, and we proceed with an analysis of b, , , and and , which significantly affect the reproduction number and LTBI reactivation.

Figure 6 plots two parameter bifurcation diagrams that consider parameter b with , , , and . For convenience in the bifurcation analyses, we denote the TB model (1) as , where , , and . Here, denotes the two chosen bifurcation parameters, which are one of the following pairs: , , , and . represents the remaining 16 parameters, which take the corresponding fixed values from Table 1. We represent equilibrium solutions as . Single-zero eigenvalue bifurcation happens if with positive Hurwitz conditions , for . Hopf bifurcation occurs when (1) , (2) , (3) (positive determinant of the corresponding Jacobian matrix), and (4) . Single-zero (saddle-node) bifurcation curves are denoted as LP in blue. Transcritical bifurcation satisfying (4) are denoted as BP in green. Hopf bifurcation curves are plotted in red. Here, LP, LP, and LP correspond to the three limit point (saddle-node) bifurcations in Fig. 3a. BP is consistent with the branching point (transcritical) bifurcation and Hopf matches with the two Hopf bifurcations. Note that LP and LP denote two different branches of equilibrium solutions in Fig. 3a. Therefore, even though the parameter curves for LP and LP intersect (see Fig. 6) in the corresponding 2-dimensional parameter space, the corresponding equilibrium solutions are different. Therefore, there is no intersection in the parameter-state space. It follows that no higher-codimension bifurcation occurs. We note that Bogdanov-Takens, generalized Hopf, and zero-Hopf bifurcations are detected for the bifurcation analyses b vs and b vs . However, the bifurcation critical points lie outside of biologically feasible regions from Table 1, and therefore, are not considered.

Fig. 6

2-dimensional bifurcation diagrams of b vs , , and

Figure 6 first shows that with a fixed value for , , , and , the order of bifurcation occurrence matches the order shown in Fig. 3a – as b increases, we have LP, Hopf, LP, Hopf, LP, and BP. Panels (a) and (b) also show that as b changes, the infection outcome does not vary widely, given change in and . Therefore, if drug therapy moderates the macrophage apoptosis rate b, very large changes in Mtb infection rate or Mtb proliferation rate are needed to see a change from a current stable equilibrium or periodic cycle. Additionally, to achieve bacteria elimination and prevent reactivation of LTBI, for a given or , b must be decreased to a value smaller than LP. As the macrophage apoptosis rate b changes, changes in cell-mediated immunity rate and average number of released bacteria can also have an effect. Here, we identify the disease clearance region (in Fig. 6 above LP in (c) or below LP in (d)) and LTBI controllable region (in Fig. 6 above LP in (c) or below LP in (d)). To enlarge the disease clearance and LTBI controllable regions, we need to increase b and increase or decrease . It means that a drug therapy that moderates the macrophage apoptosis works better with the combination of an increased killing rate of infected cells by the immune system or a decreased Mtb released number during programmed cell death.

In drug development, research on apoptosis modulation indicates therapeutic value in cancer (Reed 2004), cardiovascular disease, and other diseases (Cotter et al. 2003). Our findings enhance the idea that a drug therapy promoting macrophage apoptosis, but inhibiting the progression to necrosis induced by heavily-infected cell death, could benefit host defense against Mtb. In the mean time, we should consider the pathogen anti-apoptosis strategies, which could increase bacillary viability, then increase the number of released intracellular Mtb . For example, Mtb virulent strains are capable of inhibiting host macrophage apoptosis for intracellular reproduction by hindering the tumor necrosis factor (TNF). A negative therapeutic effect is associated with treating TB with infliximab, a TNF-neutralizing agent (Keane et al. 2001). Moreover, adaptive immunity is promoted by the effect of the efferocytosis process, in which Mtb-infected apoptotic corpses are removed by phagocytic cells. These phagocytic cells then present Mtb antigen and active CD8+ T cells. This phenomenon not only explains the positive effect of macrophage apoptosis on host, but also inspires vaccine development to promote adaptive immunity (Lee et al. 2009).

In conclusion, to prevent LTBI reactivation and to eliminate TB disease, the development of apoptosis modulation drugs must take the Mtb survival strategies into account. Furthermore, an understanding of the enhanced cell-mediated immunity gives new insight into vaccine development.

Conclusion

In this paper, we investigate the macrophage and Mycobacterium tuberculosis (Mtb) dynamics through an established 4-dimensional host-pathogen model. We pay particular attention to identify the determining mechanisms for Mtb elimination, LTBI reactivation prevention, and active TB treatment.

The first goal is to identify the determining factors for various disease outcomes that are unique features for TB infection. These include clearance, LTBI, active disease, and LTBI reactivation. Various disease outcomes in the mathematical model indicate the co-existing Mtb steady states levels (equilibriums). The basic mathematical mechanism for the creation and emergence of steady states is the occurrence of saddle-node bifurcation. Additionally, the fundamental mathematical mechanism for the disease clearance and persistence is the occurrence of transcritical bifurcation. Both saddle-node and transcritical bifurcations are single zero-eigenvalue bifurcation. The problem of identifying parameter sensitivity for TB elimination and LTBI reactivation is replaced by the problem of identifying parameter sensitivity for the occurrence of single zero-eigenvalue bifurcation. Note that pitchfork bifurcation is the other single zero-eigenvalue bifurcation, but is not common in infectious disease models due to its symmetric nature. It is therefore left out of our consideration.

We then perform a sample-based sensitivity analysis for the condition of transcritical bifurcation (the basic reproduction number) on all eleven model parameters. Our results demonstrate that eight of the eleven parameters significantly influence the basic reproduction number (). The eight parameters are the infection rate , T cell-mediated immunity rate , extracellular bacterium proliferation rate , infected macrophages apoptotic rate b, uninfected macrophages bacterial killing rate , and average number of bacterium engulfed by uninfected macrophages or released through the death of infected macrophage induced by apoptosis and cell-mediated immunity . Fixing the parameters that do not statistically significantly affect the basic reproduction number, another sensitivity analysis is carried out for the conditions of saddle node bifurcation on the eight parameters that significantly influence . These two steps of sensitivity analyses find that the infected macrophage cell death rate b is the most influential factor on both and the occurrence of saddle-node bifurcation, which determines the occurrence of various disease outcomes and, in turn, rules the fate of the disease. Biologically, the host cell death in the form of apoptosis is considered to inhibits Mtb replication and to produce bacterial antigens to activate adaptive immune responses (Behar et al. 2011; Queval et al. 2017). On the other hand, cell death in the form of necrosis is associated with spreading infection (Behar et al. 2011; Queval et al. 2017). Other proposed host factors and Mtb invading strategies, which contribute to the disease progression, include genetic diversity (Gagneux 2013), T-cell mediated immune responses (Verrall et al. 2014; Flynn and Chan 2001), the responses of two pivotal cytokines gamma interferon (IFN-) and tumor necrosis factor alpha (TNF-) (Flynn and Chan 2001), Mtb inhibiting the elimination of phagocytized bacteria Queval et al. (2017), and Mtb controlling host cell death to escape from kaltilling Queval et al. 2017. One limitation of our mathematical results is the choice of parameter values. The parameter ranges in Table 1 are based on current experimental and clinical literature, but still inherit their restrictions. Nevertheless, our mathematical results confirm the crucial and controversial role in the host immune defense against invading Mtb pathogens (Queval et al. 2017).

Sensitivity analysis dramatically narrows down parameters of interest for understanding model behaviors (including infection elimination and LTBI progression), but it can not tell how the identified parameters affect the model behaviors, or what model dynamics can be unfolded from the change of the identified parameters. Answering these two questions serves as our second goal. We employ bifurcation analysis to tackle these questions. Considering the macrophage apoptosis rate b as a bifurcation parameter, we explore the role of macrophage apoptosis in TB infection through bifurcation analyses. A spectrum of disease outcomes are verified by numerical simulations. We carry out bifurcation analyses on single-zero eigenvalues bifurcations and derive Hopf bifurcation normal forms. Our analytical results are consistent with the numerical bifurcation output from MatCont (Dhooge et al. 2003).

Regarding TB pathobiology, our findings of the mathematical model dynamics agree with the clinical and experimental macrophage-Mtb interaction outcomes (Lee et al. 2009), which consider the complex roles of macrophage apoptosis in the host defense against TB infection. Our mathematical model (1) demonstrates the following: (1) the clearance of low level Mtb load during phagocytosis is associated with low macrophage apoptosis rate (Verrall et al. 2014); (2) the TNF-induced macrophage apoptosis facilitates the host defense against intracellular Mtb proliferation (Lee et al. 2009); (3) macrophage apoptosis is suppressed by virulent Mtb strains interfering with TNF signaling pathways. Virulent Mtb then gain an intracellular reproduction environment (Kornfeld et al. 1999); and (4) heavily infected macrophage apoptosis induces macrophage necrosis, which releases a large amount of intracellular bacteria and benefits the spread of the pathogen (Lee et al. 2009; Behar et al. 2011).

Elucidating the complex roles of macrophage apoptosis played in Mtb in-host dynamics could give new insights into the treatments aiming at TB elimination and LTBI reactivation prevention. Treatment strategies regulate macrophage apoptosis in combination with the control of bacillary viability is recommended. Mtb antigen presentation after the efferocytosis of Mtb-containing apoptotic remaining enhance adaptive immunity by activating CD8+ T cells. The adaptive immunity promoting mechanism can help the development of the new vaccine and new treatments.

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 1849 KB)