Mathematical Model of the Immunopathological Progression of Tuberculosis

1 Introduction

Tuberculosis (TB) is a worldwide persistent infectious disease caused by bacteria from the Mycobacterium tuberculosis complex that mainly affects the lungs and can be fatal (Pai et al., 2016). It is one of the top 10 causes of death worldwide and approximately a quarter of the world’s population is latently infected (WHO, 2019). Its high prevalence is, among other factors, a result of the ineffectiveness of treatments for the prevention and reversal of the infectious process. Such ineffectiveness is attributed to the heterogeneous constitution of the bacterial genotypes, antibiotic-resistant strains, and expensive treatments (Smith, 2003). Efficient treatments are difficult to establish as some stages of TB progression are clinically subtle and there is a lack of understanding of the molecular mechanisms behind them (Salvatore and Zhang, 2017). Using an integrative systems biology approach, we can study the immunopathological progression of this disease, describing the dynamics of the host–pathogen interactions in a qualitative and quantitative way that goes beyond the current capacity of experimental studies (Reynolds, 2014).

The holistic approach of systems biology in medicine offers a unique way to develop personalized treatment strategies, taking full advantage of clinical and experimental data sets (Loscalzo and Barabasi, 2011). From mathematical models that predict and explain the molecular mechanisms behind skin diseases such as atopic dermatitis (Dominguez-Hüttinger et al., 2017) and upper airway infection by Streptococcus pneumoniae (Domínguez-Hüttinger et al., 2017a) to modeling the dynamics of signaling networks behind the development and progression of cancer (Du and Elemento, 2015), there is a wide range of pathologies that have already been addressed by the discipline (Barber et al., 2021; Reynolds et al., 2006). Mathematical TB models have been proposed to represent the transmission of the disease at the epidemiological level (Trauer et al., 2014) and to assess the impact of risk factors and treatments on the population (Dowdy et al., 2008). Agent-based models that describe the progression of TB infection at a tissue level have been formulated (Wigginton and Kirschner, 2001; Segovia-Juarez et al., 2004). At the cellular level, ODE models that describe the host–pathogen interactions have been formulated (Wigginton & Kirschner, 2001; Magombedze et al., 2006). However, a minimal ODE model that describes how the host–pathogen interactions are dynamically regulated through the different phases of TB progression has not been formulated until now.

In the present work, we propose an ODE model that describes the interactions that occur at cellular level throughout the immunopathological progression of active pulmonary TB in an in vivo mouse model. Specifically, we model the interactions between bacteria and macrophages, which are the main effectors of the TB progression, and how these interactions are shaped by adaptive immune responses that dynamically change as TB progresses. The model is built, calibrated, and validated by integrating a vast number of in vivo and in vitro published experimental assays (Supplementary Table S1). The model mechanisms replicate the qualitative behaviors of the macrophage and bacterial population dynamics observed in in vivo studies, with an initial slow bacterial growth accompanied by gradual macrophage recruitment (phase 1), followed by a sharp rise in both macrophages and bacterial population (phase 2) and finally macrophage population decreasing their protection efficiency and bacteria reaching its maximum values (phase 3). Furthermore, our model captures how these phases are shaped by the key adaptive immune response players, namely Th1, Th2, and dendritic cells and protective cytokines such as TNFα.

The model can predict possible disease outcomes for different genotypes of mice by simulating the interaction between a range of host genotypes with greater or lesser immune response capacity and a range of pathogen genotypes with different degrees of virulence through parameter modulation. Two disease outcomes emerged, bacterial clearance and uncontrolled bacterial growth, and this all-or-nothing behavior persisted through all the phenotype space. The probability of ending up with an uncontrolled TB infection showed up to be very high (84%). This is in accordance with the progressive pulmonary TB experimental model that our model aims to replicate, which is designed so that a successful progressive TB infection is established on the mouse, resembling the active pulmonary TB of humans. The transitions across phases of disease progression show that in all cases, adaptive immune response (phase 2) needs to be recruited and that an intense anti-inflammatory response (phase 3) might be detrimental to the organisms in terms of bacterial clearance. Additionally, we sampled host and pathogen genotypes from the three different phases and studied their long-term behavior in absence of disease progression. Almost half of these genotypes showed a bistable behavior, one being the stable behavior, bacterial clearance, and the other one is the uncontrolled bacterial growth. Genotypes sampled from phase 3 (anti-inflammatory) showed only a 7% chance of bacterial clearance; meanwhile, phase 1 (preparation) and phase 2 (proinflammatory) showed a 32% and 53% chance, respectively. Thus, genotypes sampled from phase 2 have the highest probability of mapping to an infection clearance phenotype across the three phases.

This model provides a powerful tool to test the effect of specific alterations of the host–pathogen interactions over the outcome of the infection. These in silico experiments can lead to future experimentation and help reduce the number of in vivo experiments.

2 Results

2.1 Mathematical Model of the Immunopathological Progression of Tuberculosis

Our mathematical model of the progression of a TB infection (Figures 1A, B) is a systems-level representation of the key interactions between the effector cells of the immune response—namely, macrophages—and the pathogen observed in vivo on mouse models.

FIGURE 1

FIGURE 1. Mathematical model of the immunopathological progression of tuberculosis. (A) Stages of the immunopathological progression of tuberculosis infection in mice. (B) Graphical representation of the network of interactions of the population dynamics in granulomas. (C) Three phases of disease progression. We integrated counts of TH1 (orange) and TH2 (red) cells in granuloma (Hernández-Pando et al., 1996), dendritic cells in lung (yellow) (Pedroza-González et al., 2004), and concentration of iNOS mRNA expression that represents M1 macrophages (orange) (Hernández-Pando et al., 2001). In all cases, the BALB/c mouse was inoculated with 10⁶ M. tuberculosis H37Rv bacteria. The data were normalized to their maximum percentage to allow comparison. (D) Example of how an adaptable parameter can change through the different phases of the infectious process.

The first stage of the infectious process is the primary infection (Figure 1Aa). After natural bacterial infection in humans or inoculation of the mouse, generally through the aerial route (intratracheal instillation or aerosol infection) (Hernandez-Pando et al., 1995; Hernández-Pando et al., 1996), the pathogen reaches the pulmonary alveoli where it is recognized by non-activated macrophages (Wolf et al., 2007). Macrophages then phagocyte bacteria, and because of their non-activated state, bacilli are not efficiently eliminated. Thus, the macrophages generally serve as a reservoir for the bacteria, which proliferate there, eventually causing apoptosis of these host cells (Cohen et al., 2018). During this early stage, bacteria stimulate the phagocytic cells, triggering secretion of cytokines such as tumor necrosis factor alpha (TNFα) which permeabilizes the vascular endothelium and attracts circulating monocytes to the site of infection, which then differentiate into non-activated macrophages upon entering the infected tissue (Figure 1Ab) (Shi & Pamer, 2011). Together with small populations of dendritic cells and polymorphonuclear leukocytes, this regulatory network formed between Tb and non-activated macrophages forms the early granuloma that characterizes phase 1 of TB (Figure 1Ac) (García-Romo et al., 2004). In cases of very severe and virulent infections, this can lead to uncontrolled bacterial growth, macrophage necrosis, and, in humans, signs and symptoms due to tissue damage such as coughing up blood, chest pain, difficulty in breathing, fever, fatigue, and cachexia, which is called primary active TB (Figure 1Ad) (Mayo Clinic, no date). On the other hand, infection can be mild and/or the immune system can be fit, ending up in the clearance of the bacterial population (Figure 1Ae). These two outcomes are possible through all the disease progression.

After there is a sufficient quantity of bacteria available in the area of infection, some will be phagocytosed by dendritic cells, which will migrate to the closest lymph nodes, where they present M. tuberculosis antigens through major histocompatibility complex type II (MHC-II) molecules to immature CD4⁺ T lymphocytes; this interaction is co-stimulated by the link between the dendritic cell membrane receptor B7 and the T cell-stimulating molecule CD28 that in the presence of interleukin (IL)-12 in the microenvironment, triggers the maturation of CD4⁺ T cells into “proinflammatory” TH1 cells (Figure 1Af) (Alberts et al., 2002). The TH1 cells then migrate from the lymph nodes to the lung infection sites and secrete IFNγ, which, together with co-stimulation by interactions between CD40 molecules on the surface of macrophages and CD40L on the surface of TH1 cells, activates the macrophages, differentiating them into M1 macrophages (Hernández-Pando et al., 1996), which are characterized by high production of reactive oxygen species (ROS, oxidative burst), maturation of their phagolysosomes, and production of antimicrobial molecules (RNI, NO, NOS2, defensins, etc.) (Ley, 2017). As a result of the effects of the recruited adaptive immune response, the early granuloma will become a mature granuloma, which is more efficient in eliminating bacteria thanks to its activated M1 macrophages (Figure 1Ag) (Hernández-Pando et al., 2001). Due to the high burden of bacteria that is necessary to induce experimental models of progressive pulmonary TB ($>1\times {10}^5$) as mice are naturally resistant to Mycobacterium tuberculosis infection (Rook et al., 2009), the early and mature granuloma are typically unable to contain the bacterial growth and the infection, which therefore progresses to a phase called primary progressive TB infection in which an uncontrolled infection and tissue damage are developed, resembling human active pulmonary TB. As a result, a late granuloma is formed (Figure 1Ah), in which the phenotype of lymphocytes changes, with a predominant number of TH2 cells over TH1 cells, as a response to the tissue damage. Through the production of TH2 cytokines by T cells, macrophages differentiate into the M2 phenotype that performs anti-inflammatory, immunoregulatory, and tissue preservation activities (Ley, 2017). The late granuloma turns out again to be less efficient than the mature granuloma in eliminating bacteria as a result of the phenotype change in macrophages from M1 to M2.

Although there are multiple host cells that participate in the onset and progression of TB (Davis & Ramakrishnan, 2009), macrophages are the direct effectors of the immune response (Hernández-Pando et al., 2001). Therefore, the interplay between these cells and M. tuberculosis is the key host–pathogen interactions that occur in TB at the cellular level.

Our mathematical model is a formal representation of the dynamical interplay between alveolar macrophages and M. tuberculosis that occurs in granulomas (Figure 1B). The ODE model consists of four state variables: free macrophages ($M$), free bacteria ($T$), macrophages that have phagocytosed ($Mf$), and phagocytosed bacteria ($Tf$). While we aimed to keep the model as simple as possible, we carefully represented the mechanisms of macrophage deaths and phagocytosis. The model distinguishes between many types of macrophage deaths, the natural death of $M$ and that of $Mf$, both of which occur independently of the presence of bacteria but can also be enhanced by them. The deaths are caused by free or phagocytosed bacteria. Death by apoptosis, which is an “ordered” form of death, also contributes to the elimination of $Tf$. Finally, during death by necrosis, $Mf$ discharges its cellular content to the extracellular medium, releasing $Tf$ that becomes $T$ and also stimulating phagocytosis due to damage-associated molecular patterns (DAMPs) released in the process (Roh & Sohn, 2018). Phagocytosis of $T$ by $M$ for the creation of $Tf$ and $Mf$ occurs irreversibly although necrosis allows the release of $Tf$ back to the extracellular space. Both free and phagocytized bacteria can die naturally, and $Tf$ can also die by the antimicrobial activity of $Mf$. Bacteria can proliferate up to a carrying capacity at the site of infection; since $Tf$ needs a $Mf$ to exist, its carrying capacity and proliferation is $Mf$ dependent. Macrophages are recruited from the bloodstream and do not proliferate de novo; this recruitment depends on the presence of $T$ since only $Mf$ produces the necessary cytokines for it. Further details of the reactions considered in our model given in Figure 1B can be found in the Methods section. The resulting system of equations is given in Eq. 1; the model parameters are in Tables 1 and 2.

TABLE 1

TABLE 1. Fixed nominal model parameters.

TABLE 2

TABLE 2. Modulable model parameters.

Tuberculosis is a progressive disease, with immune responses changing their phenotype over time. Our in vivo data set reflects these changes, with the three clearly distinguishable phases shown in Figure 1C. As the main role of the dendritic cells is to sense bacteria, the pathogen levels are represented through these cells. Initially, in phase 1, no adaptive immune cells are present in the system. Phase 2 shows a clear peak in TH1 cells, which play a proinflammatory and antimicrobial role. TH1 cells then decrease while the anti-inflammatory Th2 cells increase in phase 3. These three phases correspond, thus, to a preparation phase on days 0–6, an inflammatory phase on days 7–27, and an anti-inflammatory phase on days 28–120. Phase 1 corresponds to the early granuloma, where the non-active macrophages interact with the pathogen and prepare the recruitment of the adaptive immune response. Phase 2 corresponds to the mature granuloma and a large number of M1 macrophages that have been recruited and activated by the TH1 cells. Phase 3 corresponds to a late granuloma, with a higher number of TH2 cells at the site of infection, leading to a sharp decrease in the number of M1 macrophages accompanied by a dramatic increase in M2 types of macrophages.

The model represents these three phases by dynamically changing kinetic rates that are affected by macrophage differentiation in response to the accumulating bacterial population. These rates, denoted by dotted arrows in Figure 1B, are the levels of recruitment, phagocytosis, and antimicrobial capacity which change dramatically between M1, M2, and not activated macrophages. Additionally, their apoptosis and necrosis rates are modulated since both cell deaths appear in different proportions at different stages of infection. Since bacteria adapt to the challenge imposed by the immune system (Montoya-Rosales et al., 2017), their killing capacity is modulated throughout the phases. Parameters change dramatically between each phase of disease progression (Figure 1D) to emulate the abrupt changes in phenotypes that occur in in vivo studies (Figure 1C).

The optimal values of the 16 parameters of the model were obtained using experimental data, either from direct measurements, estimated from quantitative in vitro studies by implementing simplified versions of our in vivo model or by parameter optimization of the full model from in vivo studies, as detailed in Methods section. Most of the data come from our own wet lab, minimizing the batch effects and therefore the margin of error in the estimation since, despite being different studies, the experimental conditions such as mouse strain, bacterial strain, number of bacteria used for inoculation, and frequency of sampling, among others, remain constant. Tables 1 and 2 describe the parameters for each reaction.

In Figure 2A we see the result of simulating our model with optimal parameters. We plot the trajectories of the compound variables of the model $Mtot$ ($M$ + $Mf$) and $Ttot$ ($T$ + $Tf$) against the data set used to perform the optimization. The model recovers the qualitative behavior of the cells throughout the progression of the disease. We were unable to replicate these behaviors without modulating the parameters between phases (Supplementary Figure S1).

FIGURE 2

FIGURE 2. Parametric optimization and validation of the model. (A) Model replicates in vivo behavior. The model (solid line) was plotted together with the in vivo experimental data set (points) used for optimization. On the vertical axis, we see macrophages and total bacteria, and on the horizontal axis, we have time. (B) Behavior of the modulable parameters throughout the phases is biologically accurate. Normalized parameter values are on the vertical axis and time on the horizontal. The dotted lines indicate the phase change threshold, depending on the case, the parameter will increase or decrease its value when passing this threshold.

We plot the changes in the modulable parameters’ values over time (Figure 2B). Depending on their function, the parameters will decrease or increase between each phase. In phase 1, macrophage recruitment (${\alpha }_1$) and macrophage deaths caused by bacteria (${\beta }_4$, ${\beta }_5$, and ${\beta }_7$) are low while the antimicrobial activity of macrophages ($\delta$ and $\gamma$) is null. Apoptosis of $Mf$ (${\beta }_3$) is the only parameter that reaches its peak in this phase since macrophages are not activated and opt for apoptosis to deal with bacteria (Bocchino et al., 2005). Necrosis (${\beta }_6$) is a mechanism of uncontrolled death (Golstein & Kroemer, 2007), so when macrophages go for apoptosis, it will be null. In phase 2, we see peaks in antimicrobial activity ($\delta$ and $\gamma$) and recruitment of macrophages (${\alpha }_1$). Apoptosis decreases as M1 macrophages are more efficient at killing bacteria. We also see a peak in the activity of bacteria (${\beta }_4$, ${\beta }_5$, and ${\beta }_7$) that corresponds to the ability of the pathogen to respond to a more aggressive immune system (Montoya-Rosales et al., 2017). In phase 3, there is a massive drop in macrophage parameters associated with their antimicrobial ($\delta$ and $\gamma$) and inflammatory (${\alpha }_1$) activity, distinctive characteristics of M2 macrophages. It is worth noting that other macrophage populations different from the alveolar macrophages that we are modeling, like vacuolated macrophages, may be recruited to the infection site at late tuberculosis stages (Ríos-Barrera et al., 2006; Garcia-Romo et al., 2013). The only macrophage parameter that peaks in this phase is necrosis. This is probably due to extreme stress conditions caused by the uncontrolled growth of bacteria (Poon et al., 2010). Bacteria reduce all their reaction strengths in this final phase, and this decrease responds to a weaker immune system that presents a minor challenge. The biological accuracy observed through all parameter value changes supports the validity of our model.

2.2 Construction of a Disease-Transition Map for Predicting the Most Likely Outcomes of the Immunopathological Progression of Tuberculosis

To explore possible disease outcomes resulting from variations in the host and pathogen genotypes, we performed simulations of the in vivo model under parameter variations. We produced a cohort of 500 virtual mice by varying the 16 parameters around the nominal values corresponding to the optimal set of phase 1 and performed numerical integrations for each randomly generated set. As before, we allowed the modulable parameters to adapt in response to the accumulating bacterial population, which resulted in virtual disease trajectories going through a variable number of phase transitions. Figure 3A shows simulations whose last phase was phase 2 or 3. All sampled genotypes progressed from phase 1 without clearing the infection, which indicates that there is no sampled genotype capable of eliminating the bacterial population without recruiting the adaptive immune response. Since phase progression dynamically depends on bacterial numbers, simulations stationed in phase 2 represent those genotypes with an adaptive immune system capable of eradicating the infection. Additionally, we see that $Mtot$ slowly decreases upon having eradicated the bacteria, implying that the model is capable to return to macrophage population homeostasis numbers after eradicating the infection. Simulations stationed at phase 3 represent a mature granuloma that recruited an anti-inflammatory response and became a late granuloma. The model with the optimal parameter set (Figure 2A) is grouped here, however, not all the simulations in this group show the uncontrolled bacterial growth observed for the nominal parameter set. We can see this in the levels of $Mtot$, where some simulations grow above their initial condition and then begin to decrease (like our optimal set) while others decrease from the beginning and never exceed their initial condition. The trajectories of $Mtot$ that more rapidly decrease to zero correspond to simulations in which the bacteria established a hyper-aggressive infection (Supplementary Figure S2). The fact that some $Ttot$ trajectories reach $K$ before others correspond to variations in the virulence of the bacteria in combination with an immune system with varying degrees of efficacy.

FIGURE 3

FIGURE 3. Disease-transition of the immunopathological progression of tuberculosis. (A) Compound variables $Mtot$ and $Ttot$ plotted against time, the continuous line represents the mean of the total simulations performed, and the light grey area represents the standard deviation. The simulations are grouped by the final phase they reached. $Ttot$ was normalized by its carrying capacity $K$. (B) Transition times from phase 1 to 2 (critical time 1) plotted against the transition times from phase 2 to 3 (critical time 2), the red point represents the values obtained with the nominal parameters, and the blue points represent the random genotype samples. (C) Disease progression map with outcome probabilities and estimated transition times.

In Figure 3B, we analyzed how the critical times in which a phase transition occurs are distributed for the sampled genotypes. Recall that these critical times correspond to the times at which bacteria reach a population threshold and thereby activate a response from the host. While in the original dataset, the model enters phase 2 on day 7, in general, this critical time 1 is more densely distributed between days 1 and 2. As for the critical time 2, our simulations show that late granulomas can already form within a week post-bacterial challenge. These accelerated immunopathological transitions are consistent with the rapidly progressing disease dynamics observed for more aggressive bacterial strains, such as the Beijing strain (Aguilar L et al., 2010). Critical times that are higher than the nominal ones correspond to organisms that opposed greater resistance to the progression of the disease.

When analyzing the disease outcomes on all the genotype network (Figure 3A), only two phenotypes persisted: bacterial clearance and uncontrolled bacterial growth. This all-or-nothing behavior is consistent with the experimental model of progressive pulmonary TB in mice, where there is no coexistence between mice and bacteria (Garcia-Romo et al., 2013; Bini et al., 2014). Information about disease outcomes across all the genotype networks grouped with the critical times allowed us to build a predictive disease progression map (Figure 3C). We see that all the simulations transitioned from phase 1 to phase 2. All the simulations that avoided phase 3 cleared the infection, and those that transitioned to the final phase developed an uncontrolled bacterial growth. These results show that when starting from the initial conditions used in an experimental model of progressive pulmonary TB in mouse, it is 83% likely that the mice will suffer an aggressive symptomatic disease under variations of the mice and the bacteria genotypes (strains). Also, the fact that all the simulations that transitioned to phase 3 have an uncontrolled bacterial growth suggests that the anti-inflammatory response is a key determinant in the host–pathogen interaction. Only those simulated disease trajectories that effectively cleared the infection could avoid transitioning to phase 3.

2.3 Long-Term Behavior of the Different Phases Reveals All-or-Nothing Outcomes

To explore the influence of the immune response phenotype on the disease outcome, we sampled genotypes from the three different phases and studied their long-term behavior in a hypothetical scenario where the organism can start from any phase and there’s no change of phases. We varied all parameters from each optimal set and produced 500 random sets from each (1,500 total), and then equilibrium points for each set were obtained exactly (for details see Methods section).

In accordance with the dynamical simulations (Figure 3), the only long-term disease outcomes found was the all-or-nothing behavior where the system either clears the infection or allows bacteria to grow out of control; Figures 4A, B show the frequency in which simulations sampled from each phase end up in either of these two possible outcomes. Phase 2 showed the highest probability of bacterial clearance among phases (53%), suggesting that a recently stimulated immune system that still has active adaptive immune response cells has a chance of clearing even an aggressive infection, such as the experimental model of progressive pulmonary TB. The BCG vaccine, the most common prevention method for TB infection used worldwide, is believed to generate protective immunity through this mechanism (Ritz et al., 2008). Genotypes locked in phase 1 showed a 32% chance of clearing the infection; however, these genotypes needed between 10⁸ and 10¹¹ macrophages to achieve the task while phase 2 genotypes needed 10⁵–10⁹ and performed better. This demonstrates that macrophage activation is crucial for bacterial clearance and for preventing excessive inflammation, which can also result in disease (Furman et al., 2019). Finally, phase 3 phenotypes had a 7% chance of clearing the infection, and, as expected, this anti-inflammatory phase had the lowest macrophage load needed to clear the infection (10⁵–10⁷). The low chance of infection clearance may be a price too high to pay for the tissue-preserving anti-inflammatory response recruitment.

FIGURE 4

FIGURE 4. The disease outcome shows an all-or-nothing phase-dependent behavior. (A,B) Steady-state distribution of disease outcomes shows the probability of bacterial clearance and of uncontrolled growth and its associated final macrophage concentrations in a basal (phase 1), inflammatory (phase 2), and anti-inflammatory (phase 3) environment. A set of 500 genotypes was sampled for each phase. (C) Proportion of mono and bistability found in each phase, the clear part of the pie graph indicates monostability and the filled bistability. (D,E) Distributions of the differences of bacteria and macrophages between the two stable steady states for the bistable genotypes.

We also saw the appearance of bistability in the system, which is present in all phases and is more frequent in phase 2 (Figure 4C). Also herein, we observe a dramatic phenotypic segregation, with all the bistable cases displaying an all-or-nothing behavior with one stable steady state corresponding to uncontrolled bacterial growth ($\mathrm{\Delta }Ttot=K$) and the other one to bacterial clearance (Figure 4D,E). This dramatic difference between the two stable steady states is also reflected in the total macrophage load, which spans several orders of magnitude, with the highest difference for phase 1 and the lowest for phase 3. This indicates that, for genotypes presenting bistability, it is possible to achieve complete recovery in a finite period of time through treatment and that initial conditions are key in determining the disease outcome.

3 Discussion

The mechanisms behind the immunopathological progression of TB have not been fully understood, and this makes the creation of more efficient treatments difficult. Based on experimental data from the literature, we propose a minimal mathematical model that represents the key regulatory mechanisms that can account for the pathological progression of TB through the three phases observed in in vivo mouse models. The model describes not only the interactions between the main effector of the immune response and the pathogen but also how these interactions are modulated by the inflammatory microenvironment. For its construction, we integrated data from a total of 75 experimental assays from different published manuscripts. The model coherently replicates the changes that the main molecular mechanisms undergo through the progression of the infection by the dynamic modulation of its parameters. With simulation analysis, we observed how the model is segregated depending on its initial parameter conditions by qualitatively different phenotypes. This model is a tool to study the interactions between the effector cells of the immune response and M. tuberculosis and allows the proposal of hypotheses related to current therapeutic problems, offering potential future use in clinical and experimental strategies for the design and search of new TB treatments.

Our model also reproduces the abrupt phenotypic changes shown by the immune system and the bacteria throughout the progression of an experimental model of progressive pulmonary TB in mouse (Figure 1C), specifically, from an early phase where there is relatively low activity in all cells involved, going by an intermediate phase characterized by the growth of the bacterial population and the increase in the recruitment of macrophages due to the initiation of the adaptive immune response, to an advanced phase where bacteria grow to their carrying capacity and alveolar macrophages lose their efficiency due to the effect of anti-inflammatory mechanisms of tissue preservation (i.e., TH2 cytokines). Even though our mathematical model has three different phases because of the experimental model it replicates, there are other types of mice models like the latent TB model (Arriaga et al., 2002) that use a significantly lower dose of bacteria (4,000 bact.) and a different strain of mice (C57B1). This produces a slow progression infection that results in an extended phase 2 and a lack of phase 3.

The phase changes are emulated by the model through the continuous modulation of eight of its parameters (Table 2). This type of dynamic adaptation of values of parameters that underlie the different immunopathological stages of a disease was proposed for the first time by Tiemann et al. (2011) and van Riel et al. (2013). Therefore, the ADAPT algorithm is proposed for adjusting a model to different data sets corresponding to different phases of a progressive phenomenon and identifying those parameters that change during the progression. This is, therefore, a strategy similar to ours. However, by dynamically modeling the changes with Hill equations that depend on the variables of our model, we are going one step further by mathematically representing these mechanisms that underlie parametric adaptation. The model was validated by analyzing the qualitative changes of the modulated parameters through the phases; indeed, we found that all changes were biologically accurate. Although we modeled the macrophage and bacteria phenotype transition in a phenomenological way, in the future, it could be carried out by coupling a Boolean network that models this phenotypical plasticity like the macrophage differentiation model proposed by Palma et al. (2018).

We explored the effects of genotypic variations on disease outcomes and progression by simulating disease trajectories from 1,500 randomly sampled parameter sets. We observed that mice are not capable of eliminating the infection without recruiting the adaptive immune response and that its recruitment does not guarantee bacteria elimination by the immune system since only for certain combinations of the enhanced immune system and not-so-virulent bacteria, the organism cleared the infection. Accordingly, experimental assays have been carried out testing a broad spectrum of M. tuberculosis genotypes that demonstrate variability in virulence and in induced immune response (Marquina-Castillo et al., 2009; Hernández-Pando et al., 2012). In particular, the Beijing genotype is of special interest because of its association with numerous outbreaks worldwide (Glynn et al., 2002) and the hyper-aggressive infection and non-conventional immunopathological events that have been shown to evoke (Aguilar L et al., 2010; López et al., 2003). Our model sets a work frame for future studies of this and other genotypes of bacteria.

Using the random parameter sets of the three phases, we obtained the equilibrium points of the system. The hypothetical long-term analysis allows us to visualize the space of possible outcomes of the model. The results indicate that the model can exist in two mutually exclusive stable states: macrophages between 10⁴–10¹¹ and bacteria at zero (clearance) and bacteria in their carrying capacity and macrophages at zero (uncontrolled growth), which clinically correspond to a severe symptomatic TB disease and to eradication of the infection, respectively. These behaviors are consistent with the experimental model of progressive pulmonary TB in mouse, where coexistence between bacteria and macrophages is not observed. Our steady-state analysis also shows that, for a significant proportion of parametric variations, the model exhibits bistability, with the stable coexistence of the two mutually exclusive outcomes of the infection. These results suggest that not only genetic variations but also the inflammatory microenvironment (phase 1, 2 vs. 3) and the initial macrophage and bacterial loads (bistable genotypes) could have a dramatic effect on the outcome of the infection. The parametric variations around the nominal values of phase 2 are the ones that have a higher probability to end up in bacterial clearance. This suggests that, from a therapeutic intervention design point of view, it may be important to promote the conditions that lead to and maintain phase 2, and this could be achieved by therapeutic interventions that strengthen macrophage recruitment and activation. On the other hand, transition to phase 3 could be blocked by preventing the activation of the anti-inflammatory response promoted by TH2 cells.

Although our model shows that, in the long term, only two disease outcomes (clearance and uncontrolled bacterial growth) are possible, it is important to acknowledge that coexistence between the host and the pathogen has been widely observed. Global clinical data shows that for humans, latent TB infection is the most common form of the disease (WHO, 2019), and this disease state is characterized by low bacteria levels residing inside fibrotic old granulomas. As our model reproduces a progressive pulmonary TB mouse infection, coexistence is never observed. However, an in silico reproduction of a TB infection could be achieved by appropriately changing the parameters and initial conditions of our model to mimic the mouse model of latent TB (Arriaga et al., 2002) that uses a significantly lower dose of bacteria (4,000 bact.) and a different strain of mice (C57B1).

Mathematical TB models that study the disease dynamics at the epidemiological level have been proposed previously, focusing on how the disease is transmitted and how treatments affect this transmission (Dowdy et al., 2008; Trauer et al., 2014), but without providing information on the underlying cellular mechanisms. At the cellular/tissue level, several agent-based models that describe granuloma formation have been proposed (Young et al., 2008; Marino & Kirschner, 2016; Pienaar et al., 2016; Evans et al., 2020). These models represent the interplay between a large number of variables compared with our four variable minimal model, making it more difficult to estimate the resulting large number of parameters from experimental data and to assess which mechanisms are the most influential on disease outcome. Additionally, mathematical analysis for in silico design and optimization of therapeutic interventions, such as bifurcation analysis as carried out in Domínguez-Hüttinger et al. (2017a), cannot easily be performed on agent-based models. At the cellular/molecular level, ODEs models that describe the host–pathogen interactions have been formulated before (Wigginton & Kirschner, 2001; Magombedze et al., 2006), but these models do not explicitly encompass how the disease progression through different phases of a TB infection emerges, how it shapes the phenotypes of both the host and the pathogen, and it dynamically affects disease outcome.

Here we propose, for the first time, a minimal, mechanistic ODE model of the immunopathological progression of TB. The model represents the dynamic modulation of the intra-host mechanisms that underly the development of TB across different disease phases. It was constructed, calibrated, and validated by integrating a wealth of in vivo and in vitro datasets. The model analysis allowed us to predict the most probable disease trajectories under host and pathogen genotypic variations and to pinpoint the key role that the initial conditions in bacterial and macrophage loads and in the inflammatory microenvironment play in determining the outcome of the disease. It is a step forward in understanding the mechanisms that shape TB progression.

4 Methods

4.1 Data Integration

Relevant experimental data to explain the immunopathological progression of TB were collected from different bibliographic sources, the biggest source being the published work from the Laboratory of Experimental Pathology at the National Institute of Medical Sciences and Nutrition Salvador Zubirán, Mexico City, Mexico. We were interested in collecting the largest amount of data from in vivo studies that used the same bacteria strain, H37Rv, and the same mouse strain, BALBc, to ensure the greatest possible coherence and minimal batch effects when integrating the data from different assays. When necessary, the raw data were extracted directly from figures using WebPlotDigitizer 4.4 (Rohatgi, 2020). The resulting database is given in Supplementary Table S1.

4.2 Model Construction

There are multiple cells of the immune response that participate in the containment of the infection; however, macrophages are the direct effectors of these responses. Therefore, we decided to put the macrophage as the centerpiece of the model. Together with bacteria, they represent the population dynamics in granulomas throughout the progression of the disease. The following mathematical model formally represents these dynamics:

$\dot{M}=\stackrel{\mathbf{R}\mathbf{e}\mathbf{c}\mathbf{r}\mathbf{u}\mathbf{i}\mathbf{t}\mathbf{m}\mathbf{e}\mathbf{n}\mathbf{t}}{\overbrace{Mf\cdot {\mathit{F}}_{\mathit{\alpha }\mathbf{1}}}}-\stackrel{\mathbf{P}\mathbf{h}\mathbf{a}\mathbf{g}\mathbf{o}\mathbf{c}\mathbf{y}\mathbf{t}\mathbf{o}\mathbf{s}\mathbf{i}\mathbf{s}}{\overbrace{M\cdot \frac{T}{PTf}\cdot {\mathit{F}}_{\mathit{\delta }}\cdot \left(1+Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{6}}\cdot \sigma \right)}}-\stackrel{\mathbf{D}\mathbf{e}\mathbf{a}\mathbf{t}\mathbf{h}}{\overbrace{M\cdot {\beta }_1}}-\stackrel{\mathbf{D}\mathbf{e}\mathbf{a}\mathbf{t}\mathbf{h}\mathbf{b}\mathbf{y}\mathit{T}}{\overbrace{M\cdot T\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{4}}}},(1)$

$\dot{Mf}=\stackrel{\mathbf{P}\mathbf{h}\mathbf{a}\mathbf{g}\mathbf{o}\mathbf{c}\mathbf{y}\mathbf{t}\mathbf{o}\mathbf{s}\mathbf{i}\mathbf{s}}{\overbrace{M\cdot \frac{\mathrm{T}}{PTf}\cdot {\mathit{F}}_{\mathit{\delta }}\cdot \left(1+Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{6}}\cdot \sigma \right)}}-\stackrel{\mathbf{D}\mathbf{e}\mathbf{a}\mathbf{t}\mathbf{h}}{\overbrace{Mf\cdot {\beta }_2}}-\stackrel{\mathbf{D}\mathbf{e}\mathbf{a}\mathbf{t}\mathbf{h}\mathbf{b}\mathbf{y}\mathrm{T}\mathbf{a}\mathbf{n}\mathbf{d}\mathbf{T}\mathbf{f}}{\overbrace{Mf\cdot \left(T\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{5}}+Tf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{7}}\right)}}-\stackrel{\mathbf{Apoptosis}}{\overbrace{Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{3}}}} -\stackrel{\mathbf{N}\mathbf{e}\mathbf{c}\mathbf{r}\mathbf{o}\mathbf{s}\mathbf{i}\mathbf{s}}{\overbrace{Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{6}}}},(2)$

$\dot{T}=\stackrel{\mathbf{Proliferation}}{\overbrace{{\alpha }_2\cdot T\cdot \left(1-\frac{T}{K}\right)}}+\stackrel{\mathbf{Necrosis}}{\overbrace{Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{6}}\cdot PTf}}-\stackrel{\mathbf{Phagocytosis}}{\overbrace{M\cdot \frac{T}{PTf}\cdot {\mathit{F}}_{\mathit{\delta }}\cdot \left(1+Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{6}}\cdot \sigma \right)}},(3)$

$\dot{Tf}=\stackrel{\mathbf{Proliferation}}{\overbrace{{\alpha }_3\cdot \mathit{T}\mathit{f}\cdot \left(1-\frac{Tf}{\left(1+\tilde{K}\cdot Mf\cdot \xi \right)}\right)}}+\stackrel{\mathbf{Phagocytosis}}{\overbrace{M\cdot \frac{T}{PTf}\cdot {\mathit{F}}_{\mathit{\delta }}\cdot \left(1+Mf\cdot {\mathit{F}}_{\mathit{\beta }\mathbf{6}}\cdot \sigma \right)}}-\stackrel{\mathbf{Death}\mathbf{by}\mathbf{M}\mathbf{f}}{\overbrace{Mf\cdot Tf\cdot {\mathit{F}}_{\mathbf{\gamma }}}}.(4)$

Most terms of the system of the equation were built following the Guldberg and Waage Law of Mass Action, except for the bacterial growth rates for which we assumed a logistic rate. Terms with a negative sign decrease and positive terms increase the rate of change of the affected variable. The system parameters quantitatively represent the weight of each reaction in the network of interactions (Figure 1B).

Equation 1 represents the dynamics of macrophages. The arrival of $M$ is stimulated by $Mf$, and the phagocytosis of $T$ by $M$ can be enhanced by the cell interior released by $Mf$ upon necrosis. The third term corresponds to the natural death of $M$ and finally, we have the death of $M$ caused by $T$. Eq. 2 represents the macrophages that have phagocytosed. The individual terms correspond to $M$ phagocytizing $T$, the natural death of $Mf$, the death of $Mf$ caused by $T$ and $Tf$, apoptosis, and necrosis, respectively. Eq. 3 represents the dynamics of M. tuberculosis. The first term corresponds to the logistic growth of $T$, where ${\alpha }_2$ is a term composed by the proliferation and death rates and represents the reproduction rate of $T$ and $K$ is the carrying capacity of the system. The second term is the release of $T$ to the extracellular medium due to the necrosis of $Mf$, and the last one is $M$ phagocytizing $T$. The number of $T$ released during necrosis is controlled by the scaling factor $PTf$ that represents the average bacteria phagocytized by macrophages in vivo. Finally, Eq. 4 represents the bacteria that have been phagocytosed. The first term corresponds to the logistic growth of $Tf$, where ${\alpha }_3$ represents the net reproduction rate of $Tf$ and $\tilde{K}$ is the carrying capacity of the system, which is conditioned by the population of $Mf$ since $Tf$ cannot proliferate outside of it. The second term is $M$ engulfing $T$, and the last one is the death of $Tf$ caused by $Mf$.

The immune system and the bacteria show different phenotypes throughout the disease progression. Our model captures these changes along the phases of progression through the modulation of the parameters marked with an $F$ in Eq. 1. To reflect the abrupt and immediate change in macrophage phenotype, this modulation is modeled in a phenomenological way by two Hill functions (one for each phenotypic transformation) that depend on the integral or the history of bacterial load:

$f_i\left(X\right)=(P_{max}^i-P_{min}^i)+P_{min}^i\left(\frac{X^{n_i.j_i}}{X^{n_i.j_i}+K_{m^i}^{n_i.j_i}}\right),$

where $i=\left\{\mathrm{1,2}\right\}$ and $X={\displaystyle \underset{0}{\overset{t}{\int }}\left(T+Tf\right)dx}$. The constants $K_m^i$ represent the threshold values of $X$ where the system changes phases. We set K_m¹∼10⁷ and K_m²∼ 2 × 10⁸ CFU, corresponding to the total bacterial loads at times t_crit¹ = 6 and t_crit² = 27 days inferred from the experimental data (Figures 1C,2A; Table 3). The Hill coefficient, $n_i$, defines the steepness of the Hill function; the higher it is, the more abrupt its change will be. We chose n_i → ∞. $P_{min}^i$ and $P_{max}^i$ are the minimal and maximal parameter values, respectively, between the phase changes. $j$ can be either 1 or -1 depending on if the Hill function increases or decreases. As shown in Figure 1D, we can see three behaviors: up-down ($j_1$ = 1, $j_2$ = -1), up-up ($j_1$ = $j_2$ = 1), and down-down ($j_1$ = $j_2$ = -1). A big assumption of the model is that in each phase, only one type of macrophage and of bacteria phenotypes are present, that is, that the phenotypic conversion occurs for 100% of cell types.

TABLE 3

TABLE 3. In vivo data used for indirect parameter estimation.

4.3 Model Simulation

The solution of the model (Eq. 1) was approximated by numerical integration using ode15s. As the initial condition for $T$, we used the first available data, corresponding to day 1 (Table 3). For $M$, we used counts from Garcia-Romo et al. (2013) corresponding to a mouse under control conditions (without undergoing a bacterial challenge), and we assumed that the average of these counts is equivalent to the levels of macrophages in homeostasis. We set $Mf\left(0\right)=0$ and $Tf\left(0\right)=0$ to reflect that the dynamics for the phagocytic bacteria and macrophages do not change significantly within one day after the first bacterial exposure.

4.4 Parameter Optimization

We estimated the 16 model parameters following these three approaches:

1. Direct extraction of three parameters from the literature.

2. Direct estimation of five parameters using three different sets of in vitro data while fixing the parameters estimated in the previous step.

3. Indirect estimation of the remaining eight parameters by fitting the model results to two in vivo data sets simultaneously, fixing the previously estimated parameters.

4.4.1 Direct Extraction From Literature

4.4.1.1 Carrying Capacity K

The carrying capacity $K$ of the logistic growth of $T$ was obtained from Bini et al. (2014), where they measured the pulmonary bacilli loads (CFU) on different treatments on experimental models of progressive pulmonary TB in mice. From the counts corresponding to the control, we extracted the maximum value attained, and since it was the highest value of bacilli loads extracted so far from in vivo data, we decided to use it as $K$. For $Tf$, we assumed that $\tilde{K}=K$.

4.4.1.2 Average Bacteria Phagocytosed by Macrophage PTf

The average number of bacteria phagocytosed by macrophages $PTf$ was obtained from an in vitro study where monocyte-derived macrophages were put together with different strains of M. tuberculosis for 2 hours and then the phagocytosis index was counted (Montoya-Rosales et al., 2017). We selected the H37Rv strain and extracted the number. The number was a rounded average of six counts.

4.4.2 Direct Estimation Using in Vitro Data

Leaving the previously estimated parameters fixed, we estimate the nominal values of five other parameters by integrating in vitro data into three sub-models of our main model. These sub-models represent in vitro experimental conditions in which not all the variables of our complete in vivo model are present.

4.4.2.1 Bacterial Growth ${\mathit{\alpha }}_2$

To estimate the growth of $T$, ${\alpha }_2$, we modeled the logistic growth of isolated bacteria in the absence of macrophages or other host components as $\dot{T}={\alpha }_2·T·\left(1-\frac{T}{K}\right)$, which has the solution $T\left(t\right)=\frac{K·T_0·e^{{\alpha }_2·t}}{K+T_0·(e^{{\alpha }_2·t}-1)} .$ Solving for ${\alpha }_2$ retrieves ${\alpha }_2=\frac{\mathbf{ln}(T\left(t\right)·T\left(t_0\right)-K)}{t·\ln \left(T\left(t\right)·T\left(t_0\right)-K·T\left(t_0\right)\right)}$ . We extracted and normalized bacteria counts measured from an in vitro axenic culture of M. tuberculosis H37Rv from Yokobori et al. (2018) corresponding to hour 300 for $T\left(t\right)$ and hour 200 for $t_0$, and counts for $T\left(t\right)$ were subtracted with a sufficiently small $\epsilon$ to account for the small distance between $T\left(t\right)$ and $K$ at $t$. By substituting the resulting $K=100$, $T\left(t\right)=100-\epsilon$, $t_0=32.6$, and $t=4.17$, we obtain ${\alpha }_2=0.3423 \frac{1}{d} .$ For $Tf$, we assumed that ${\alpha }_3={\alpha }_2$.

4.4.2.2 Macrophage Death ${\mathit{\beta }}_1$

To estimate the natural death of $M$, ${\beta }_1$, we represent an isolated macrophage population, with nutrients to survive but without arrival of new macrophages, as $\dot{M}=-M·{\beta }_1$ which has the solution $M\left(t\right)=M_0·e^{-{\beta }_1·t}$ from which we can solve for ${\beta }_1=\frac{\mathbf{log}\left(\frac{M\left(t\right)}{M_0}\right)}{t}$ . According to Janssen et al. (2011), the half-life of a macrophage is approximately 365 days, thus ${\beta }_1=0.0019\frac{1}{d} .$ For $Mf$, we assume that ${\beta }_2={\beta }_1$.

4.4.2.3 Phagocytosis Rate $\mathit{\delta }$

To estimate $\delta$, the capacity of $M$ to phagocytose $T$, we represent an in vitro killing assay, which consists in putting bacteria in contact with macrophages under controlled conditions to determine the antimicrobial capacity of macrophages. Therefore, in this model, we have macrophages ($M$), bacteria ($T$), macrophages that have phagocytosed bacteria ($Mf$), and bacteria that have been phagocytosed ($Tf$) as our system variables:

$\dot{M}=-M·T·{\mathit{F}}_{\mathit{\delta }}-M·{\beta }_1-M·T·{\mathit{F}}_{{\mathit{\beta }}_4}, (5)$

$\dot{T}={\alpha }_2·T·\left(1-\frac{T}{K}\right)-M·T·{\mathit{F}}_{\mathit{\delta }}, (6)$

$\dot{Mf}=M·T·{\mathit{F}}_{\mathit{\delta }}-Mf·{\beta }_2-Mf·\left(T·{\mathit{F}}_{{\mathit{\beta }}_5}+Tf·{\mathit{F}}_{{\mathit{\beta }}_7}\right),(7)$

$\dot{Tf}={\alpha }_3·Tf·\left(1-\frac{Tf}{\left(1+\tilde{K}·Mf\right)}\right)+M·T·{\mathit{F}}_{\mathit{\delta }}-Mf·Tf·{\mathit{F}}_{\mathit{\gamma }}.(8)$

To find the value of $\delta$, we make some simplifications to Eq. 2, and these are justifiable by the following assumptions made from the conditions in which the killing assays were carried out (Bongiovanni et al., 2015):

• It is assumed that no parameter adaptations occur, since the in vitro assay contains no components of the adaptive immune system.

• It is assumed that there is no natural death of macrophages (${\mathit{\beta }}_1\to 0$ and ${\mathit{\beta }}_2\to 0$) since their half-life is 365 days (Janssen et al., 2011), and the killing assays were performed within 3 h. It is also assumed that there is not enough time for there to be a significant number of macrophage deaths caused by bacteria, so ${\mathit{F}}_{{\mathit{\beta }}_4}\to 0$, ${\mathit{F}}_{{\mathit{\beta }}_5}\to 0$, and ${\mathit{F}}_{{\mathit{\beta }}_7}\to 0$.

• It is assumed that there is no proliferation of bacteria (${\alpha }_2\to 0$ and ${\alpha }_3\to 0$). Bacteria divide every 16–20 h (Sharma & Singh, 2011), which is well above the 3 h that the experiment lasts.

• In Eq. 8.4, we also assume that there is no proliferation of bacteria (${\alpha }_3\to 0$) and that there is no death of bacteria by macrophages (${\mathit{F}}_{\gamma }\to 0$).

With these assumptions, Eq. 2 can be reduced to

$\dot{M}=-T·M·{\mathit{F}}_{\delta },(9)$

$\dot{T}=-T·M·{\mathit{F}}_{\delta },(10)$

$\dot{Mf}=+T·M·{\mathit{F}}_{\delta },(11)$

$\dot{Tf}=+T·M·{\mathit{F}}_{\delta }.(12)$

where Eq. 9 is uncoupled from Mf and Tf. If we additionally assume that the variable $T\gg Tf$ in this time interval, that is, that $T$ does not change significantly, then $M$ only depends on itself.

$\dot{M}=-M·\breve{T}·\delta ,(13)$

where $\breve{T}$ is a constant that represents the invariable $T$ population. Eq. 13 has the solution $M\left(t\right)=M\left(t_0\right)·e^{-{\widehat{\mathbf{F}}}_{\delta }·\breve{T}·t},$ which can be solved for ${\widehat{\mathbf{F}}}_{\delta }$ as ${\widehat{\mathbf{F}}}_{\delta }=\frac{\mathbf{log}\left(\frac{M\left(t\right)}{M_0}\right)}{\breve{T}·t}$. Substituting the experimental values for $M\left(t_{exp}\right)=0.79$, $M\left(t_0\right)=1$, $\breve{T}=25\times {10}^6$, and $t_{\mathit{exp}}=3$ reported in (Bongiovanni et al., 2015), we obtain $\delta =3.14\times {10}^{-9} \frac{1}{bac\times h} .$

4.4.3 Parameter Optimization From In Vivo Data

Leaving the previously estimated parameters fixed, we estimated the remaining eight parameters using global optimization algorithms, fitting the model simultaneously to two in vivo data sets detailed below by minimizing the cost function:

$\mathit{Cost}\left(\overline{P}\right)\mathrm{=}{\displaystyle \sum \sqrt{{\mathit{(log(Mto}{t}_{\mathit{pred}}\mathit{+1)-log(Mto}{t}_{\mathit{exp}}\mathrm{))}}^{\mathrm{2}}}}\mathrm{+}{\displaystyle \sum \sqrt{{\mathit{(log(Tto}{t}_{\mathit{pred}}\mathit{+1)-log(Tto}{t}_{\mathit{pred}}\mathrm{))}}^{\mathrm{2}}}},$

which quantifies the distance between the model ($Mto{t}_{pred}$ and $Tto{t}_{pred}$) and the experimental points ($Mto{t}_{exp}$ and $Tto{t}_{exp}$). We apply a logarithm since the experimental data ranges over several orders of magnitude and therefore small errors could otherwise generate a lot of noise.

Minimization was performed using MATLAB (version R2020a, SCR: 001622) using the GlobalSearch, a heuristic designed to find the global optimum for nonlinear problems.

For the system parameters that are modulated, their optimization was carried out three times, one for each phase of the disease. The parameter ${\widehat{\mathbf{F}}}_{\delta }$, directly estimated from in vitro data as detailed earlier, is also modulable, so we assume that the directly estimated value corresponds to phase 1, and we estimated the optimal values for the rest of the phases with the optimization algorithm.

The in vivo data sets used were a CD11c+ alveolar cell (macrophages) count from Garcia-Romo et al. (2013) and a total bacteria count from Bini et al. (2014) (Table 3); in both studies, the BALB/c mouse was put through the same bacterial challenge of 10⁶ M. tuberculosis H37Rv. We assume the same uncertainties for all the data used.

According to our previous classification (Figure 1C), days 1, 3, and 7 are considered part of phase 1, days 14, 21, and 28 of phase 2, and days 60 and 70 of phase 3.

For phase 3, there were no readings for CD11c+ alveolar cell day 60 count, so we made an interpolation with days 28 and 70 to obtain an approximate count. For bacteria, there was no reading for day 70, but we used day 60 count since this is its carrying capacity (the maximum value that we assume it can reach).

4.5 Simulation and Analysis of Genotypic Variants

4.5.1 Parameter Sampling

Using the statistical method Latin hypercube sampling in MATLAB (version R2020a, SCR: 001622), we varied the 16 nominal parameters of the model in a multidimensional uniform distribution from 0.1 to 10. We produced 500 sets of random parameters from the nominal sets of each phase (1,500 sets in total).

4.5.2 Numerical Simulations

We performed numerical integrations with the 500 sets of random parameters sampled from phase 1 parameters using the ode15s function in MATLAB (version R2020a, SCR:001622), using a $t = 100$ and initial conditions as previously described.

4.5.3 Computation of Stable Equilibrium Points

For each parameter set, we obtained all the steady-state solutions using the vpasolvefunction in MATLAB (version R2020a, SCR:001622) and collected only those solutions that had only real and strictly negative eigenvalues in their Jacobian matrix.

Data Availability Statement

The code for the model construction and analysis presented in the study is publicly available in the GitHub (SCR:002630) repository https://github.com/eliezerflores/Tuberculosis-model. The database generated in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Author Contributions

EF-G, MZ, RH-P, and ED-H conceived the project, EF-G and ED-H designed the project, EF-G and ED-H did the theoretical and the numerical analysis, RH-P and MZ supervised the biological accuracy of the model formulation and data integration, and EF-G and ED-H wrote the manuscript. All authors approved the final manuscript.

Funding

EF-G was a master’s student from Posgrado en Ciencias Biológicas, Universidad Autónoma de México (UNAM) and recieved fellowship no. 736434 from the Consejo Nacional de Ciencia y Tecnología (CONACyT). ED-H acknowledges the Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (PAPIIT) UNAM IA100119 and IA207822, CONACyT Ciencia de Frontera, project number 319600 and the College of Life Sciences of the Wissenschaftskolleg zu Berlin.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fsysb.2022.912974/full#supplementary-material

References

Aguilar L, D. D., Hanekom, M., Mata, D., Gey Van Pittius, N. C., Van Helden, P. D., Warren, R. M., et al. (2010). Mycobacterium tuberculosis Strains with the Beijing Genotype Demonstrate Variability in Virulence Associated with Transmission. Tuberculosis 90 (5), 319–325. doi:10.1016/j.tube.2010.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Alberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K., and Walter, P. (2002). Helper T Cells and Lymphocyte Activation.

Google Scholar

Arriaga, A. K., Orozco, E. H., Aguilar, L. D., Rook, G. A. W., and Pando, R. H. (2002). Immunological and Pathological Comparative Analysis between Experimental Latent Tuberculous Infection and Progressive Pulmonary Tuberculosis. Clin. Exp. Immunol. 128 (2), 229–237. doi:10.1046/j.1365-2249.2002.01832.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Barber, J., Carpenter, A., Torsey, A., Borgard, T., Namas, R. A., Vodovotz, Y., et al. (2021). Predicting Experimental Sepsis Survival with a Mathematical Model of Acute Inflammation. Front. Syst. Biol. 1, 5. doi:10.3389/FSYSB.2021.755913

CrossRef Full Text | Google Scholar

Bini, E. I., Mata Espinosa, D., Marquina Castillo, B., Barrios Payán, J., Colucci, D., Cruz, A. F., et al. (2014). The Influence of Sex Steroid Hormones in the Immunopathology of Experimental Pulmonary Tuberculosis. PLoS ONE 9 (4), e93831–10. doi:10.1371/journal.pone.0093831

PubMed Abstract | CrossRef Full Text | Google Scholar

Bocchino, M., Galati, D., Sanduzzi, A., Colizzi, V., Brunetti, E., and Mancino, G. (2005). Role of Mycobacteria-Induced Monocyte/macrophage Apoptosis in the Pathogenesis of Human Tuberculosis. Int. J. Tuberc. Lung Dis. 9 (4), 375–383.

PubMed Abstract | Google Scholar

Bongiovanni, B., Mata-Espinosa, D., D'Attilio, L., Leon-Contreras, J. C., Marquez-Velasco, R., Bottasso, O., et al. (2015). Effect of Cortisol And/or DHEA on THP1-Derived Macrophages Infected with Mycobacterium tuberculosis. Tuberculosis 95 (5), 562–569. doi:10.1016/j.tube.2015.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, S. B., Gern, B. H., Delahaye, J. L., Adams, K. N., Plumlee, C. R., Winkler, J. K., et al. (2018). Alveolar Macrophages Provide an Early Mycobacterium tuberculosis Niche and Initiate Dissemination. Cell Host Microbe 24 (3), 439–446. e4. doi:10.1016/j.chom.2018.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, J. M., and Ramakrishnan, L. (2009). The Role of the Granuloma in Expansion and Dissemination of Early Tuberculous Infection. Cell 136 (1), 37–49. doi:10.1016/j.cell.2008.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Domínguez-Hüttinger, E., Boon, N. J., Clarke, T. B., and Tanaka, R. J. (2017a). Mathematical Modeling of Streptococcus Pneumoniae Colonization, Invasive Infection and Treatment. Front. Physiol. 8 (115), 1–14. doi:10.3389/fphys.2017.00115

CrossRef Full Text | Google Scholar

Domínguez-Hüttinger, E., Christodoulides, P., Miyauchi, K., Irvine, A. D., Okada-Hatakeyama, M., Kubo, M., et al. (2017). Mathematical Modeling of Atopic Dermatitis Reveals "Double-Switch" Mechanisms Underlying 4 Common Disease Phenotypes. J. Allergy Clin. Immunol. 139, 1861–1872. doi:10.1016/j.jaci.2016.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Dowdy, D. W., Chaisson, R. E., Maartens, G., Corbett, E. L., and Dorman, S. E. (2008). Impact of Enhanced Tuberculosis Diagnosis in South Africa: A Mathematical Model of Expanded Culture and Drug Susceptibility Testing. Proc. Natl. Acad. Sci. U.S.A. 105 (32), 11293–11298. doi:10.1073/pnas.0800965105

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, W., and Elemento, O. (2015). Cancer Systems Biology: Embracing Complexity to Develop Better Anticancer Therapeutic Strategies. Oncogene 34 (25), 3215–3225. doi:10.1038/onc.2014.291

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, S., Butler, J. R., Mattila, J. T., and Kirschner, D. E. (2020). Systems Biology Predicts that Fibrosis in Tuberculous Granulomas May Arise through Macrophage-To-Myofibroblast Transformation. PLoS Comput. Biol. 16 (12), e1008520. doi:10.1371/JOURNAL.PCBI.1008520

PubMed Abstract | CrossRef Full Text | Google Scholar

Furman, D., Campisi, J., Verdin, E., Carrera-Bastos, P., Targ, S., Franceschi, C., et al. (20192019). Chronic Inflammation in the Etiology of Disease across the Life Span. Nat. Med. 2525 (1212), 1822–1832. doi:10.1038/s41591-019-0675-0

CrossRef Full Text | Google Scholar

Garcia-Romo, G. S., Pedroza-Gonzalez, A., Aguilar-Leon, D., Orozco-Estevez, H., Lambrecht, B. N., Estrada-Garcia, I., et al. (2004). Airways Infection with Virulent Mycobacterium tuberculosis Delays the Influx of Dendritic Cells and the Expression of Costimulatory Molecules in Mediastinal Lymph Nodes. Immunology 112 (4), 661–668. doi:10.1046/j.1365-2567.2004.01904.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia-Romo, G. S., Pedroza-Gonzalez, A., Lambrecht, B. N., Aguilar-Leon, D., Estrada-Garcia, I., Hernandez-Pando, R., et al. (2013). Mycobacterium tuberculosis Manipulates Pulmonary APCs Subverting Early Immune Responses. Immunobiology 218 (3), 393–401. doi:10.1016/j.imbio.2012.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Glynn, J. R., Whiteley, J., Bifani, P. J., Kremer, K., and Van Soolingen, D. (2002). Worldwide Occurrence of Beijing/W Strains of Mycobacterium tuberculosis: A Systematic Review. Emerg. Infect. Dis. 8 (8), 843–849. doi:10.3201/eid0805.020002

PubMed Abstract | CrossRef Full Text | Google Scholar

Golstein, P., and Kroemer, G. (2007). Cell Death by Necrosis: towards a Molecular Definition. Trends Biochem. Sci. 32, 37–43. doi:10.1016/j.tibs.2006.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernández-Pando, R., Marquina-Castillo, B., Barrios-Payán, J., and Mata-Espinosa, D. (2012). Use of Mouse Models to Study the Variability in Virulence Associated with Specific Genotypic Lineages of Mycobacterium tuberculosis. Infect. Genet. Evol. 12 (4), 725–731. doi:10.1016/j.meegid.2012.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernandez-Pando, R., Orozco, H., Honour, J., Silva, P., Leyva, R., and Rook, G. A. W. (1995). Adrenal Changes in Murine Pulmonary Tuberculosis; a Clue to Pathogenesis? FEMS Immunol. Med. Microbiol. 12 (1), 63–72. doi:10.1111/J.1574-695X.1995.TB00176.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernández-Pando, R., Orozcoe, H., Sampieri, A., Pavón, L., Velasquillo, C., Larriva-Sahd, J., et al. (1996). Correlation between the Kinetics of Th1, Th2 Cells and Pathology in a Murine Model of Experimental Pulmonary Tuberculosis. Immunology 89 (1), 26–33.

PubMed Abstract | Google Scholar

Hernández-Pando, R., Schön, T., Orozco, E. H., Serafin, J., and Estrada-García, I. (2001). Expression of Inducible Nitric Oxide Synthase and Nitrotyrosineduring the Evolution of Experimental Pulmonary Tuberculosis. Exp. Toxicol. Pathology 53 (4), 257–265. doi:10.1078/0940-2993-00182

CrossRef Full Text | Google Scholar

Janssen, W. J., Barthel, L., Muldrow, A., Oberley-Deegan, R. E., Kearns, M. T., Jakubzick, C., et al. (2011). Fas Determines Differential Fates of Resident and Recruited Macrophages during Resolution of Acute Lung Injury. Am. J. Respir. Crit. Care Med. 184 (5), 547–560. doi:10.1164/rccm.201011-1891OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Ley, K. (2017). M1 Means Kill; M2 Means Heal. J. I. 199 (7), 2191–2193. doi:10.4049/jimmunol.1701135

CrossRef Full Text | Google Scholar

López, B., Aguilar, D., Orozco, H., Burger, M., Espitia, C., Ritacco, V., et al. (2003). A Marked Difference in Pathogenesis and Immune Response Induced by Different Mycobacterium tuberculosis Genotypes. Clin. Exp. Immunol. 133 (1), 30–37. doi:10.1046/j.1365-2249.2003.02171.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Loscalzo, J., and Barabasi, A. L. (2011). Systems Biology and the Future of Medicine. WIREs Mech. Dis. 3 (6), 619–627. doi:10.1002/wsbm.144

CrossRef Full Text | Google Scholar

Magombedze, G., Garira, W., and Mwenje, E. (2006). Modelling the Human Immune Response Mechanisms to mycobacterium Tuberculosis Infection in the Lungs. Math. Biosci. Eng. MBE 3 (4), 661–682. doi:10.3934/MBE.2006.3.661

PubMed Abstract | CrossRef Full Text | Google Scholar

Marino, S., and Kirschner, D. (2016). A Multi-Compartment Hybrid Computational Model Predicts Key Roles for Dendritic Cells in Tuberculosis Infection. Computation 4 (4), 39. doi:10.3390/COMPUTATION4040039

PubMed Abstract | CrossRef Full Text | Google Scholar

Marquina-Castillo, B., García-García, L., Ponce-de-León, A., Jimenez-Corona, M.-E., Bobadilla-Del Valle, M., Cano-Arellano, B., et al. (2009). Virulence, Immunopathology and Transmissibility of Selected Strains ofMycobacterium Tuberculosisin a Murine Model. Immunology 128 (1), 123–133. doi:10.1111/j.1365-2567.2008.03004.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Montoya-Rosales, A., Provvedi, R., Torres-Juarez, F., Enciso-Moreno, J. A., Hernandez-Pando, R., Manganelli, R., et al. (2017). lysX Gene Is Differentially Expressed Among Mycobacterium tuberculosis Strains with Different Levels of Virulence. Tuberculosis 106, 106–117. doi:10.1016/j.tube.2017.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Pai, M., Behr, M. A., Dowdy, D., Dheda, K., Divangahi, M., Boehme, C. C., et al. (2016). Tuberculosis. Nat. Rev. Dis. Prim 2, 16076. doi:10.1038/nrdp.2016.76

PubMed Abstract | CrossRef Full Text | Google Scholar

Palma, A., Jarrah, A. S., Tieri, P., Cesareni, G., and Castiglione, F. (2018). Gene Regulatory Network Modeling of Macrophage Differentiation Corroborates the Continuum Hypothesis of Polarization States. Front. Physiol. 9 (November), 1–19. doi:10.3389/fphys.2018.01659

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedroza-González, A., García-Romo, G. S., Aguilar-León, D., Calderon-Amador, J., Hurtado-Ortiz, R., Orozco-Estevez, H., et al. (2004). In Situ analysis of Lung Antigen-Presenting Cells during Murine Pulmonary Infection with Virulent Mycobacterium tuberculosis. Int. J. Exp. Pathology 85 (3), 135–145. doi:10.1111/j.0959-9673.2004.00381.x

CrossRef Full Text | Google Scholar

Pienaar, E., Matern, W. M., Linderman, J. J., Bader, J. S., and Kirschner, D. E. (2016). Multiscale Model of Mycobacterium tuberculosis Infection Maps Metabolite and Gene Perturbations to Granuloma Sterilization Predictions, Infect. Immun., 84. 1650–1669. doi:10.1128/IAI.01438-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Poon, I. K. H., Hulett, M. D., and Parish, C. R. (2010). Molecular Mechanisms of Late Apoptotic/necrotic Cell Clearance. Cell Death Differ 17, 381–397. doi:10.1038/cdd.2009.195

PubMed Abstract | CrossRef Full Text | Google Scholar

Reynolds, A., Rubin, J., Clermont, G., Day, J., Vodovotz, Y., and Bard Ermentrout, G. (2006). A Reduced Mathematical Model of the Acute Inflammatory Response: I. Derivation of Model and Analysis of Anti-inflammation. J. Theor. Biol. 242 (1), 220–236. doi:10.1016/j.jtbi.2006.02.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Reynolds, N. J. (2014). One Hundred and Twenty-Five Years and Counting: into an Era of Systems Dermatology. Br. J. Dermatol 171 (6), 1279–1281. doi:10.1111/bjd.13447

PubMed Abstract | CrossRef Full Text | Google Scholar

Ríos-Barrera, V. A., Campos-Peña, V., Aguilar-León, D., Lascurain, L. R., Meraz-Ríos, M. A., Moreno, J., et al. (2006). Macrophage and T Lymphocyte Apoptosis during Experimental Pulmonary Tuberculosis: Their Relationship to Mycobacterial Virulence. Eur. J. Immunol. 36 (2), 345–353. doi:10.1002/eji.200535202

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritz, N., Hanekom, W. A., Robins-Browne, R., Britton, W. J., and Curtis, N. (2008). Influence of BCG Vaccine Strain on the Immune Response and Protection against Tuberculosis. FEMS Microbiol. Rev. 32 (5), 821–841. doi:10.1111/J.1574-6976.2008.00118.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Roh, J. S., and Sohn, D. H. (2018). Damage-associated Molecular Patterns in Inflammatory Diseases. Immune Network 18. doi:10.4110/in.2018.18.e27

CrossRef Full Text | Google Scholar

Rohatgi, A. (2020). WebPlotDigitizer User Manual.Version 4.3.

Google Scholar

Rook, G. A. W., Hernández-Pando, R., and Zumla, A. (2009). Tuberculosis Due to High-Dose Challenge in Partially Immune Individuals: A Problem for Vaccination? J. Infect. Dis. 199 (5), 613–618. doi:10.1086/596654

PubMed Abstract | CrossRef Full Text | Google Scholar

Salvatore, P. P., and Zhang, Y. (2017). “Tuberculosis: Molecular Basis of Pathogenesis,” in Reference Module in Biomedical Sciences (Elsevier). doi:10.1016/b978-0-12-801238-3.95697-6

CrossRef Full Text | Google Scholar

Segovia-Juarez, J. L., Ganguli, S., and Kirschner, D. (2004). Identifying Control Mechanisms of Granuloma Formation during M. tuberculosis Infection Using an Agent-Based Model. J. Theor. Biol. 231 (3), 357–376. doi:10.1016/J.JTBI.2004.06.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, B., and Singh, V. (2011). Textbook of Pulmonary and Critical Care Medicine. 1st ed. New Delhi, India: Jaypee Brothers Medical Publishers (P) Ltd.

Google Scholar

Shi, C., and Pamer, E. G. (2011). Monocyte Recruitment during Infection and Inflammation. Nat Rev Immunol. 11, 762–774. doi:10.1038/nri3070

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, I. (2003). Mycobacterium tuberculosis Pathogenesis and Molecular Determinants of Virulence. Clin Microbiol Rev. 16, 463–496. doi:10.1128/CMR.16.3.463-496.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiemann, C. A., Vanlier, J., Hilbers, P. A., and van Riel, N. A. (2011). Parameter Adaptations during Phenotype Transitions in Progressive Diseases. BMC Syst. Biol. 5 (1), 174. doi:10.1186/1752-0509-5-174

PubMed Abstract | CrossRef Full Text | Google Scholar

Trauer, J. M., Denholm, J. T., and McBryde, E. S. (2014). Construction of a Mathematical Model for Tuberculosis Transmission in Highly Endemic Regions of the Asia-Pacific. J. Theor. Biol. 358, 74–84. doi:10.1016/j.jtbi.2014.05.023

PubMed Abstract | CrossRef Full Text | Google Scholar

van Riel, N. A. W., Tiemann, C. A., Vanlier, J., and Hilbers, P. A. J. (2013). Applications of Analysis of Dynamic Adaptations in Parameter Trajectories. Interface Focus. 3, 20120084. doi:10.1098/rsfs.2012.0084

PubMed Abstract | CrossRef Full Text | Google Scholar

Who, (2019). Global Tuberculosis Report. Availableat: https://www.who.int/tb/publications/global_report/tb19_Exec_Sum_12Nov2019.pdf?ua=1.

Google Scholar

Wigginton, J. E., and Kirschner, D. (2001). A Model to Predict Cell-Mediated Immune Regulatory Mechanisms during Human Infection withMycobacterium Tuberculosis. J. Immunol. 166(3), 1951–1967. doi:10.4049/JIMMUNOL.166.3.1951

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, A. J., Linas, B., Trevejo-Nuñez, G. J., Kincaid, E., Tamura, T., Takatsu, K., et al. (2007). Mycobacterium tuberculosis Infects Dendritic Cells with High Frequency and Impairs Their Function In Vivo. J. Immunol. 179 (4), 2509–2519. doi:10.4049/jimmunol.179.4.2509

PubMed Abstract | CrossRef Full Text | Google Scholar

Yokobori, N., López, B., Monteserin, J., Paul, R., Von Groll, A., Martin, A., et al. (2018). Performance of a Highly Successful Outbreak Strain of Mycobacterium tuberculosis in a Multifaceted Approach to Bacterial Fitness Assessment. Int. J. Med. Microbiol. 308 (3), 349–357. doi:10.1016/J.IJMM.2018.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, D., Stark, J., and Kirschner, D. (2008). Systems Biology of Persistent Infection: Tuberculosis as a Case Study. Nat. Rev. Microbiol. 6 (7), 520–528. doi:10.1038/nrmicro1919

PubMed Abstract | CrossRef Full Text | Google Scholar