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].

In this paper we develop a spatial mathematical model of the granuloma. The model is based on the network shown in Fig 1. After infection with Mtb, alveolar macrophages secrete MCP-1 [9], which attracts monocytes from the blood migrate to the site of infection. The monocyte differentiate into macrophages M0 that are polarized to become either M2 macrophages under stimulation by IL-4, IL-13 or GM-CSF [10-12], or M1 macrophages under stimulation by TNF-α and IFN-γ, which is regulated by IL-10 [13]. M1 macrophages produce TNF-α, IL-12 and GM-CSF [14, 15]. M2 macrophages produce IL-1β[14], IL-10 [13], IL-13 [16] and IFN-α/β[14] while infected M2 macrophages (M2) produce also TNF-α[14] but not IL-13. Dendritic cells play a critical role. After digesting external bacteria, Be, they travel to lymph nodes where they induce naive T cells to differentiate to either Th1 cells under IL-12 environment [17] or Th2 cells under IL-4 environment [18]. The newly formed Th1 cells become activated by contact with M1 macrophages and IL-12 [17], a process resisted by IL-10 [19], whereas Th2 cells become activated by contact with M2 macrophages and IL-4 [18], a process resisted by Th1 cells [18, 20, 21]. TNF-α, IL-1β and Th1-produced IFN-γ inhibits the growth of bacteria in infected macrophages M2[13, 22]. On the other hand Th2 cells secrete IL-4 and IL-13 [11, 23] which enhance the polarization of M0 toward M2 [11, 12]. Th1 cells inhibit Th2 cells [18].

Fig 1

Schematic network d(See 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 anti-inflammatory 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.

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 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-29]. The concept of M1/M2 polarization ratio, R, was introduced in [30] as the ratio of pro-inflammatory signals (STAT1, NF-κB) to anti-inflammatory signals (STAT3). 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–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.

Table 1

The variables of the model; concentration and densities are in units of g/cm3.

M0:	Interstitial macrophage density	M1:	M1 macrophage density
M2:	M2 macrophage density	M2i:	Infected M2 macrophage density
T1:	Th1 cell density	T2:	Th2 cell density
I2:	IL-2 concentration	I10:	IL-10 concentration
I12:	IL-12 concentration	I13:	IL-13 concentration
Tα:	TNF-α concentration	Iγ:	IFN-γ concentration
G:	GM-CSF concentration	u:	cell velocity
D:	dendritic cell density	Iα	IFN-α and β
T0:	naive CD4+ T cells density	D0	Interstitial dendritic cells density
Be:	extracellular bacteria	Bi	intercellular bacteria
I1β:	IL-1β concentration

Equations for macrophages (M1, M2)

We assume that there is a constant source of macrophages (M0) 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 I10, whereas the function ε2 is increasing in GM-CSF, I4 and I13. 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 I10[13], whereas affinity toward M2 differentiation is provided by GM-CSF [10] and I4, I13[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 I10[13]. Hence

Equations for T cells (T1, T2)

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 T0, 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, B)

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 MCP-1

MCP-1 is produced by M2i macrophages [9]:

Equation for Interleukins

IL-1β is produced by M2 and M2i macrophages and resisted by IFN-α[14, 38]:

IL-2 is produced by Th1 cells [15]:

IL-4 is produced by Th2 cells and M2 macrophages [23, 39], hence

IL-6 is produced by infected M2 macrophages [40, 41], hence

IL-10 is produced by M2 and M2i macrophages and by DCs [13]. IL-12 is produced by M1 and M2i macrophages and by DCs, a process resisted by IL-10 [13] and IL-13 [42]. Hence IL-10 and IL-12 satisfy the equations:

IL-13 is produced by alveolar macrophages (M2) [16] and (mostly) by Th2 cells [11], so that where λ is larger than λ.

Equation for TNF-α

TNF-α is produced by M2i and M1 macrophages, processes resisted by IL-10 [43] and IFN-α[14]:

Equations for IFN

IFN-α/β is produced by M2i and M2 macrophages, a process resisted by IL-1β[14]: and IFN-γ is produced by Th1 [15], a process resisted by IFN-α[14],

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/cm3[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 , where v is the common radial cells velocity:

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, T0, contributes to Th1 and Th2 cells not only within the granuloma, as in Eqs (5) and (6), but also through the boundary conditions: where , and .

Similarly we assume that the source of migrating macrophages M0 from the lymph nodes also contributes to M1 macrophages through the boundary condition where .

We take 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:

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/cm3 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/cm3[15]:

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 g/cm3. We assume that at the time when the granuloma was formed there were 102 Be and 103 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 when α = 1, where α is the parameter representing the immune strength. The profiles of all of the model’s variables in the first 100 days are shown in Fig 2.

Table 2

Parameters’ description and value.

Parameter	Description	Value
DM	diffusion coefficient of macrophage	8.64 × 10−7 cm2 day−1[15, 55, 56]
DD	diffusion coefficient of dendritic cell	8.64 × 10−7 cm2 day−1[15, 55, 56]
DT	diffusion coefficient of T cell	8.64 × 10−7 cm2 day−1[15, 55, 56]
DBe	diffusion coefficient of Be	1 × 10−6 cm2 day−1[15]
DBi	diffusion coefficient of Bi	8.64 × 10−7 cm2 day−1[15]
DIγ	diffusion coefficient of IFN-γ	1.08 × 10−2 cm2 day−1[15]
DTα	diffusion coefficient of TNF-α	1.29 × 10−2 cm2 day−1[15]
DP	diffusion coefficient of MCP-1	1.78 × 10−2 cm2 day−1[15]
DI1β	diffusion coefficient of IL-1β	1.08 × 10−2 cm2 day−1[15]
DI2	diffusion coefficient of IL-2	1.08 × 10−2 cm2 day−1[15]
DI4	diffusion coefficient of IL-4	1.08 × 10−2 cm2 day−1[15]
DI10	diffusion coefficient of IL-10	1.08 × 10−2 cm2 day−1[15]
DI12	diffusion coefficient of IL-12	1.08 × 10−2 cm2 day−1[15]
DI13	diffusion coefficient of IL-13	1.08 × 10−2 cm2 day−1[15]
DG	diffusion coefficient of GM-CSF	1.728 × 10−2 cm2 day−1[15]
DTα	diffusion coefficient for TNF-α	1.29 × 10−2 cm2 day−1[15]
λM0	differentiation rate of M0 to M1/M2	9.3 × 10−3/day estimated
λM1	maximal rate at which M2 is activated to become M1	6 × 10−3/day [13] & estimated
λMIr	production rate by IFN-γ	10−3/day estimated
λMTα	production rate by TNF-α	10−3/day estimated
λMG	production rate by GM-CSF	10−3/day estimated
λMI4	production rate by IL-4	10−3/day estimated
λMI13	production rate by IL-13	10−3/day estimated
λM2i	maximal rate at which M2 phagocytose (extracellular) bacteria	0.28/day [13]
λD	production rate of DCs	0.06/day [13]
λBe	burst rate of M2i	1/day [13]
λT1M1	production rate of Th1 cells by M1 macrophages and IL-12	0.23/day [13] & estimated
λTI2	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 Th1 cells	2.87 × 10−5 day−1[15] & estimated
λI12M1	production rate of IL-12 by M1 macrophages	9.64 × 10−2 day−1[15] & estimated
λI12M2i	production rate of IL-12 by M2i macrophages	9.64 × 10−3 day−1[13, 15] & estimated
λI12D	production rate of IL-12 by DCs	9.64 × 10−4 day−1[13, 15] & estimated
λTαM2i	production rate of TNF-α by M2 macrophage	1.07 × 10−3 day−1[13, 15] & estimated
λTαM1	production rate of TNF-α by M1 macrophage	1.07 × 10−1 day−1[15] & estimated
λI2T1	production rate of IL-2 by Th1 cells	4.2 × 10−4 day−1[15] & estimated
λGM1	production rate of GM-CSF by M1 macrophages	4.81 × 10−3 day−1[15] & estimated
λGM2	production rate of GM-CSF by M2 macrophages	4.81 × 10−4 day−1[15] & estimated
λBi	production rate of bacteria	0.8 day−1[13]
λI1βM2	production rate of IL-1β by M2 macrophages	7.86 × 10−4 day−1 estimated
λI1βM2i	production rate of IL-1β by M2i macrophages	7.86 × 10−3 day−1 estimated
λI6M2i	production rate of IL-6 by M2i macrophages	1.74 × 10−5 day−1 estimated
λI13T2	production rate of IL-13 by Th2 cells	2.24 × 10−4 day−1[15, 16] & estimated
λI13M2	production rate of IL-13 by macrophages	5.94 × 10−4 day−1[15, 16] & estimated
λI10M2i	production rate of IL-10 by M2i macrophages	6.67 × 10−3 day−1[13, 15] & estimated
λI10M2	production rate of IL-10 by M2 macrophages	6.67 × 10−4 day−1[15] & estimated
λI10D	production rate of IL-10 by DCs	10−4 day−1[13, 15] & estimated
λI4T2	production rate of IL-10 by Th2 cells	5.96 × 10−4 day−1[15, 57, 58] & estimated
λI4M2	production rate of IL-10 by M2 macrophages	2.38 × 10−3 day−1[15, 57, 58] & estimated
λPM2i	production rate of MCP-1 by M2i macrophages	4.86 × 10−3 day−1[15]
λIαM2	production rate of IFN-α by M2 macrophages	7.72 × 10−8 day−1 estimated
λIαM2i	production rate of IFN-α by M2i macrophages	7.72 × 10−7 day−1 estimated

Table 3

Parameters’ description and value.

Parameter	Description	Value
dM1	death rate of M1 macrophage	0.02 day−1[13]
dM2	death rate of M2 macrophage	0.008 day−1[13]
dM2i	death rate of infected M2 macrophage	0.02 day−1[13]
dT1	death rate of Th1 cell	1.97 × 10−1 day−1[15, 59, 60]
dT2	death rate of Th2 cell	1.97 × 10−1 day−1[15, 59, 60]
dD	death rate of dendritic cell	0.1 day−1[15]
dIγ	degradation rate of IFN-γ	2.16 day−1[13]
dTα	degradation rate of TNF-α	55.45 day−1[13, 15]
dI1β	degradation rate of IL-1β	6.65 day−1[15]
dI2	degradation rate of IL-2	2.376 day−1[15]
dI4	degradation rate of IL-4	50 day−1[57, 61]
dI6	degradation rate of IL-6	0.173 day−1[62]
dI10	degradation rate of IL-10	8.32 day−1[63]
dI12	degradation rate of IL-12	1.38 day−1[64]
dI13	degradation rate of IL-13	12.47 day−1[15]
dG	degradation rate of GM-CSF	4.16 day−1[15]
dP	degradation rate of MCP-1	1.73 day−1[15]
M0	source of M1 macrophages	0.5 g/cm3 estimated
T0	Source of naive T cells	0.2 g/cm3 estimated
D0	inactive DC density	5 × 10−2 g/cm3[65]
KT1	Th1 cell saturation	1 × 10−1 g/ml [66, 67]
KBe	bacterial saturation	2 × 10−11Â gcm−3[13] & estimated
KIγ	IFN-γ saturation	2 × 10−7Â gcm−3[15]
KTα	TNF-α saturation	1 × 10−6 g/cm3[15]
KG	GM-CSF saturation	1 × 10−6 g/cm3[15]
KI2	IL-2 saturation	5 × 10−7 g/cm3[15]
KI4	IL-4 saturation	2 × 10−7 g/cm3[15] & estimated
KI10	IL-10 saturation	2 × 10−7 g/cm3[15]
KI12	IL-12 saturation	1.5 × 10−5 g/cm3[15]
KI13	IL-13 saturation	2 × 10−7 g/cm3[15]
δ	small parameter for numerical purpose	10−4[13] & estimated
N	burst size	50 [68]
A	source term of M2	0.05 g/cm3[15]
KM1	M1 saturation	0.5 g/cm3[13]
KM2	M2 saturation	1 g/cm3[13]
KI1β	IL-1β saturation	10−8 g/cm3 estimated
χC	chemotaxis rate	10 g/cm3[13]

Fig 2

Simulation results over 100 days with α = 1 and initial inhalation of 102 Be and 103 Bi; and all of the units are in g/cm3.

The symbols for cells and cytokines are listed in Table 1, and all the parameters are listed in 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 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 3

“Switching time” of macrophages; all the parameters and initial conditions are as in Fig 2.

M1 macrophages become dominant over infected M2 macrophages at day 47.

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 I2 in Fig 2 is a consequence of the increase in Th1, and the decrease in I13 and I4 is a consequence of decrease in T2 and M2. Similarly, the decrease in I6 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.

If we take smaller initial values in Eq (31), the granuloma may decrease: for example, if Be(0) = 10−6 g/cm3, and Bi(0) = 10−5 g/cm3 then the radius decreases continuously from 0.01 cm to 0.0095 cm (not shown here).

As demonstrated experimentally in [48], there is great variability in the nature of non-human primate granulomas. For example, the number of bacteria may vary from 0 to 106 CFU/granuloma [48]. In Fig 2, dthe average density of B approaches 0.027, which means (see the paragraph following Eq (30)) that the number of B in the granuloma at day 100 is 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).

Fig 4

“Switching time” v.s. α; all the parameters (except α) and initial conditions are as in Fig 2.

Since switching time cannot be shorter than at least two weeks [44], the biologically range of α is 0 < α < 50.

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.

Fig 5

The total bacteria load (in g/cm3) v.s. α; all the parameters (except α) and initial conditions are the same as in Fig 2.

The bacterial load decreases as the immune strength parameter α increases.

Fig 6 profiles the radius of the granuloma at day 100, as a function of the parameter α in the biologically relevant range, where 0 < α < 50. As the immune response strength increases (i.e., as α increases) more immune cells arrive to, and will remain in, the granuloma, so R(100) will increase. While with increased α the total number of bacteria decreases (see Fig 5), the granuloma radius actually increases (see Fig 6), and hence there are increased inflammatory signals, and the resulting inflammation may damage the host tissue. We conclude that a very strong immune response to infection with Mtb is not necessarily beneficial.

Fig 6

The radius of the granuloma at day 100 v.s. the immune strength parameter α; all the parameters (except α) and initial conditions are as in Fig 2.

The radius at day 10 is monotonically increasing with respect to the parameter α.

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 λ, λ and λ are reduced by a factor 1/10 during the 30 days of drug treatment. We sought to determine the optimal time to begin treatment in order to maximize the reduction in total bacterial load (B and B).

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 7

The total bacteria under IL-10 Ab treatment for α = 1; the treatment begins at different weeks and goes on for 30 days.

The IL-10 Ab is assumed to reduce the production of IL-10 by 90%, i.e., the parameters λ, λ in Eq (14) and λ are reduced by factor 1/10 during the 30 days of the treatment. All the other parameter and initial values are the same as in Fig 2.

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, α = 10(1+i) (i = 0, ⋯, 5), and divided the treatment starting-week axis by 7 equidistant points, i.e., w = j (j = 2, ⋯, 8). For each pair (α, w), we computed the total bacteria, B(α, w), after 100 days, and determined the total bacteria load reduction percentage by the formula where B0 is the total bacteria concentration without treatment. The legend of Fig 8 scales the percentage of the total bacteria reduction (B and B) with anti IL-10 Ab.

Fig 8

The percentage of total bacterial load reduction with anti IL-10 when the Ab is administered at different initial times (for 30 days) to patients with different immune strength α.

The week when treatment begins is noted on the vertical axis. The color column indicates the percentage of bacterial reduction after 100 days from the initial time when the granuloma was formed. The parameter λ, λ and λ in Eq (14) are reduced by factor of 1/10, and all the other parameters (except α) and initial values are the same as in Fig 2.

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).

Fig 9 simulates the percent bacteria reduction using anti IL-13 Ab. The results are similar to those with anti IL-10 Ab.

Fig 9

The percentage of total bacteria load reduction with anti IL-13 when the Ab is administered at different initial times (for 30 days) to patients with different ‘immune strength’ α.

The parameter λ and λ in Eq (16) are reduced by factor of 1/10, and all the other parameters (except α) and initial values are the same as in Fig 2.

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 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 pro-inflammatory 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 M0 → M1 is larger than the rate M2 → M1, taking λ larger than λ. 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/cm3, and T0 = 0.2 g/cm3. The estimation of λ and λ and the fitting of the parameter λ 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 by the same factor for all of these components, i.e., we shall sample rather than each components of separately. We have 12 equations for the cytokines, Eqs (9)–(20), and hence 12 variables to be sampled,

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,

Fig 10

The sensitivity analysis for the cytokine production rates.

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 (which decreases B). Similarly, I2 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 I13 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 λ: λ 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 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

Let , denote numerical approximations of i-th grid point and respectively when t = nτ, where τ is the time stepsize. The discretization is derived by the explicit Euler finite difference scheme, i.e., where , , and . Then the mesh is moving by , where is solved by the velocity equation. In order to make the Euler method stable, we take .