A METAPOPULATION MODEL OF GRANULOMA FORMATION IN THE LUNG DURING INFECTION WITH MYCOBACTERIUM TUBERCULOSIS

The immune response to Mycobacterium tuberculosis infection (Mtb) is the formation of unique lesions, called granulomas. How well these granulomas form and function is a key issue that might explain why individuals experience different disease outcomes. The spatial structures of these granulomas are not well understood. In this paper, we use a metapopulation framework to develop a spatio-temporal model of the immune response to Mtb. Using this model, we are able to investigate the spatial organization of the immune response in the lungs to Mtb. We identify both host and pathogen factors that contribute to successful infection control. Additionally, we identify specific spatial interactions and mechanisms important for successful granuloma formation. These results can be further studied in the experimental setting.

1. Introduction.Tuberculosis (TB) is an aerosol-transmitted infectious disease that accounts for approximately 1.5 million deaths per year.It has been estimated that there were 8.4 million new cases of TB in 1999 [34] and that currently one-third of the world's population is infected with Mycobacterium tuberculosis (Mtb) [9].
An important and striking aspect of tuberculosis, however, is that exposure to the bacteria rarely leads to active disease.In a proportion of exposed individuals, innate immunity is effective in clearing the pathogen.Of those individuals unable to clear the bacteria, only a small proportion progress to active disease (∼ 5%) [4].In the vast majority (∼ 95%), adaptive immunity succeeds in containing the pathogen via the formation of lesions called granulomas, resulting in latent infection.Approximately 5% to 10% of individuals who initially achieve latency later develop active tuberculosis as a result of reactivation.Reactivation can occur if the immune system is compromised in some way, for example, through infection with human immunodeficiency virus (HIV-1), by aging, or by drug and alcohol abuse [8].A key aspect of tuberculosis is understanding differences in the host immune response that 536 A METAPOPULATION MODEL OF GRANULOMA FORMATION contribute to these different disease trajectories: infected individuals who achieve and maintain latency versus those who progress to active disease, either directly or from reactivation of latent infection.In particular, formation and maintenance of a unique immunological structure, called a granuloma, is thought to play a central role in the ability of the host immune response to achieve and maintain latency.
Here we develop a spatio-temporal model of the adaptive immune response to Mtb infection that allows us to simulate the dynamics of granuloma formation.
It has been conjectured that how well these granulomas form, how well they function, and how well they are maintained likely play a central role in determining which disease trajectory an infected individual follows [5,7].Some evidence suggests that smaller solid granulomas (< 3 mm in size) form in patients whose immune system controls infection.Larger, necrotic granulomas (> 5 mm in size) are less likely to effectively contain bacterial growth and spread [28].Thus, granuloma size and structure may correspond and contribute to infection outcome.Some researchers have begun to investigate the spatial distribution of immune cells within granuloma and to speculate on how their organization contributes to granuloma function [3,12].Clearly, further explication of these and other questions regarding granuloma formation is crucial to understanding the pathogenesis of tuberculosis.Granuloma formation and maintenance, however, involves a complex series of spatial and temporal interactions among a large number of components (bacteria, macrophages, chemokines, T cells, etc.) [8].It is difficult to study this integrated system experimentally, and therefore much remains unknown about specific details.We believe that mathematical modeling can give key insights into this process.
A first step in this direction was the model of Wigginton and Kirschner that examined the temporal dynamics of the host response to Mtb infection by modeling the interactions of various populations of bacteria and immune cells (e.g., macrophages and T cells) as well as effector molecules (e.g., cytokines) with a system of nonlinear ordinary differential equations (ODEs) [35].By isolating a small number of bifurcation parameters, three trajectories corresponding to different disease outcomes (clearance, latency, and active disease) were obtained.However, the host immune response to Mtb infection has an important spatial structure defined by these granulomas.
In this paper, we use a metapopulation approach (i.e., a coarse-grid discrete spatial model) to extend the original temporal model [35] to a spatio-temporal model tracking the dynamics of granuloma formation and maintenance in the lung during Mtb infection.A metapopulation approach essentially means that a spatial domain is discretized into an n × n lattice of compartments.Each compartment within the lattice contains subpopulations of each of the various cell and protein types.ODEs are used to model the local interactions of these subpopulations within each compartment as well as the migration of subpopulations between adjacent compartments.
The metapopulation framework has been used extensively in mathematical biology as a way of incorporating spatial heterogeneity into population models, usually in ecological or epidemiological settings (c.f.[18,19,25,2,13,14,16,17]). Whereas these works have been performed at the level of populations of individuals, we apply the metapopulation framework in a novel way to study within-host dynamics of host-pathogen interactions at the cellular level.
2. Mathematical model.Here, we develop a spatial version of a mathematical model of Mtb infection based on the original model of Wigginton and Kirschner [35].That model incorporated three macrophage classes, two types of bacteria, three types of T cells, and four effector molecules (i.e., cytokines).For simplicity, we reduce the number of variables in that model in a number of ways.First, we choose to remove the direct action of cytokines.Instead, we include their effects indirectly via the cells that produce them, namely, macrophages and T cells.This is reasonable, as there exists a proportionality between the two.The original model in [35] also included three subsets of T cells.Here, we collapse those to a general T cell population that incorporates the effects of all three different T cell types.Following [35], we keep three macrophage classes and two types of bacteria (discussed in sections 2.1.1 and 2.1.3).Chemokines were not directly included in [35].Since chemokines direct cell movement and enhance cell recruitment, they are important when considering spatial aspects of granuloma formation.Therefore, we introduce a new variable, C, to capture the effects of a macrophage-produced chemokine.Thus, the following types of variables are included in the present model: extracellular bacteria (B E ), intracellular bacteria (B I ), resting macrophages (M R ), infected macrophages (M I ), activated macrophages (M A ), T cells (T ), and chemokine (C).
The metapopulation framework requires discretization of the spatial domain into an n × n lattice of compartments, indexed by (i, j) for 0 ≤ i, j < n − 1 (see Fig. 1).We refer to this entire domain as the site of infection.For each compartment (i, j), there are corresponding subpopulations , and C (i,j) (t).The dynamics of each subpopulation are governed by a nonlinear ODE.Thus, the model consists of a system of 7 n 2 coupled ODEs.These ODEs capture local within-compartment interactions of these subpopulations, as well as their movement between different compartments.It is conceptually useful, however, to view the model as a collection of n 2 subsystems of ODEs, where the seven differential equations for B I (i,j) (t), B E (i,j) (t), M R (i,j) (t), . . .associated with the single compartment (i, j) comprise a subsystem.Subsystems for adjacent compartments are coupled via terms that represent movement between the compartments.In this model, we allow resting macrophages, activated macrophages, and T cells to move.Hence, the differential equations for M R (i,j) (t), M A (i,j) (t), and T (i,j) (t) include movement terms, which we denote by M ov w (i,j) (t) for w = M R , M A , T (described in detail in the next section).All other terms in the equations capture interactions between variables within each compartment, and so the subsystems are internally coupled as well.
2.1.Model.In this section, we discuss the interactions between the various cell types present in the model.These interactions give rise to terms in the ODEs governing the dynamics of cell and chemokine subpopulations.The forms of these terms are adapted from [35], and we present the equations below.
2.1.1.Macrophage dynamics.We assume that resting macrophages (M R (i,j) (t)) take up extracellular bacteria (B E (i,j) (t)) within their local environment, that is, within the same compartment.If these resting macrophages cannot clear their bacterial load, they become chronically infected macrophages (M I (i,j) (t)) at a maximal rate k 2 .This is represented by the infection term in equations ( 1) and (2).Chronic infection may result in the death of a macrophage [33].We assume that this process is determined by the number of intracellular bacteria (B I (i,j) (t)) relative to the maximal carrying capacity of the infected macrophage population (N M I (i,j) (t)), where N is the average maximal bacterial capacity of a macrophage).This process of chronically infected macrophages bursting because of excessive numbers of intracellular bacteria occurs at a rate k 17 and appears in the infected macrophage equation (2).Direct killing of infected macrophages occurs at a maximal rate, k 14 , by T cells inducing death of the macrophage, and is dependent on the ratio of T cells to infected macrophages [30,32].Mtb are capable of interfering with this process [1,26,15]; this is represented in the T cell lysis term in equation (2).Natural cell death of infected macrophages is assumed to occur at a rate of µ I (the death term in equation ( 2)).
On the other hand, resting macrophages can become activated at the site of infection (i.e., in the presence of bacteria) in response to T cell (T (i,j) (t)) signals at a maximal rate of k 3 (the activation term in equations ( 1) and (3)).These activated macrophages are important to the immune response because they kill extracellular bacteria (see section 2.1.3).Activated macrophages (M A (i,j) (t)) that are no longer receiving activation signals from T cells in response to bacteria may revert to the resting state at a maximal rate of µ da .Resting and activated macrophages are assumed to have natural death rates of µ R and µ A , respectively.
At the boundary, we assume that resting macrophages are recruited into the spatial domain through two separate mechanisms.First, a baseline level of resting macrophages are required in the lung to patrol tissue constantly for inhaled particles [24,31].This is captured in the constant source term, S M R in equation (1).Additionally, macrophages are recruited to the site of infection through chemokines released by other macrophages [33].We assume macrophages are recruited by a generic chemokine, C (i,j) (t), at a maximal rate s ch and this is modified by a Michaelis-Menten saturation function.
Finally, macrophage migration between adjacent compartments of the lattice is captured by a set of "movement terms," M ov w (i,j) (t) (for w = M R , M A ), that are included in the equations for both resting and activated macrophages (equations ( 1) and (3), respectively).Chronically infected macrophages are known to be less effective at a variety of immune functions.For example, infected macrophages have diminished phagocytic potential, a decreased level of antigen presentation ability (which also plays an important role in macrophage activation) [6], and produce less chemokine than do activated macrophages [22,27].Therefore, we expect that infected macrophages have a reduced ability to sense chemokine gradients and to move.Hence, we assume that infected macrophages cannot move out of their current compartment.A description of the movement terms, M ov w (i,j) (t), is given section 2.1.5.
These assumptions give rise to the following equations for M R (i,j) (t), M I (i,j) (t), and M A (i,j) (t).Henceforth, for simplicity, we suppress the (t) notation: (1) where f (B B I (i,j) +N M I (i,j) ; and ( In equation (1), δ i and δ j define the compartments in the spatial domain into which resting macrophage recruitment occurs.As resting macrophages are assumed to enter the site of infection only at the boundaries (i = 0 or n − 1, j = 0 or n − 1), we define δ i and δ j as follows: 2.1.2.T cell dynamics.We assume T cells are recruited into the boundary compartments because of a gradient of chemokine, C.This occurs at a maximal rate s chT which is modified by a Michaelis-Menten saturation function.Also, T cells proliferate in the presence of activated macrophages at a maximal rate of α 2 and have a natural death rate of µ T .Finally, as T cells migrate in response to chemokine signalling, movement terms (M ov T (i,j) ) are included.Hence, suppressing (t) notation, the governing equations for the T cell subpopulations within the lattice are given by where δ i and δ j are as in equation ( 4) and restrict recruitment to the boundary compartments.
2.1.3.Bacterial dynamics.Extracellular bacteria becomes intracellular when their host macrophage becomes chronically infected, at a rate k 2 .The rate of extracellular bacterial growth is α 20 and these bacteria are killed by activated and resting macrophages at rates k 15 and k 18 , respectively.When a chronically infected macrophage dies naturally, it releases bacteria at a rate µ I .In addition, we assume that intracellular bacteria growth is limited to the maximal rate α 19 .As discussed above in section 2.1.1,infected macrophages burst (and die) because of excessive amounts of intracellular bacteria (i.e., greater than the average maximal capacity N ).When these macrophages die, they release their bacteria.Thus, these bacteria become extracellular at the same rate, k 17 , at which their host macrophages die.Equations ( 6) and (7) capture the dynamics of the extra-and intracellular bacteria subpopulations, respectively, where the (t) notation has been suppressed: where 2.1.4.Chemokine dynamics.In this model, we consider a generic chemotactic signalling protein, or chemokine, C (i,j) (t).We assume that chemokine is secreted by infected macrophages at a maximal rate c I and by activated macrophages at a rate c A .Chemokine decays naturally at rate Γ C .We model chemokine diffusion by two parameters: a rate of diffusion µ c , and a parameter α C corresponding to the proportion of chemokine that leaves the compartment (i, j).This leads to the following equation for C (i,j) , where the (t) notation has been suppressed: where δ R , δ L , δ D ,and δ U are defined in equation ( 9) to indicate that there is no chemokine movement into the boundary compartments from outside the lattice: 2.1.5.Movement terms.We assume movement of certain subpopulations from a given compartment (i, j) to adjacent compartments, that is, up (U ) to (i, j − 1), down (D) to (i, j + 1), left (L) to (i − 1, j), and right (R) to (i + 1, j) (see Fig. 1).The movement terms (M ov w (i,j) ) that occur in equations ( 1), (3), and (5) (i.e., for w = M R , M A , or T , respectively) are defined as ) .The coefficients δ l ( l = R, L, D, U ), as defined in equation ( 9), are included so that only the appropriate movement terms appear in the differential equations for the boundary compartments.The parameter χ w represents the rate of movement (for each motile cell type w); while α w i,j,S , α w i,j,R , α w i,j,L , α w i,j,U , and α w i,j,D are a set of "movement coefficients."The values of these movement coefficients are intended to capture the chemotactic movement of the motile cells.The movement coefficients change over time, but they always sum to 1, with α w i,j,R , α w i,j,L , α w i,j,U , α w i,j,D , and α w i,j,S representing, respectively, the percentage of cells that are moving right (R), left (L), up (U ), and down (D) out of compartment (i, j), and that remain stationary (S) within compartment (i, j).In the absence of chemokine, we assume cells move via diffusion only, in which case the percentage of cells moving up is the same as the percentage of cells moving down (and left and right), that is, α w i,j,R = α w i,j,L = α w i,j,U = α w i,j,D (= (1 − α w i,j,S )/4).To capture chemotaxis, these percentages are altered as a function of the chemokine gradient between compartment (i, j) and the neighboring compartments.For example, if the chemokine gradient is greatest between compartment (i, j) and the one above it (i, j − 1), the percentage of cells moving up (α w i,j,U ) would be adjusted proportionally so that a greater proportion of the subpopulations move in that direction.The function used for calculating the movement coefficients is described in detail in the appendix.
Note that chemokine levels, C i,j , will change continuously as the system evolves, and thus the chemotactic movement of cells should adapt accordingly to the changing chemokine environment.Thus, we recalculate the movement coefficients at every step of the numerical solver used to solve our system of differential equations.
2.1.6.Boundary conditions.At the boundary, we assume that chemokine levels outside the lattice (site of infection) are negligible and that no cells move into the lattice by way of the movement processes described here: cells enter the lattice only through recruitment terms; however, cells can leave the lattice by movement processes (defined in section 2.1.5).

2.2.
Parameter estimation and initial conditions.Tables 1 and 2.2 summarize model parameters and the references used to estimate their values.Kinetic parameters for the within-grid dynamics are taken from the ODE model of Wiggington and Kirschner [35]; these are listed in Table 1.The parameters listed in Table 2.2 arise from two new aspects of this model: chemokine and movement terms.Movement parameters were described above.As a key aspect of this work is to explore the role of spatial contributions to granuloma formation, we vary movement parameters over a wide range to examine their effect.Approximations for chemokine parameters can be estimated from experimental data, as referenced in Table 2.2.
To establish initial conditions, we first define the spatial domain.For simplicity, we study a 5 × 5 lattice of compartments.(We discuss the implications of this choice in the discussion section.)Next, we determine the spatial dimensions of each compartment (i, j).Well-formed granulomas are of radius size < 2 mm, whereas poorly formed granulomas have radii > 2 mm (see [29]).Since we assume that bacteria and infected macrophages do not move, we take each compartment to be large enough to hold a well-formed granuloma.Therefore, we assume that each compartment on the lattice is of dimension 5 mm × 5 mm.
We assume that prior to any infection there is a background level of macrophages present in alveolar tissue that patrol for foreign particles, but that there are no T cells present prior to infection.Based on estimates of the total number of alveolar macrophages in the human lung and alveolar surface area [24,31], we estimate that there are between 10 and 100 alveolar macrophages per mm 2 of alveolar tissue.Hence, each compartment should hold approximately 250 to 2, 500 macrophages when no bacteria are present.For simplicity, we assume that the center compartment (compartment (2, 2) in our 5 × 5 lattice) has a baseline level of 2 × 10 3 macrophages, that is, M R(2,2) (0) = 2 × 10 3 .Two parameters govern baseline initial conditions of M R(i,j) (0): s M R , the baseline recruitment of resting macrophages; and χ R , the rate of movement between compartments.Using a baseline value of χ R = 0.1, we take S M R = 119.47so that the initial conditions of M R(i,j) coincide with the steady-state values in the absence of infection.

3.
Results.We will focus on a metapopulation model consisting of an n×n lattice of compartments for most of this paper.For brevity, we will present results only for n = 5 but discuss the implications for larger n in the discussion section.To test the model, we studied the model for n = 1 and found results comparable to the original ODE model of [35].To capture any of the new spatial dynamics, one must consider n > 1.
Here, we give a brief outline of the stages of the immune response to Mtb as viewed through our model.Initially we start with 2,000 resting macrophages and 25 extracellular bacteria in the center compartment (i.e., compartment (2, 2)).Resting macrophages engulf bacteria and become infected.Infected macrophages produce chemokine recruiting additional resting macrophages and T cells.Chemokine must diffuse to the boundary, and then it acts to recruit immune cells and aids their directional movement towards the center compartment (where bacteria and infected macrophages reside).Newly recruited T cells activate macrophages and begin killing infected macrophages.Activated macrophages also kill extracellular bacteria.Simultaneously, intracellular bacteria replicate within their host macrophages.These infected macrophages can burst because of excessive intracellular bacteria that are then released as extracellular bacteria.How well the immune system responds to the initial infection dictates disease progression.
3.1.Model outcomes.We obtained three disease outcomes: clearance, latency, and active disease.We distinguish between these different disease outcomes according to the bacterial load within the center compartment of the lattice.We achieve Half-sat, bacteria on chronic infection 10 6

B E c15
Half-sat, MA on T cell proliferation 10 5 MA a Max., maximal.3.2).In this situation, the adaptive immune response reduces the number of bacteria to levels B E 10 −2 .The bacterial load remains at such low levels for 150 − 400 days, depending on specific parameter values (discussed in section 3.3).This causes the immune response to relax, allowing the small but nonzero bacterial numbers to grow again.This retriggers adaptive immunity, reducing the bacterial load to negligible levels.The cycle repeats with a period ranging between 1 to 2 years.In fact, this may be the most realistic scenario; once the innate immune system fails to clear bacteria initially, there is no evidence that the immune system ever clears Mtb after infection takes hold.
This pseudoclearance trajectory arises because the chemokine level decreases as the bacterial load is reduced.This not only causes decreased recruitment of immune cells but also means previously recruited immune cells are more likely to leave the center compartment.Thus, it appears that the spatial aspects of the adaptive immune response lead to this pseudoclearance.We will discuss the possible interpretations and implications of this result in section 4.
We have numerically observed bifurcations from latency to disease, and from latency to pseudoclearance.Hence, all clearance results in this section refer to pseudoclearance, unless otherwise stated.Figures 3, 4, and 5 show the results of the model for three different values of α 20 (the extracellular bacterial growth rate) that lead to clearance, latency, and disease, respectively.(Note: Fig. 3 shows the first period of a pseudoclearance trajectory.)3.2.Surface plots and spatial distributions.Figures 3, 4, and 5 show levels of various populations in the center compartment versus time only.A key aspect of this model is that we can also examine the spatial distribution of the various populations over the lattice.We visualize the spatial distributions of certain populations using surface plots.The domain of such a surface plot consists of the twodimensional 5 × 5 lattice of compartments, as shown in Figure 2. The surface plot for a population W at time t is formed by plotting the value of W (i,j) (t) as the height of the surface over compartment (i, j), where W ∈ {M R , M A , M I , B E , B I , T, C}.
For example, Figure 6 shows the spatial distributions of B E (i,j) (t), M R (i,j) (t), M I (i,j) (t), M A (i,j) (t), T (i,j) (t), and C (i,j) (t) (from top to bottom) at days t = 0, t = 150, t = 300, t = 450, and t = 600 (from right to left) in the simulation using the baseline (latency) parameter values (as given in Table 1). 1 The data for the center compartment from this simulation is given in Fig. 4, indicating that these parameter values yield latency.Note that we have not included the surface plots for B I , since they are qualitatively similar to those for M I .
This sequence of surface plots (Fig. 6) shows the spatial aspects of how infection is contained and latency is achieved in this model.The first column shows initial conditions.All variables except M R and B E are initially zero in every compartment of the lattice.For M R we have a background level of resting macrophages across the lattice; B E (2, 2) is initially 25, representing a small inoculum with Mtb.
Infection develops as resting macrophages engulf bacteria.By times t = 150 days and t = 300 days, shown in the second and third columns, infection has started to develop, while the immune response is just beginning.Substantial numbers of infected macrophages and extracellular bacteria appear in the center compartment: B E (2,2) (t) has grown from its initial value of 25 bacteria to approximately 400 at t = 150 days and 700 bacteria at t = 300 days, while M I (2,2) (t) grows from zero initially to 70 infected macrophages (at t = 150 days) and then to 500 (at t = 300 days).
Infected macrophages, however, start to secrete chemokine, which diffuses evenly from the source in the center compartment.Thus, a chemokine "cone" forms, with its peak at the center compartment.This can be seen in the surface plots for chemokine C (i,j) (t) in the last row of Figure 6.The diffusion of chemokine to the boundaries of the lattice leads to the recruitment of M R and T cells into the boundary compartments, and the chemokine gradient represented by the cone directs their movement to the center compartment.Hence, the surface plots for M R (i,j) (t) and T (i,j) (t) at t = 300 days (3rd column, 2nd and 6th rows) exhibit elevated levels in the boundary compartments and exhibit peaks in the center.
The M R peak in the center compartment collapses because resting macrophages become either infected or activated when they reach the center compartment (as seen in the transition from t = 300 days to t = 450 days in the surface plots of M R (i,j) (t)).This can be inferred from the sharp increases in the center compartment levels of M I and M A , going from t = 150 days to t = 450 days (Figure 6, 3rd and 4th rows).The M R (i,j) (t) values soon achieve the the steady-state spatial distribution shown for days t = 450 and t = 600 (2nd row, 4th and 5th columns).There continue to be elevated levels of M R (i,j) (t) at the boundaries of the lattice due to chemokine-driven recruitment, but a valley at the center appears because of the infection and activation dynamics mentioned previously.
T cells also approach a steady state with a distinct spatial pattern.As with M R , chemokine-driven recruitment leads to elevated levels of T cells at the boundaries.T cells, however, maintain a sharp peak at the center, because they move to the center and are held there by the chemokine environment.There is correspondingly a relative valley in T cell levels in the ring of "intermediate" compartments, that is, those between the boundary and the center (compartments {(i, j) : 1 ≤ i, j ≤ 3} − {(2, 2)}; see Fig. 2).This follows, because the majority of T cells that are recruited into the boundary compartments pass through these compartments quickly and then accumulate in the center compartment.
By t = 600 days (the last column of surface plots in Fig. 6), all variables have nearly reached their steady-state values: M R and T exhibit the patterns discussed above, while C levels maintain a conelike gradient.Recall that B E , M I , and B I are nonzero only in the center compartment; we showed in Figure 4 that B E (2,2) , M I (2,2) , and B I (2,2) go to steady-state values.Even though M A are allowed to move, they are heavily concentrated in the center compartment, since that is the only compartment in which resting macrophages can become activated, and since the chemokine gradient works to hold the vast majority of these activated macrophages in the center.
These surface plots demonstrate how latency is achieved in our model through the development and maintenance of a functioning granuloma.Steady-state levels of various populations within the center compartment represent the composition of a stable granuloma.Meanwhile, notice that there are nonzero steady-state values of M R , T and C outside the center compartment-both in the boundary compartments and the "intermediate ring" compartments.This highlights an important aspect of granuloma maintenance: continuous spatial recruitment of M R and T , through chemokine, is necessary to maintain the steady-state of the granuloma within the center compartment.The steady-state chemokine cone on the lattice represents the continuous production of chemokine by infected and activated macrophages in the center compartment.The diffusion of chemokine over the lattice leads to a steady influx of M R and T into the boundary compartments and then to the center compartment by chemotactically directed movement.This constant recruitment balances the internal dynamics of the elements within the center compartment: natural decay, infection and activation of M R , and so on.Because of these internal dynamics, there is a continuous turnover of the elements of the granuloma within the center compartment.However, overall a spatio-temporal equilibrium is achieved and maintained through continuous recruitment of additional immune cells.The metapopulation framework allows us to examine the spatial distributions and interactions of immune cells that lead to this equilibrium.
The metapopulation model, with its (5 × 5) lattice of compartments, introduces a delay between secretion of chemokine in the center compartment and the arrival of immune cells.This delay results from the new spatial aspects included here.It takes time for the secreted chemokine to diffuse out from the center compartment to the boundary of the lattice, and it again takes time for recruited immune cells to move there through chemotactically directed movement.This delay gives the initial infection time to develop, and so higher levels of immune cells are needed to control it.We empirically find that values of s chM = 10 3 and s chT = 10 4 suffice to control infection.
Increasing values for recruitment of M R and T (s chM and s chT ) by this amount, however, overwhelms infection and leads to clearance.To achieve a balance between infection and the immune response, we needed also to alter kinetic parameters.From the results of [35], we focus on the effects of k 14 (the rate of T cell killing of M I ) and k 15 (the rate of M A killing of B E ).We find that decreasing k 14 from 0.7 With this set of baseline values in place, we numerically investigated bifurcation dynamics for each parameter (Table 3.2).By varying each parameter individually while holding all other parameters constant at the baseline levels, we observed which parameters lead to bifurcations from latency to either active disease or to clearance in the model.Both host and bacterial factors play a key role in disease outcome.

3.3.
Numerical bifurcation dynamics in the model.We begin by noting that many of the parameters have very similar bifurcation dynamics as compared to the original model [35].Varying any of these parameters produces the same bifurcations here in the metapopulation model, although the specific values at which the bifurcations occur differs.
In this paper, we discuss parameters that have a markedly different bifurcation properties compared with the original model of [35].Interestingly, we find that it is predominantly parameters related to T cell and M I dynamics that demonstrate different behavior in the current model.This includes both cell-interaction parameters and motility parameters.In this discussion, we elaborate on the implications regarding the role of key spatial aspects of the immune response.3.2 (hence, α 20 = 0.03).

Growth rate of T cells (a 2
) and M R activation rate (k 3 ).With respect to parameters a 2 (growth rate of T cells) and k 3 (M R activation rate), we find that increasing either parameter leads to clearance.Increasing a 2 leads to an increased proliferation of T cells.This contributes to increased lysis of M I and also increased macrophage activation within the center compartment.These two mechanisms limit both the intra-and extracellular bacterial load.
But increasing a 2 also means increased T cell proliferation outside the center compartment.This effect occurs not only from increasing a 2 , but also from the increased levels of activated macrophages, M A , in the center compartment.An increase in levels of M A in the center compartment implies that more M A diffuse out to neighboring compartments, where they interact with T cells in those compartments and lead to further T cell proliferation.Furthermore, higher T cell levels oppose M A deactivation, which would otherwise tend to occur rather quickly outside the center compartment.(The significance of the spatial interactions between T cells and M A will also be important when we discuss χ A , the rate of diffusion of M A .) Thus, increasing a 2 contributes to a positive feedback loop between macrophage activation and T cell proliferation.

Chronic macrophage infection rate (k 2 )
. In [35], a bifurcation to disease was observed when k 2 was lowered from its baseline value.With the metapopulation model, we instead find a bifurcation from latency to clearance.In the metapopulation model, the greater numbers of resting macrophages on the lattice initially mean  3.2 with the exception of α 20 = 0.06.that there are higher initial levels of M I , leading to a stronger chemokine signal.The higher levels of T cells induced by the chemokine signal balances the higher rate of infection.Moreover, the higher levels of M A limit B E growth.We believe this results from a spatial representation and is a more intuitive and biologically realistic outcome.

Rate of M death due to bursting (k 2 )
. The rate at which infected macrophages burst, k 17 , leads to clearance for values of k 17 higher (k 17 > 0.61) and lower (k 17 < 0.08) than the value required for latency (k 17 = 0.1).It is not intuitively clear why this should be.In Table 4, we attempt to elucidate this phenomena.When k 17 is small, the rates of extracellular bacterial killing by activated and resting macrophages (k 18 and k 15 , respectively) are no longer bifurcation parameters (within the range of values tested).However, decreasing T cell killing of infected macrophages, k 14 , leads to active disease.Biologically, small values of k 17 imply that intracellular bacteria have a longer time to multiply.Therefore, unless the intracellular bacterial levels are controlled, disease progression will occur because of large amounts of bacteria becoming extracellular.Thus, the rate at which infected macrophages are killed is extremely important.However, when k 17 is large, infected macrophages burst more rapidly, and therefore intracellular bacteria have less time to multiply.Therefore, although decreasing k 14 can lead to latent infection, k 14 is no longer a disease bifurcation parameter in that range of values.The rate at which activated macrophages kill extracellular bacteria is now the most important factor.This follows, since intracellular bacteria are rapidly becoming extracellular bacteria.
3.3.4.Spatial parameters.The three spatial parameters χ T , χ C , and χ A (the rates of diffusion of T cells, chemokine, and activated macrophages, respectively) demonstrate similar bifurcation dynamics to each other: increasing any one of them from its baseline value eventually leads to clearance.Conversely, decreasing any one leads to disease.This indicates that spatial movement across the lattice is key to determining disease outcome, as measured by bacterial load in the center compartment.That χ T and χ C exhibit these bifurcation dynamics is strong evidence that timing of the immune response is essential to determining the outcome of infection, as we conjectured at the beginning of this section.Increasing χ C sufficiently leads to clearance, since chemokine diffuses more quickly from the center compartment to the boundary of the lattice, and so the immune system responds to infection more quickly.On the other hand, if χ C is too small, chemokine diffusion is too slow and hence the immune response is too late, allowing infection to take hold.This leads to active disease.Similarly, χ T determines how quickly T cells recruited into the boundary compartments move by way of chemotaxis to the center compartment.Max.Killing of B E by M I never never a bifurcation to latency when k14 < 0.45 b bifurcation to latency when k 18 < 9.38e −7   That χ A is also a bifurcation parameter is very interesting, because one might think that activated macrophages play a role only within the center compartment (the only place macrophages become activated), and thus their rate of diffusion is irrelevant.However, χ A is a bifurcation parameter: increasing χ A leads to clearance, and decreasing χ A leads to disease.We believe there are two mechanisms by which χ A contributes to the immune response, which may explain these bifurcations.Outside the center compartment, activated macrophages contribute to (a) T cell proliferation, and (b) chemokine production.
Since activated macrophages can move between compartments, some migrate out of the center compartment.Increasing χ A has the effect of increasing the number of activated macrophages that move out of the center compartment.We now examine  (2,2), an adjacent "ring compartment" (1, 2), and a boundary compartment (0, 2) (a sketch of these compartments is given in Figure 2).Shown are the results for three different values of χ A : χ A = 0.01 (-) leads to disease; χ A = 0.39 (• • • ) leads to clearance; χ A = 0.2 (-) leads to latency.Other parameters are at baseline values shown in Table 3.2.the spatial distributions of certain variables to explicate the two mechanisms discussed above.We note that when bacteria initially land in the center compartment, the spatial distribution of cells and chemokine are symmetric.So, examining the variables at the center compartment, an adjacent "intermediate ring" compartment (see Fig. 2) and the adjacent boundary compartment will indicate the overall spatial distribution of the variables.Hence, we look at the values of certain variables in the three compartments (2, 2), (1, 2) and (0, 2).
Figures 7 and 8 show the temporal distributions of M A , C, and T populations for the three compartments (2, 2), (1, 2) and (0, 2). Figure 7 shows the data through t = 200 days, while Figure 8 shows the data through t = 500 days.In both figures, each row corresponds to a different variable (from top to bottom: M A , C, and T ), and each column corresponds to a different compartment (from right to left: (2, 2), (1, 2), and (0, 2)).Each plot shows data for three different trajectories, corresponding to values of χ A = 0.01, 0.2, and 0.39.In the long term, these three different values of χ A lead to the three distinct outcomes: disease, latency, and  3.2.clearance, respectively.We use Figures 7 and 8 to observe effects of different values of χ A in the short term, as we suggest that these early dynamics play an important role in final disease outcome.
The M A (1,2) population (see the top row, middle column of Figs. 7 and 8) in the three simulations diverge starting at approximately t = 140 days, according to the value of χ A .The smallest value of χ A = 0.01 leads to much lower values of M A, (1,2) (up to t ≈ 300 days; see Fig. 8), since the rate of movement of activated macrophages from compartment (2, 2) to compartment (1, 2) is much smaller.Conversely, the highest value of χ A = 0.39 yields higher levels of M A (1,2) .The pattern is similar for the levels of M A (0,2) (top row, left column of Figs. 7 and 8), since the movement of activated macrophages from compartment (1, 2) to (0, 2) is the primary source of M A (0,2) .Such differences in the propagation of activated macrophages from the center compartment to the outer compartments have important consequences for the immune response to infection within the center compartment through two mechanisms.First, activated macrophages produce chemokine efficiently, so higher levels of M A in the boundary compartments lead to higher levels of chemokine in those compartments.This can be observed in the plots of C (0,2) (middle row, left column of Figs. 7 and 8).This is significant because chemokine levels in the boundary compartments determine the amount of M R and T recruitment.Increased levels of chemokine leads to increased M R and T recruitment onto the lattice, and hence to the center compartment.
Increased recruitment explains some of the divergence in the T cell levels.But there is a second mechanism by which increased M A migration away from the center compartment strengthens the immune response: T cell proliferation.The majority of T cells move across the lattice from the boundary compartments toward the center compartment.Higher levels of M A in outer compartments leads to greater proliferation of T cells in those compartments as they migrate from the boundary to the center.
In summary, our model indicates that activated macrophage migration is an essential part of the immune response.Activated macrophages that migrate away from the center compartment ultimately contribute to the immune response, by two mechanisms-chemokine production and T cell proliferation.Activated macrophage migration thus serves as a way of carrying the "signal" of infection to the periphery.4. Discussion.In this paper, we have developed a spatio-temporal model of the immune response to Mycobacterium tuberculosis infection.Using a metapopulation framework, we have extended the temporal model of Wigginton and Kirschner [35] to explore spatial aspects of the immune response, namely granuloma formation.Our model uses a discretized spatial domain consisting of an n × n lattice of compartments, where each compartment in the lattice corresponds to an area of tissue.Within-compartment dynamics are governed by ODEs and movement between compartments are based on diffusion and chemotaxis mechanisms.
The metapopulation model introduces a delay in the initiation of an immune response, because the chemokine needs to diffuse to the boundary before additional immune cells are recruited into the lattice.We believe this delay leads to the novel pseudoclearance trajectories seen in section 3.As the immune response controls and nearly clears the bacterial load (B E 10 −2 ), chemokine levels decrease as activated macrophages are deactivated and infected macrophages are eliminated.This has the effect of decreasing the recruitment of immune cells, and this allows the bacterial load to rebound.This then triggers the immune response, which decreases the bacterial load.This cycle repeats periodically.This pseudoclearance could be interpreted in two ways.It could correspond to true clearance, as the bacterial load is reduced to negligible levels.Alternatively, a biological interpretation could be that a very small number of bacteria remain at the center of a granuloma, and as the adaptive immune response relaxes the bacterial load increases.This could thus be a mechanism for reactivation of latent infection.
One of the most interesting results is for the parameter governing activated macrophage movement, χ A , where we found that: increasing this parameter from the baseline level leads to clearance, whereas decreasing it leads to disease.This suggests that movement away from the center compartment, where the bacteria reside, aids the immune response.By examining the spatial distribution of immune cells and chemokine, we proposed that two mechanisms are responsible for this: (a) activated macrophage initiation of T cell proliferation and (b) chemokine production by activated macrophages outside the center compartment.This result emphasizes that movement kinetics and the spatial distribution of cells are key to successful control of bacterial growth.
A natural extension of our work in this paper would be to increase our lattice size to values n > 5.This would imply we are: (a) studying a larger spatial domain; or (b) creating a finer mesh leading to a higher resolution model of the same area of tissue.Studying a larger spatial domain probably would not lead to any more interesting dynamics, unless we incorporate cell movement between the lymph nodes and lung tissue (as in [20]).Using a higher resolution model would mean that we could no longer assume that the granuloma grows inside a single compartment.Therefore, we would have to account for movement of bacteria and infected macrophages.Also, we would need to recalculate all movement and source terms to obtain a new set of initial conditions.This is left to future work.
Marino et al. [20], [21] developed a two-compartmental extension of [35] to study the importance of spatial trafficking of dendritic cells and T cells between lung and lymph node in response to Mtb infection.Gammack et al. [10] have developed a reaction-advection-diffusion model to capture the spatial structure of the granuloma only for the innate immune response to infection.Using a so-called internal states formulation, they tracked the position and density of bacteria throughout the forming granuloma.Approaching this problem from a more stochastic and discrete perspective, Segovia-Juarez et al. [29] used an agent-based model to explore granuloma formation.In an effort to understand the strengths and weaknesses of each of these mathematical approaches to the same biological problem (namely, the immune response to Mtb), we compared and contrasted the results obtained with each method and summarized them in [11].
In conclusion, our model uses a metapopulation framework to develop a spatial model of the immune response to M. tuberculosis.This framework can be easily extended to account for other cell types and cytokines known to have important effects on the progression of tuberculosis infection [35].However, as a first approximation, this model indicates the relative importance of key spatial parameters relating to chemotaxis and diffusion of immune cells and the effects of spatial distribution of cells at the site of infection.The mechanisms and interactions we have identified serve as important places for bench scientists working on this topic to explore further.

Appendix. Movement terms and chemotaxis algorithm
Recall that some of the cell types in the model (namely, resting macrophages, activated macrophages, and T cells) migrate between compartments of the lattice in a chemokine-dependent manner.As described in section 2, this is accomplished by incorporating movement terms into the ODEs for these subpopulations.These movement terms include a set of movement coefficients (see Fig. 1) that are calculated at every step of the numerical solver as a function of the chemokine gradients between a given compartment and its neighboring compartments.In this appendix, we describe the algorithm used for computing these movement coefficients.
Recall that for each compartment (i, j) and each motile cell type w, we have five movement coefficients that represent the percentage of cells that move in a given

Figure 1 .
Figure1.Sketch of the lattice framework for n = 5.Shown is a generic cell population W (i,j) , representing the subpopulation for cell type W in the compartment with coordinates (i, j).If W is one of the cell types that move, W (i,j) transfers to the four adjacent compartments, as indicated: in the up (U), down (D), left (L), and right (R) directions.Also indicated is the recruitment of cells onto the lattice via the boundary compartments (here shown from two sides only).

Figure 3 .
Figure 3. Plots showing solution curves in the center compartment of the 5×5 metapopulation lattice, for parameter values that lead to pseudoclearance (shown is the first period).Parameter values are the baseline values shown in Table3.2, with the exception of α 20 = 0.005.

Figure 4 .
Figure 4. Plots showing solution curves in the center compartment of the 5×5 metapopulation lattice, for parameter values that lead to latency.All parameter values are the baseline values shown in Table3.2(hence, α 20 = 0.03).

Figure 5 .
Figure 5. Plots showing solution curves in the center compartment of the 5×5 metapopulation lattice, for parameter values that lead to disease.Parameter values are the baseline values shown in Table3.2 with the exception of α 20 = 0.06.

Figure 6 .
Figure 6.Full model: n = 5.Surface plots of the spatial distribution of B E , M R , M I , M A , T , and C variables at times t = 0, 100, 150, 300, 450, 600 (days) during latency.The domain of each plot is the 5 × 5 lattice of compartments.(The compartments correspond to the points in the grid that make up the domain.)Parameter values are the baseline values shown in Table 3.2.

Figure 7 .
Figure 7. Plots showing the temporal distribution (up to t = 200 days) of M A , C, and T populations at the center compartment (2, 2), an adjacent "ring compartment" (1, 2), and a boundary compartment (0, 2) (a sketch of these compartments is given in Figure2).Shown are the results for three different values of χ A : χ A = 0.01 (-) leads to disease; χ A = 0.39 (• • • ) leads to clearance; χ A = 0.2 (-) leads to latency.Other parameters are at baseline values shown in Table3.2.

Figure 8 .
Figure 8. Plots showing the temporal distribution (up to t = 500 days) of M A , C, and T populations at the center compartment (2, 2), an adjacent "ring compartment" (1, 2), and a boundary compartment (0, 2) (a sketch of these compartments is given in Figure2) Shown are the results for three different values of χ A : χ A = 0.01 (-) leads to disease; χ A = 0.39 (• • • ) leads to clearance; χ A = 0.2 (-) leads to latency.Other parameters are at baseline values shown in Table3.2.

Table 2 .
(2,2)kine and movement parameters.Chemokine production estimates derived from[23]and[27]true clearance with this model (i.e., the extracellular bacterial load equals zero) only by innate immunity, in which the initial resting macrophage population kills off infection before T cells arrive on the scene.Using a simple stability analysis on equations (1)-(6), we can show this occurs only if the rate of resting-macrophage killing of extracellular bacteria, k 18 , is sufficiently large-specifically, k 18 ≥ α 20 /M ss R(2,2), where α 20 is the extracellular bacterial growth rate and M ss R (2,2) is the steady-state resting macrophage population in the center compartment.The model does not yield true clearance by adaptive immunity.It does, however, demonstrate what we define as pseudoclearance for a variety of parameters (see Table

Table 3 .
Numerical bifurcations in the metapopulation model: varying each parameter individually.pseudoclearance: cell levels approach zero at the trough of each cycle; may correspond to clearance.b innate clearance c when a 19 < 0.01 : oscillations a

Table 4 .
Bifurcation table for the bursting parameter k 17