Modeling Granulomas in Response to Infection in the Lung

Alveolar macrophages play a large role in the innate immune response of the lung. However, when these highly immune-regulatory cells are unable to eradicate pathogens, the adaptive immune system, which includes activated macrophages and lymphocytes, particularly T cells, is called upon to control the pathogens. This collection of immune cells surrounds, isolates and quarantines the pathogen, forming a small tissue structure called a granuloma for intracellular pathogens like Mycobacterium tuberculosis (Mtb). In the present work we develop a mathematical model of the dynamics of a granuloma by a system of partial differential equations. The ‘strength’ of the adaptive immune response to infection in the lung is represented by a parameter α, the flux rate by which T cells and M1 macrophages that immigrated from the lymph nodes enter into the granuloma through its boundary. The parameter α is negatively correlated with the ‘switching time’, namely, the time it takes for the number of M1 type macrophages to surpass the number of infected, M2 type alveolar macrophages. Simulations of the model show that as α increases the radius of the granuloma and bacterial load in the granuloma both decrease. The model is used to determine the efficacy of potential host-directed therapies in terms of the parameter α, suggesting that, with fixed dosing level, an infected individual with a stronger immune response will receive greater benefits in terms of reducing the bacterial load.


Introduction
The lung encounters frequent challenges from inhaled particulate and microbes, the latter include intracellular pathogens of macrophages such as M. tuberculosis (Mtb). Although this organ must effectively combat these invasions, it must do so without causing excessive inflammation. Alveolar macrophages play a large role in the innate immune response of the lung. However, when these unique, immune regulatory cells are unable to eradicate pathogens, the adaptive immune system, which includes activated macrophages and lymphocytes, particularly T cells, is called upon to control the pathogens. This collection of immune cells surrounds, isolates and quarantines the pathogen, forming a small tissue structure called a granuloma.
Granulomas are well-organized structures with immune cells at various stages of differentiation [1]. Granulomas form shortly after infection [1]. Once formed, granulomas evolve [2]. In the present work, we develop a mathematical model of granuloma evolution in the context of tuberculosis (TB).
TB is an infectious disease caused by inhalation of Mtb. The bacteria disseminate inside and outside of the lung during infection, and form granulomas in various organs of the body [3]. Granulomas are the histologic hallmark of infection with Mtb. It is estimated that one third of the human population is infected with Mtb, predominantly in latent state. Only 10% of those infected with Mtb develop active TB disease over their lifetime, while 90% of otherwise healthy humans remain in a latent state following infection [4,5]; however the percentage of those with latent infection that develop active disease increases if they are immunocompromised, e.g., HIV or diabetes [1]. Inhaled bacteria that enter the lung airway and are ingested by alveolar macrophages which are highly regulated in their immune responses, and their ability to destroy digested bacteria is limited. Alveolar macrophages as well as dendritic cells (DCs) that ingest bacteria migrate to the thoracic lymph nodes where T cells are primed. T cells that arrive at infected regions of the tissue prime macrophages to become classically activated, which helps control the bacteria. Interaction between the immune cells (macrophages, DCs and T cells) and bacteria results in granuloma formation [6]. A granuloma varies in size and activity during the course of latent and active TB, and may disappear [2,7,8].
Neutrophils are typically first responders in host defence, but reports on their mycobactericidal capacity in tissue are conflicting, perhaps because of their inherent variability and short life span [24,25]; we shall therefore not include them in the model.
More details on the various links in the network of Fig 1 will be given in the section "Mathematical model".
The signaling network within a granuloma reflects the number of bacteria in the granuloma, the ability to grow within macrophages, and the more complex efforts of the immune system to kill the bacteria while keeping inflammation under control. Thus resident infected macrophages, as well as M1 and M2 macrophages, produce all classes of cytokines discussed above, but vary in the quantity of each, creating differential balancing of the cytokines and associated signaling networks. The "pro-inflammatory" versus "anti-inflammatory" response is also seen within T cell activities. Th1 cells produce IFN-γ and TNF-α, as well as IL-2, while Th2 cells produce IL-4 and IL-13; these two populations of T cells are mutually inhibiting.
Our mathematical model focuses on the evolution of a granuloma, once formed; we do not address the early mechanisms involved in granuloma formation. The evolution of the granuloma is modeled by a system of partial differential equations (PDEs) that takes place within the granuloma, and the boundary of the granuloma varies in time. Several mathematical models of TB granulomas represented by PDEs were developed and studied in [19,26]. These models include several populations of macrophages classified according to their capacity to  Table 1 for notation). Macrophage M0, from blood monocytes, polarize into either M1 or M2 macrophages under stimulation by different cytokines. Dendritic cells attract T cells from the lymph nodes into the granuloma. When an M2 macrophge phagocytoses extracellular bacteria Be, it becomes an infected macrophage M2i. M1, M2, M2i, D, Th1 and Th2 cells produce a variety of proinflammatory and antiinflammatory cytokines, which upregulate or downregulate these cells. The growth of bacteria Bi inside macrophages M2i is inhibited by IL-1β, TNF-α and IFN-γ (modulated by IL-6). When an infected macrophage M2i bursts, its Bi are released to become extracellular bacteria. phagocytose and/or kill bacteria. The immigration of immune cells from the thoracic lymph nodes into the granuloma is represented by flux conditions at the boundary of the granuloma. Simulations of the model show how the number of bacteria that enter the granuloma affects the growth of the structure. We also simulate the effect of introducing potential host-directed therapies on the growth of the granuloma. These in-silico experiments explore how potential treatments reduce the bacteria load and limit granuloma growth, which could have important implications for the treatment of TB.
The correlation between granulomas in nonhuman primates and granulomas generated in silico was studied in [27] The in silico granuloma was derived by a mathematical model that describes, in terms of ODEs, the interaction between bacteria, immune cells, and the cytokines they secrete [27][28][29]. The concept of M1/M2 polarization ratio, R MP , was introduced in [30] as the ratio of pro-inflammatory signals (STAT1 R , NF-κB R ) to anti-inflammatory signals (STAT3 R ). This ratio was assumed to be linked to the immune response to infection with Mtb. A different approach to gauge the immune response is the "switching time", which is the time (from initial infection) that it takes for M1 macrophages to exceed the number of infected M2 alveolar macrophages. This concept was introduced in [13] and used also in [27]. As shown in [13], a switching time of 55 days, for example, can be reduced to 34 days by injection of IFN-γ, which strengthens the immune response. The mathematical model in [13] uses a one compartment ODE model to describe the infection among bacteria, immune cells and the cytokines they produce, whereas the model in [27] accounts, in more detail, for the T cell response by using a two compartment model, one in the lung and the other for lymph nodes, both described by systems of ODEs.
In the present paper, we describe the evolution of granulomas by a system of PDEs, with bacteria, immune cells and cytokines that vary in time and space within the granuloma; as does the granuloma volume. Importantly, the immune response is represented phenomenologically in boundary conditions where a parameter α quantifies the rate by which macrophages and T cells enter into the granuloma; these 'flux' boundary conditions may be viewed as the contribution to the immune response from the lymph node tissue compartment. It has been shown in previous models ( [13,26,27,31] and the references therein) that recruitment of immune cells to the lung is a key factor in granuloma outcome. Herein we express this recruitment by a parameter α. The parameter α plays a critical role in our paper, in gauging the strength of the immune system. Indeed, there is a negative correlation between the flux rate α and the switching time: As α increases the switching time decreases. The biologically relevant values of α can then be determined by the fact that the switching time is typically from a few weeks to two months [13].
Therapeutic approaches to TB have been considered in [28,29,[32][33][34] using both in silico experiments and nonhuman primate experiments. In these papers multi-scale and dynamical system approaches were developed and correlated with in vivo experiments; in particular, an anti IL-10 drug was considered. In the present paper, we use the PDE model to evaluate the efficacy of anti IL-10 and anti IL-13. As seen in Fig 1, IL-10 blocks an effective anti-TB response by the immune system. Indeed, it blocks polarization of M0 into M1 macrophages and, consequently, the production of IL-12 and TNF-α and the activation of Th1 cells. IL-10 has therefore been proposed to be a good target for TB therapy. Computational models for TB treatment with IL-10 Ab were already developed in [28,32,33] showing agreement with nonhuman primate experiments. We assume that IL-10 Ab is administered continuously (daily) for a period of 30 days, and assume, for illustration, that its effect is to reduce the production of the cytokine by 90%. We then proceed to consider with our model the following questions, applicable to personalized medicine: Does the reduction in granuloma volume depend on the immune response (i.e., on the flux rate parameter α) of the infected person, and does the reduction in granuloma volume depend on the time (relative to initial infection) when the treatment began? The answers to these questions will help determine the optimal therapeutic window for the infected individual. As seen from Fig 1, IL-13 enhances the polarization of M0 to M2 macrophages and it inhibits the activation of Th1 cells. Hence, like IL-10, IL-13 compromises the immune response to TB. We use our model, as in the case of IL-10, to explore the efficacy of treatment with IL-13 Ab.

Mathematical model
The mathematical model is based on the diagram shown in Fig 1. The variables of the model are listed in Table 1.

Equations for macrophages (M 1 , M 2 )
We assume that there is a constant source of macrophages (M 0 ) migrating into the airspace and then differentiating into either alternatively (M2) activated macrophages or into activated M1 macrophages, as they enter the alveolar space [35] depending on cues they receive from several cytokines. The equations for macrophages are given as follows: where The function ε 1 is increasing in I γ and T α and decreasing in I 10 , whereas the function ε 2 is increasing in GM-CSF, I 4 and I 13 . When interstitial macrophages enter the lung (from the blood vessels) they differentiate into either M1 or M2. I γ and T α increase the affinity to differentiate into M1 macrophages, and this is regulated by I 10 [13], whereas affinity toward M2 differentiation is provided by GM-CSF [10] and I 4 , I 13 [11,12]. These facts are accounted for by the first terms on the right-hand side of Eqs (1) and (2). The second term in these equations represents the transition from M2 to M1 induced by I γ and T α [13]. The third term in Eq (3) accounts for infection of M2 by external bacteria, and the second term on the right-hand side of Eq (3) represents the burst of internal bacteria as infected macrophages undergo necrosis [13]; N is the average number of bacteria at the time of necrosis. The last term in each of Eqs (1)-(3) is death by apoptosis.

Equations for DCs (D)
Dendritic cells are activated by contact with infected macrophages, a process resisted by I 10 [13]. Hence

Equations for T cells (T 1 , T 2 )
Resting naive T cells in the thoracic lymphoid tissue are primed while in contact with antigen presented by DCs, and develop into an intermediate stage called Th0 cells, which we denote by T 0 , The Th0 cells migrate into the lung and are activated into Th1 by contact with M1 macrophages in an IL-12 environment [17], a process down-regulated by IL-10 [13], or into Th2 cells by contact with M2 macrophages under an IL-4 environment [18]. IL-2 induces proliferation of Th1 cells [15]. Both activation and proliferation of Th1 cells are antagonized by IL-13 [36], while the activation of Th2 cells is antagonized by Th1 [18,20,21]. Hence, the equation of Th1 and Th2 cells is given as follows: Equations for bacteria (B e , B i ) The extracellular and intercellular bacteria satisfy the following equations: Extracellular bacteria emerge from bursting of infected macrophages, with N bacteria per macrophage, and they are phagocytosed by DCs and macrophages. The intracellular bacteria population increases when M2 macrophages phagocytose extracellular bacteria, but also by proliferation within infected macrophages: the proliferation is modeled by logistic growth with carrying capacity limited by the macrophage capacity to combat the bursting of cells, and this rate is inhibited by I γ , T α [13] and IL-1β [22]; inhibition by IFN-γ is modulated by IL-6 [14,37].

Equation for GM-CSF (G)
GM-CSF is secreted by macrophages [15]; it satisfies the equation Equation for u We assume uniform density of the total cell population within the granuloma, so that and take θ = 1 g/cm 3 [15]. We assume that the granuloma occupies a spherical region where the radius R(t) may change in time {r < R(t)} and assume no-flux boundary conditions at r = R(t). We take all the functions in Eqs (1)-(16) to be radially symmetric. By adding all the equations for the variables in Eq (17) and assuming that all their diffusion coefficients are equal, we get an equation for r Á u ¼ 1 r 2 @ @r ðr 2 vÞ, where v is the common radial cells velocity: 1 r 2 @ @r ðr 2 vÞ ¼ X Right À hand sides of ð1Þ À ð6Þ: ð22Þ

Boundary conditions
The dynamics of free boundary is given by the kinetic equation We assume that the source of migrating naive T cells from the lymph nodes, T 0 , contributes to Th1 and Th2 cells not only within the granuloma, as in Eqs (5) and (6), but also through the boundary conditions: where a T 1 ðD; , and a T 2 ðD; Similarly we assume that the source of migrating macrophages M 0 from the lymph nodes also contributes to M 1 macrophages through the boundary condition where a M ðI g ; T a Þ ¼ a D DþD 0 We take @I 4 @n þ a 0 I 4 ¼ 0; @I 12 @n þ a 0 I 12 ¼ 0; @P @n þ a 0 P ¼ 0; @T a @n þ a 0 T a ¼ 0 and @I g @n þ a 0 I g ¼ 0; where α 0 > 0 since some of these cytokines are consumed at the boundary. Similarly, we assume that on the boundary. The dimension of the parameters α 0 , α is 1/day. We take α 0 = 1, but α will vary depending on the strength of the immune system; α is viewed as a parameter that gauges the immune response. Finally since all other variables (including M2 and M2i) are not directly involved with the migration of immune cells from the lymph nodes, we assume no-flux boundary conditions: @X @n ¼ 0; for the remaining variables: ð28Þ

Initial conditions
We estimate that in the initial phase of granuloma formation 80% of the cells are macrophages, and most of them are M2 macrophages; at this phase the number of T cells is small (only few have arrived from the lymph nodes [44]) and so is the number of activated DCs. Accordingly, we take the density of cells in units of g/cm 3 as follows: Because of presumed similarities in the granulomas produced in TB and sarcoidosis, at least in the early stages [45,46], we assume the following concentration of cytokines in units of g/ cm 3 [15]: P ¼ 8:56 Â 10 À7 ; I 4 ¼ 1:7 Â 10 À8 ; I 12 ¼ 2:49 Â 10 À6 ; I 1 ¼ 2:47 Â 10 À8 ; I 2 ¼ 1:29 Â 10 À8 ; I 6 ¼ 2:03 Â 10 À9 ; I 10 ¼ 7:45 Â 10 À10 ; I 13 ¼ 1:63 Â 10 À8 ; T a ¼ 1:24 Â 10 À7 ; I a ¼ 8:53 Â 10 À9 ; I g ¼ 9:87 Â 10 À10 ; G ¼ 8:51 Â 10 À8 : We next assume that when the granuloma is formed, it occupies a sphere of radius 0.01 cm [47], that is, R(0) = 0.01 cm. The size of Mtb is 2-4 μm in length and 0.2-0.5 μm in width. Accordingly we take the mass of one Mtb to be 4 × 10 −13 g. Then dif one Mtb in the granuloma dis replaced by a uniform mass, then the cocentration of this mass will be 4Â10 À13 g 4 3 p10 À6 cm 3 % 1 Â 10 À7 g/cm 3 . We assume that at the time when the granuloma was formed there were 10 2 Be and 10 3 Bi. Hence, initially, the simulations, given below, do not change qualitatively if we use other initial conditions for the bacteria.

Results
We first explore the PDE system Eqs (1)- (22), boundary conditions Eqs (23)-(28), initial conditions Eqs (29)-(31) and the parameters of Tables 2 and 3  As discussed in our previous work, [13], alveolar macrophages are specialized to recognize and eliminate most pathogens without causing inflammation. However when a strong λ D production rate of DCs 0.06/day [13] λ Be burst rate of M 2i 1/day [13] λ T1M1 production rate of Th1 cells by M1 macrophages and IL-12 0.23/day [13] & estimated λ TI 2 production rate of Th1 cells by IL-2 1/day [15] λ T2 production rate of Th2 cells 0.8/day [15] λ I γ T1 production rate of IFN-γ by inflammatory response is needed to eradicate pathogens, proinflammatory macrophages M1 migrate into the lung alveoli, and they eventually become dominant over the infected alveolar macrophages M2i. The time when this occurs we defined as the "switching time" [13]. From  Fig 2, but more clearly from Fig 3, we see that, with the immune strength parameter α = 1, the switching time is around 47 days, which is in agreement with the typical switching time response established in Day et. al [13]. Thus α = 1 will be considered to be the parameter value of the immune strength for a normal healthy person. The switching time represents the balance between M1/M2i. On the other hand the parameter α represents a more comprehensive system that includes Th1 and Th2 cells and a larger network of cytokines. In this sense, the parameter α provides a more robust gauge of the immune strength and response of an individual to infection with Mtb. Fig 2 shows that with initial bacterial load as in Eq (31), the granuloma continues to grow: R increases from 0.01 to 0.022 cm in 100 days. The densities of bacteria, Be and Bi, and the densities of M1 and Th1 increase while the densities of Th2 (after 10 days) and of M2, M2i decrease. We may conclude that the initial infection (Eq (31)) leads to a productive TB infection. The increase in I γ and I 2 in Fig 2 is a consequence of the increase in Th1, and the decrease in I 13 and I 4 is a consequence of decrease in T 2 and M 2 . Similarly, the decrease in I 6 and P follows the decrease in M2i. The fact that G is monotonically increasing means that M1 in Eq (20), which is increasing (in Fig 2) in time, dominates M2 which is decreasing in time.
As demonstrated experimentally in [48], there is great variability in the nature of nonhuman primate granulomas. For example, the number of bacteria may vary from 0 to 10 6 CFU/granuloma [48]. In Fig 2, dthe average density of B i approaches 0.027, which means (see the paragraph following Eq (30)) that the number of B i in the granuloma at day 100 is 0:027 10 À7 Â Gð100Þ Gð0Þ ¼ 2:7 Â 10 5 Â Gð100Þ Gð0Þ where G(t) is the volume of the granuloma at time t. An individual with a stronger immune strength, that is, with a larger α, should have a shorter switching time, that is, the switching time should be a monotonically decreasing function of α. This is illustrated in Fig 4, which shows how the switching time decreases monotonically from 47 days to 4 days as the flux rate parameter α increases from 1 to 300. However, if the initial number of bacteria is increased or decreased, the corresponding profile of Fig 4 will slightly change (not shown here). The switching time can range from a few weeks to two months [13]. From Fig 4, we see that under the initial conditions of Eq (30), the biologically relevant values of α lie in the interval 0 < α < 50. Fig 5 shows how the total number of bacteria grows in the first 100 days for individuals with different immune strength, as defined by the parameter α. We see that the total population of bacteria decreases as α increases.
We next use the model to determine efficacy of potential host-directed therapies for reducing the bacterial load. We begin with anti IL-10 Ab which is administered continuously (daily) for 30 days. We first assume that anti IL-10 reduces the production of IL-10 by 90%, i.e., the parameters λ I 10 M2i , λ I 10 D and λ I 10 M2 are reduced by a factor 1/10 during the 30 days of drug  Table 1, and all the parameters are listed in Tables 2 and 3. doi:10.1371/journal.pone.0148738.g002 treatment. We sought to determine the optimal time to begin treatment in order to maximize the reduction in total bacterial load (B e and B i ). Fig 7 simulates the case α = 1 when anti IL-10 Ab is injected at the beginning of week j, where j = 2, 3, Á Á Á, 8. We see that there is a reduction in the bacterial load at day 100 when the Ab is administered as early as possible. Since most individuals infected with Mtb do not show symptoms for the first few weeks, what Fig 7 says is, simply, that treatment should start as soon as TB is diagnosed, for example, ideally early with skin test conversion. Fig 8 shows the contour of the total bacteria load reduction percentage for general values of α. To create the percentage, we divided the α axis by 6 equidistant points, i.e, α i = 10(1+i) (i = 0, Á Á Á, 5), and divided the treatment starting-week axis by 7 equidistant points, i.e., w j = j (j = 2, Á Á Á, 8). For each pair (α i , w j ), we computed the total bacteria, B(α i , w j ), after 100 days, and determined the total bacteria load reduction percentage by the formula where B 0 is the total bacteria concentration without treatment. The legend of Fig 8 scales the percentage of the total bacteria reduction (B i and B e ) with anti IL-10 Ab. We see that delay in treatment always reduces its benefits in terms of reduction in the bacteria load. At the same time, if one patient has a stronger immune response than another patient, and both are treated with the same delay, the benefits in terms of the percent of bacteria reduction are less for a patient with a stronger immune response.
We assumed above that anti IL-10 Ab reduces IL-10 production by 90%. Results similar to Figs 7 and 8 can be established for different percentages of reductions (not shown here).

Conclusion
A granuloma is a structure consisting of a collection of immune cells that surrounds, isolates and quarantines pathogens. Granulomas are the histologic hallmark of infection with Mtb. The dynamics of a granuloma depend not only on the number of inhaled bacteria, but also on the immune response of the host, most importantly on the immune activities of macrophages and lymphocytes.
In the present work we developed a mathematical model of the dynamics of a granuloma by a system of PDEs. Although the process of growth or size reduction of a granuloma is probably similar in all granulomas, there could be differences in granulomas in different organs of the body. There may also be differences which depend on the type of Mtb. In this paper we focused on Mtb granulomas in the lung, where M2 macrophages are used to represent alveolar macrophages [13], but the model could be modified to describe granulomas in other parts of the Since switching time cannot be shorter than at least two weeks [44], the biologically range of α is 0 < α < 50. body arising from infection with other Mtb and other intracellular pathogens of macrophages that form granulomas.
In previous work, Hao et al. [15] considered granulomas that arise in sarcoidosis (a disease with unknown etiology but suspected to be of mycobacterial origin). They developed a mathematical model that was validated by tissue data, and used this model to explore potential therapeutic agents. The mathematical model of the present paper shares some common features with the one introduced by Hao et. al. [15]. However it also has several important differences. The main new aspect is that for infection with Mtb we need to differentiate between resident alveolar macrophages (alternatively activated macrophages, or M2 macrophages) and proinflammatory macrophages (M1 macrophages, or classically activated macrophages) that arise with activation of the adaptive immune system. Furthermore, as a result of ingesting bacteria, some cytokines play different roles than they do in sarcoidosis.
The 'strength' of the adaptive immune response, in our model, is represented by a parameter α, the flux rate by which T cells and M1 macrophages that immigrated from the lymph nodes enter into the granuloma through its boundary. The parameter α is negatively correlated with the 'switching time' (ST) concept introduced in [13], namely, the time it takes for the number of M1 macrophages to surpass the number of infected alveolar (M2) macrophages. Indeed, extending the concept of ST from the ODE model in [13] to the present more comprehensive PDE model, we show that ST is a monotone decreasing function of α; the immune strength parameter α increases as the ST decreases (Fig 4).
Simulations of the model show how the bacterial load decreases and the radius of the granuloma increases as the immune strength increases (Figs 5 and 6). This implies that there is a risk associated with increased immune strength, since while the bacterial load is kept under control the inflammatory signals and the inflammation caused by excessive response of T cells may be harmful to the patients' lung tissue. The model was used to determine the efficacy of potential host-directed therapies, such as anti IL-10 Ab and anti IL-13 Ab. We show that with the same dosing level for a 30 day period, delay in the start of the treatment always reduces the benefits in terms of the residual bacterial load (Fig 7). We also show that, under the same treatment, an individual with a stronger immune strength, α, will get greater benefits (Figs 8 and 9). This suggests that dosing regimen could be adjusted based on the patient's immune system: An individual whose 'adaptive immune strength' is somehow known to be very strong will require a lower dose of the drug.
It is known [49,50] that Mtb granulomas eventually undergo characteristic calcification, which has the same density as bone, making it visible on diagnostic X-rays. Calcification is generally assumed to occur in old granulomas [51]. Since we are concerned in this paper with granuloma formation during the first few months from initial infection, we did not include calcification effects in our model. In developing the mathematical model we ignored some of the important populations of immune cells. For simplicity we did not include CD8+ T cells in the model whose role can be similar to that of Th1 cells, and B cells which also regulate immune responses. We also did not include Treg cells whose main function is to regulate T cell function in order to control inflammation. The role of Th-17 cells in human granuloma is not known [44], and they were also not included in the model. We also did not include neutrophils whose role as "double edge sword" needs further experimental evaluation [52]. Finally we did not include in our model the oxygen transport within a granuloma [53]. The present paper can be considered as a first step towards a more comprehensive and predictive model of the dynamics of granulomas in response to infection.

Parameter values
Although the etiology of sarcoidosis remains unclear, there are similarities in the immune response between TB and sarcoidosis, at least during their initial phases [45,46]. For this reason we have taken parameter values from sarcoidosis [15] whenever such parameters were not available for Mtb granulomas. We have also taken parameters from earlier work by Day et. al. [13] which dealt with Mtb, and from several experimental papers. There were still several parameters for which we were unable to find any values, and these were estimated as follows.  We assumed that the rate M 0 ! M 1 is larger than the rate M 2 ! M 1 , taking λ M 0 larger than λ M 1 . We also assumed that the source of M1 macrophage (M0) is larger than the source of naive T cells (T0), taking M0 = 0.5 g/cm 3 , and T0 = 0.2 g/cm 3 . The estimation of λ D and λ T 1 M1 and the fitting of the parameter λ M 0 were done in order to ensure that granulomas can have a steady state, i.e., that the radius of granuloma, R, which is increasing under the initial condition (31) will be decreasing if the initial bacterial load is significantly smaller.

Sensitivity analysis
In the system of Eqs (1)- (20) there appear 36 production rates which were estimated in Table 2. This number is too large for sensitivity analysis. We therefore lumped together the production rates of each cytokine. For instance, in Eq (14) we shall vary the vector l I 10 ! ¼ ðl I 10 M 2 ; l I 10 D ; l I 10 M 2i Þ by the same factor for all of these components, i.e., we shall sample y I 10 l I 10 ! for 1 2 y I 10 2  rather than each components of l I 10 ! separately. We have 12 equations for the cytokines, Eqs (9)- (20), and hence 12 variables to be sampled, y P ; y I 1b ; Á Á Á ; y G : We followed the sensitivity analysis method described in [54]. We performed Latin hypercube sampling and generated 100 samples to calculate the partial rank correlation coefficient (PRCC) and p-values with respect to • the granuloma radius, R, at day t = 100; • the total bacteria, B, at day t = 100.
The results are shown in Fig 10, where, for simplicity, we replaced the variable Eq (32) by the corresponding cytokines, P; I 1b ; Á Á Á ; G:  We see that T α and I γ are positively correlated with R and negatively correlated with B. Indeed, T α and I γ bring M1 cells into the granuloma (which increases R), and they also inhibit production of B i (which decreases B). Similarly, I 2 increases proliferation of Th1 cells and hence the production of I γ , so it is positively correlated with R and negatively correlated with B. On the other hand I 13 blocks Th1 proliferation, and is negatively correlated with R and positively correlated with B. The other components in Fig 10 can similarly be explained.
From Fig 10 we conclude that any parameter that is positively (negatively) correlated with R is at the same time, negatively (positively) correlated with B.
We carried out the same sensitivity analysis with respect to the production parameters of the cells and found the same phenomenon, with the exception of the parameter λ M 0 : λ M 0 is highly positively correlated with R and mildly positively correlated with B. This means that if more macrophages enter the granuloma, the granuloma will grow at a faster rate than the bacteria, so that the bacteria concentration will decrease.

Computational method
In order to illustrate our numerical method, we consider the following convection-diffusion equation: where F X accounts for all the 'active' terms. Since the model we consider is a free boundary problem, we employ moving mesh method to compute it. Then Eq (33) can be written in the total derivative form dXðrðtÞ; tÞ dt þ divðvÞX ¼ Dr 2 XðrðtÞ; tÞ þ F X : Let r n i , X n i denote numerical approximations of i-th grid point and Xðr n i ; tÞ respectively when t = nτ, where τ is the time stepsize. The discretization is derived by the explicit Euler finite difference scheme, i.e., where X r ¼ h 2 À1 X n iþ1 Àh 2 1 X n iÀ1 Àðh 2 1 Àh 2 À1 ÞX n i h 1 ðh 2 À1 Àh 1 h À1 Þ , X rr ¼ 2 h À1 X n iþ1 Àh 1 X n iÀ1 þðh 1 Àh À1 ÞX n i h 1 ðh 1 h À1 Àh 2 À1 Þ , h À1 ¼ r n iÀ1 À r n i and h 1 ¼ r n iþ1 À r n i . Then the mesh is moving by r nþ1 where v n i is solved by the velocity equation. In order to make the Euler method stable, we take t minfh 1 ;h À1 g 2D .

Author Contributions
Conceived and designed the experiments: WH LSS AF. Performed the experiments: WH LSS AF. Analyzed the data: WH LSS AF. Contributed reagents/materials/analysis tools: WH LSS AF. Wrote the paper: WH LSS AF.