A pressure model of immune response . . .

Mycobacterium tuberculosis (Mtb) is a widely diffused infection. However, in general, the human immune system is able to contain it. In this work, we propose a mathematical model which describes the early immune response to the Mtb infection in the lungs, also including the possible evolution of the infection in the formation of a granuloma. The model is based on coupled reaction-diffusion-transport equations with chemotaxis, which take into account the interactions among bacteria, macrophages and chemoattractant. The novelty of this approach is in the modeling of the velocity field, proportional to the gradient of the pressure developed between the cells, which makes possible to deal with a full multidimensional description and efficient numerical simulations. We perform a linear stability analysis of the model and propose a robust implicit-explicit scheme to deal with long time simulations. Both in one and two-dimensions, we find that there are threshold values in the parameters space, between a contained infection and the uncontrolled bacteria growth, and the generation of granuloma-like patterns can be observed numerically.


Mathematics Subject Classification
Tuberculosis (Tb) is a common and deadly infectious disease that is caused by Mycobacterium tuberculosis.The World Health Organization (WHO) has estimated that 1.6 million deaths resulted from Tb in 2005.Both the highest number of deaths and the highest mortality per capita are in the Africa Region [25].Tb is an aerosol-transmitted infectious disease, but, in most cases, it is a curable disease.An important and striking aspect of Tb is that exposure to the bacteria rarely leads to active disease.However, it has been estimated that currently, one third of the world's population is infected with Mtb, but, for most people, innate immunity is effective in clearing the pathogens [8], [10].Of those individuals unable to clear the bacteria, only a small proportion (≈ 5%) progress to active disease [5].In the vast majority, adaptive immunity succeeds in containing the pathogen via the formation of lesions called granulomas, resulting in latent infections [6].Formation and maintenance of a granuloma is thought to play a central role in the ability of the host-immune response to achieve and maintain latency [11].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 HIV, drugs, aging, or alcohol abuse.
The organism is transmitted primarily by the respiratory route, and although it can cause disease in most organs, pulmonary tuberculosis is most common.Initial infection occurs in alveolar region of the lungs.Once infection is present in alveoli, alveolar macrophages interact with bacteria by phagocytosis process [10].In the case of Mtb, bacteria duplicate themselves in 1 − 4 days [13], i.e. they grow very slowly.Doubling times of 24 − 96 hours have been reported for Mtb also by [24] [26], and this is striking when one considers that E. coli has a doubling time of as low as 20 minutes.The bacteria growing rate is a key factor that can determine whether the system achieved latent infection or active Tb.
The initial site of infection is in the alveoli.Once in the lungs, bacteria release a chemokine, see [21], that it is an attractant for macrophages.For this reason, macrophage movements to the site of infection occurs via both random motion (diffusion) and direct cell movement (chemotaxis).Bacteria replicate also within the macrophage and induce cytokines that initiate the inflammatory response in the lungs.Macrophages and lymphocytes migrate to the site of infection and form a granuloma [6].The function of the granuloma is to segregate the infection to prevent spread to the remainder of the lung and to other organs, as well as to concentrate the immune response directly at the site of infection.The granuloma is maintained in a persistently infected host, probably due to chronic stimulation of the immune cells, and forms the basis for a tuberculous lesion.Living bacilli have been isolated from granulomas or tubercles in the lungs of persons with clinically inactive tuberculosis, indicating that the organism can persist in a granulomatous lesion for many years [16], [19].At the most fundamental level, latent tuberculosis can be viewed as an equilibrium between host and bacillus, where the number of bacteria is stable and low, while, on the other hand, in active Tb state the number of bacteria grows exponentially [13].
In this work, we propose a simple model, which describes the innate immune response to the Tb infection.Our model can be seen as a modified diffusive prey-predator system with chemotaxis.Unlike some recent works on spatial models of infection with Mtb, see for instance [10] and [11], where the so-called "no-void" condition was used, we define a velocity field using the Darcy law, as done in [1], [7].This way, we develop our model in several space dimensions, while this was not achievable by using no-void assumption.
The paper is organized as follows: in Section 2 we outline the dynamics of bacteria, macrophages and chemoattractant.Since the velocity field, in the transport term, is given by a gradient of pressure, we indicate this model as a pressure model.Then, we write the system in a dimensionless form (Section 3).In section 4, we discuss the linear stability of the system in one-dimension, by considering homogeneous steady-state solutions and their stability under small non-homogeneous perturbations.In Section 5, we describe a numerical schemes to simulate the model.We developed explicit and implicit-explicit (IMEX) numerical schemes, the second one are useful to have a better numerical stability, necessary to obtain simulations for a long time period (months or years), see [20].Hence, we investigate numerically the stable and unstable regions in the nonlinear case.After that, we compare simulations of pressure model with a no-void model (only in one space dimension), section 7, which has been developed in our framework to mimic the more complex model proposed in [10].Finally, in section 8, we extend the model in 2-dimensions.Here, we propose simulations of the model under different conditions.In the long time simulations, we consider also the strengthening (adaptive immune response) and the weakening of the immune system (drug or alcohol abuse, immune system diseases etc.) choosing a variable killing rate of bacteria and a variable death rate of macrophages.Then, we identify a range of parameters, which correspond to a persistent pattern formation in the solutions.This pattern is given by an equilibrium state between bacteria and macrophages, which is equivalent in the mathematical model to the formation of granuloma.

The Mathematical model
In this section, we introduce and develop a mathematical model able to describe, in a simple way, the innate response of human immunity system to Mtb infection.To describe this process, we assume that a droplet of Mtb is inhaled.Once bacteria are in the lungs, they begin to duplicate and produce chemoattractant, so that macrophages, randomly distributed, are attracted toward bacteria by chemotaxis movement.When macrophages meet bacteria, they can kill bacteria with more or less effectiveness.
The simplicity of our approach, is due to the fact that we deal with bacterium-macrophage interaction as a sort of prey-predator modified system.For this reason, we use reaction-advection-diffusion equations with cou-pled chemotaxis terms, describing as these unknowns interact, move, die and proliferate.These equations are obtained by a mass balance of bacteria, macrophages and chemoattractant.The transport term is due to the the pressure of macrophages and bacteria, which "push each other".This assumption seems to be more realistic, and yields a model in several space dimensions.
We propose a second model too, using the no-void condition, which is a simplified version of the model proposed in [10].This model works only in one-dimension, and it will be useful as a comparison term with the pressure model (PM model).
In order to describe the governing equations of the system, the concentration of chemoattractant [g/cm 3 ] is indicated by C, the number of bacteria and macrophages per volume [#/cm 3 ] is indicated by B and M respectively (knowing that the mass density of a single bacterium ρ B and macrophage ρ M can be considered constant).

Bacteria
Bacteria can move via diffusion and passive transport, they reproduce and can be killed by macrophages.We consider, with B(x, t), an average behavior of all bacteria, using as reproducing parameter (α) an average rate.By a mass balance of bacteria, we obtain the following equation where ∇ • (v B B) is the advection term due to the average velocity v B , which is to be determined.D B ∆B represents the random diffusion of bacteria, see also [10], α is the replication rate and λ is the death rate due to the interaction with macrophages.The last one, λ, is a sort of efficiency indicator of immune system, for this reason we choose its value in each simulation.

Macrophages
Here, we consider only the innate immunity response, i.e. we do not include activated macrophages, which require T cells for activation.Macrophages can move via diffusion, transport and chemotaxis, and they die with a rate µ.Making a mass balance, the equation governing the dynamics of the macrophage population is given by where ∇ • (v M M ) is the advection term due to the average velocity v M to be determined, D M ∆M represents the random diffusion of macrophages, χ∇ • (M ∇C) is the chemotactic term, for a detailed description see [12].µ is the death rate of macrophages, due to the natural decay and to the interactions with bacteria.

Chemoattractant
The

Initial and boundary conditions
In lungs there is a background level of macrophages, which are present in alveolar tissue and patrol for foreign particles.The total number of alveolar macrophages in the human lung, see [11], is estimated to be between 10 and 100 alveolar macrophages per mm 2 of alveolar tissue.Furthermore, we can guess that the initial infection of bacteria can be given by even 10 bacteria or more.On the boundary of the region considered, we assume to have no-matter exchange, as first approach.It means that the normal flux is zero, indicating by n the normal versor, the boundary conditions are For the initial conditions we take C(0, x, y) = C 0 (x, y), B(0, x, y) = B 0 (x, y), M (0, x, y) = M 0 (x, y).(5)

A pressure equation for the velocity field
The key unknowns of system (1, 2, 3) are the chemoattractant C, the macrophages M , the bacteria B, and the velocity fields v B and v M .Then, we need to define the two velocities to close the system.Several constitutive equations can be formulated, but unfortunately no experimental data are available on the mechanical characteristics of cells.In absence of experimental evidence, in the following we adopt the simplest constitutive equation possible: i.e. the transport velocity is proportional to the gradient of pressure exercised by the cells among them.It comes from a momentum balance and a complete description can be found in [2], [3], [1] and [7].In this paper we adopt the description previously proposed in [1], where a cell-to-cell interaction occurs in a tumor spheroid.We adapt this concept to our problem, in which bacteria interact with macrophages.First, let us consider a single type of cell, characterized by the volume fraction of cells F .The action of the surrounding cells will be described by the scalar quantity P (F ), depending only on their local density.P (F ) represents the value of pressure given by the number of particles by unitary volume, i.e., in a fixed volume.The function P will be monotone increasing in F ; i.e.: if we have a greater number of particles, the corresponding pressure will be greater.When we multiply each particle by its single volume, we get the fraction of volume filled by one type of particle, and the corresponding pressure is proportional to the fraction volume.To make this, we suppose that there exists a threshold value of volume fraction called F ; below this threshold value, P (F ) = 0, above this value we have P (F ) > 0. It means that the threshold value is the undeformed state, and corresponds to the density of cells, such that no action is exerted on their neighbors.For larger volume ratio, cells are compressed and the gradient of pressure is greater than zero, so that a repulsive force is increasing and cells tend to move toward a region with a smaller density.Following [7], [1], we model this behavior using Darcy's law, i.e.: assuming that the velocity field is given by: where K is the motility coefficient depending on the medium in which particles move, and we assumed v M = v B = v.Now, we have to choose a functional form for P (F ).We know that the pressure increases with the fraction of volume F of particles interacting, thus we can assume that the F has a maximum at F max = 1 (the fraction of volume has to be ≤ 1).Set F = F B + F M , we have neglected the volume contribution of chemoattractant C. For the sake of simplicity, we set the threshold value F = 0, and we choose a linear form for the pressure equation, so that we get Now, F B = B and F M = M , where B and M are bacterium and macrophage fractions of volume, as we will see better in the next section.So, we can write the velocity equation

Parameters
The choice of parameters is a very hard problem to solve.For this reason we followed the values found in the literature, when available, and by rough estimations otherwise.We know that the volume of alveolar macrophage is . The first parameter to know is the replication rate of bacteria: α.Mtb duplicate in 1 − 4 days [13], [24] and [26].Then, if the duplication time is t * = 20 hours (i.e. 7, 2000 • 10 4 seconds), or t * = 96 hours (i.e. 3, 456 • 10 5 seconds), we find α 20 or α 96 respectively, given by The diffusion coefficients regarding bacteria (D B ), macrophages (D M ) and chemoattractant (D C ), are given in [17], [18], [10].The chemotactic coefficient is difficult to estimate, because it is a problem to find in literature how infection can affects the chemotactic movement of the macrophages.Following the indications of [10], we choose the simplest option of assuming that uninfected macrophages have a constant chemotactic coefficient χ.Some indications are in [14], [17], [18], [22].Regarding Γ M , Γ C , σ b and b 0 we follow indications of [10].The macrophage death coefficient µ is considered smaller than the value given in [10], we make this assumption because in our model we assume no-flux of matter from external boundary.To estimate the value of K, the motility coefficient, we follow the discussion in [1] and references therein.Finally, the coefficient describing the efficiency λ of clearing bacteria is considered as a parameter, typical of a given situation, and it is given from time to time in each simulation.Parameters are listed in table 1.
10 −8 -10 −9 ˆm2 /sec ˜Macroph.death rate estimate [10] 3 Non-dimensional system It is useful to re-write our system in a non-dimensional form, assuming to stay in a two-dimensional space.We know that x, y where we choose The derivatives are as follows: B and M are the number of bacteria in a unitary volume and the number of macrophages in a unitary volume respectively.Then, to obtain nondimensional unknowns, we put where V B and V M are the volume of a single bacterium and macrophage, respectively.The non-dimensional constants are Therefore, using equations ( 1), ( 2), ( 3) and ( 8), the non-dimensional system is given by Bτ Mτ The boundary conditions are and the initial conditions are given by the nonnegative functions A list of non-dimensional parameters is given in table 2.
Remark 1 From now on, we will use always the non dimensional system; non-dimensional variables will be directly referred to as (x, y, t), and nondimensional constants, previously denoted by "hat" ( ˆ) quantities, will be indicated without "hat", unless a different notation is explicitly introduced.So, our system in the non-dimensional form rewrites where The boundary conditions are and the initial conditions are 4 Discussion on linear stability analysis in 1-D In this section, we study the linear approximation of our system in onedimension using the non-dimensional system.As first step, we find homogeneous steady state solutions.Then, we study their stability, and, finally, we analyze what happens to the linear unstable area using the nonlinear system.
Here, we assume to have an immune system able to maintain its efficiency for a long time period.Since we have no-flux of matter from the boundaries, we choose µ = 0 to maintain constant the number of macrophages.
The homogeneous steady state is indicated by V = B, M , C .Then, using equations ( 19)-( 25) and assuming µ = 0, we obtain: Thus, we indicate the homogeneous solutions by where V (t) = V + w(t), and w(t) = (b(t), m(t), c(t)) represents the homogeneous perturbations.Let V (0) = V 0 = V + w(0).Then, equations of time evolution for homogeneous solutions are We find out If m > 0, B is decreasing exponentially; if m < 0, B increases exponentially.Also C(t) is upper limited, see eq. ( 28).
This means that bacteria are contained if M + m > α λ , otherwise (m < 0) bacteria increase exponentially.

Non-homogeneous perturbations
Now, we want to see what happen considering non-homogeneous perturbations.In this case, small perturbations are ǫ, δ and γ, such that: B = B + ǫ, M = M + δ and C = C + γ.Linearizing our system about the homogeneous steady state, we find out where ω B = b 0 ( B + b 0 ) 2 .Notice that we used the approximation (1

Linear stability
We now look for solutions of (29) in the form where n indicates the number of the mode.
Proposition 1 There exists a positive value B 0 ≈ 8 • 10 −6 , such that if B < B 0 , then the solutions of (29) are linearly stable for n ≥ 1.If B > B 0 , there exists a critical index n( B), such that the n-mode n is linearly unstable for 1 ≤ n ≤ n( B).The other modes are always stable.The critical index n( B) is monotone increasing with respect to B, and it tends to a value n 1 ≈ 22.5 as B → 1.
Proof We assume p = n 2 π 2 > 0, then we have These solutions are stable if all roots Λ of the characteristic polynomial (31) lie in the left-hand complex plane: ReΛ < 0. To this goal, we use the Routh-Hurwitz conditions [15] Here, we want to see when a 1 a 2 − a 3 > 0 is satisfied.Using equations (32), ( 33) and (34), and simplifying them, we find Taking into account parameter values in table 2, λ ≈ 10 4 , and p = n 2 π 2 (n = 1, 2...), we obtain that the (36) is greater than zero for any p, making the realistic approximation that B < 10 −2 .

a 3 > 0
To satisfy a 3 > 0, we write where Now, we search for the positivity region of a 3 .First, we find out The two roots of eq. ( 37) are Calling r( B) = β 2 2 − 4β 1 β 3 , we obtain that r( B) < 0 (figure 1), for Thus we have obtained an unstable region (a 3 < 0) defined by p + ( B) and for B > B 0 .Also we observed that for B → 1, p + → αD c = 5000, this means that n + ≈ 22.5.Thus we can conclude that for n ≥ 23 the region is stable for any B.

Numerical schemes
To solve our model we make use of Implicit-Explicit (IMEX) schemes, see [20], [4].For the sake of simplicity, we show the one-dimensional scheme, which can be easily extended in two-dimensions:

IMEX scheme
We represent the system as where We integrate explicitly H(U ), while G(U ) is a stiff term which will be integrated implicitly to avoid excessively small time steps.
A general IMEX-DIRK Runge-Kutta scheme is given by, for t = n∆t The matrices Ã = (ã ik ), where ãik = 0 for j ≥ i and A = (a ik ) are ν × ν matrices such that the resulting scheme is explicit in H and implicit in G.
The DIRK formulation requires a ik = 0 for j > i [4].Now, we write H and G in the discrete numerical form.To this goal, we develop each term in the following way (where n = 1, ..., N is the time index and j = 1, ..., J is the space index): For G, we use a three-points centered approximation 6 Numerical simulations

Two examples: winning and losing bacteria
By an initial infection of 50 bacteria, we choose 170 macrophages initially diffused in the area.As an example of the variation of bacteria concentration, we present two simulations in which bacteria prevail and are contained.
In the first simulation, we choose a low efficiency parameter λ = 2.65 • 10 4 (non-dimensional), then we obtain that macrophages are not able to contain infection.In figure 3, we can see the bacteria concentration in the first two months, and the variation of the number of bacteria in the same time period.In the second example, we choose λ = 3.5 • 10 4 , and we obtain that macrophages are able to contain infection.In figure 4 we can see bacteria concentration (up), and their number evolution in the first two months (down).

Simulations of small perturbations
Here, we consider perturbations of the following homogeneous solutions: B = 0, C = 0 and M ≈ z • α λ .We analyze three cases:z = 0.1, 1, 10.The perturbation of bacteria (ǫ) is a distribution of about 100 bacteria in the center of the domain, the other perturbations are δ = 0 and γ = ǫ/(ǫ + b 0 ).Simulations are done using the non linear system ( 19)- (25).The case M ≈ 10 α λ is on the left side of figure 5, where bacteria are cleared as expected.The case M ≈ 0.1 α λ is on the right side of of figure 5: here, we observe a bacteria growth.The third case, M ≈ α λ , is in figure 6.This case is interesting because if the perturbation ǫ is homogeneously distributed in space, we observe that the number of bacteria remains the same (on the left side of figure 6 ); instead if the perturbation is non-homogeneous, and concentrated in the center of the domain, we have a mild reduction of bacteria, shown on the right side of figure 6.This means that, in our model, the geometry of the initial distribution of bacteria can influence their evolution.

Long time behavior of linear instability
In Section 4, we observed that for B > B 0 solutions are linearly unstable.In this paragraph, we study the influence of non-linear terms on the linear instability.Using the non-linear system ( 19)-( 25), we choose as initial condi- , and two different initial conditions for bacteria: 1. B(0) = B a = Ba + ǫ(0), where B a = 10 −10 < B 0 and ǫ is a small perturbation.2. B(0) = B b = Bb + ǫ(0), where Bb = 10 −5 > B 0 .
We can see the behavior of initial perturbation ǫ(0) on the left side of figure 7. Simulating ten years evolution of these two cases, by an implicit scheme, we obtain that the perturbation of B a is more quickly diffused, while perturbation of B b (linearly unstable) does not blow up, and it remains contained.We can see both simulations on the right side of figure 7.
To the light of these results, we can observe that whereas there is linear stability, perturbation tends to flatten out to its initial value; and whereas there is linear instability, perturbation tends to hold its shape.This is due to the non linear damping effect which contrasts the linear instability effect.For these reasons, results agree with our linear analysis.

Bacterium evolution for different macrophage distribution
Here, assuming an infection of 50 bacteria, we present three simulations.In the first one, we suppose to have 170 macrophages distributed in the area, in the second one 430 and in the third one 690.In each one of these cases, we adopt different values of the efficiency parameter λ, in order to estimate the efficiency required to contain or not the infection.Remark 2 It is interesting to compare these simulations with the results obtained in section 4 with homogeneous solutions.There, M ≈ α λ , was the limit between bacteria contained or not.Here, in the case I (170 macrophages) we have α/λ c ≈ 0.0167 and MI = 0.0170, in the case II we have α/λ c ≈ 0.045 and MI = 0.043, in the case III we have α/λ c ≈ 0.075 and MI = 0.069.It is interesting to observe that, both in the homogeneous solutions and full model, we obtained the same order of magnitude for the limit value of λ.

Comparison between no-void and pressure model
In this section we present a different mathematical model describing the early immune response to Mtb.We use the same unknowns, but different assumptions: no-void between bacteria and macrophages, i.e.: the same hypothesis done in [10].The goal is to make a comparison with the previous model.The adaptive immune response to Mtb infection can lead to the formation of a granuloma [9].Following the evolution of the granuloma boundary, we can be able to predict if the infection is progressing or being contained.The formation of a granuloma is mathematically similar to that of a solid tumor: both can be considered as multicellular spheroids.
It is assumed that this spheroid is well packed (i.e. it has no-voids) and that any cell proliferation, death or, in our case, phagocytosis, contributes to a volume increase (or decrease) that causes the cells to move.It means that the model is based upon a growing spheroid, in which the volume is filled by bacteria and macrophages, the chemoattractant does not contribute to the volume of granuloma.
For simplicity, we employ a one-dimensional Cartesian geometry corresponding to the growth of a slab of bacteria and macrophages rather than a spheroid.

Equations of the Model
The equations of the model are very similar to the previous equations, in fact: where (x, t) In order to complete our model, it remains to determine the advection velocity v(x, t).To this goal, we have assumed that there are no voids within the granuloma and hence, throughout the spheroid, the condition is fulfilled, where V B (constant) is the volume of a single bacterium and V M (constant) is the volume of a single macrophage.By the no-void condition (53) we have a closed system of partial differential equations.Then, we assume that the cellular velocity, v(x, t) at the granuloma boundary, x = R(t), describes the radial growth of the granuloma and thus the equation governing R is In general, for three-dimensional (two-dimensional) spheroid growth, the equations given with no-void condition are insufficient to determine all the unknowns v, B, M and C: two (one) additional equations are needed to close the system.In this model, we avoid these issues by assuming that the spheroid undergoes only one-dimensional growth.

Boundary and Initial Conditions
The point x = 0 corresponds to the origin of granuloma, here we assume that there is no incoming flux of matter: At the granuloma boundary x = R(t), we assume that the chemoattractant is released so that macrophages are attracted into the granuloma [10].Following [10], we assume that external macrophages can go into granuloma attracted by chemokines on the boundary.Therefore the boundary conditions are given by The last boundary condition, B x , is given by the no-void applied to the external boundary.

Initial condition
At the initial time t = 0 we suppose that the spheroid is filled by bacteria and, possibly by macrophages.That is The initial radius is given by the size of initial bacteria and macrophages: R(0) = R 0 > 0.

Non-Dimensional Equations
It can be useful to write our system in dimensionless form.The non-dimensional system obtained is where v is given by the no-void condition B+M = 1, which we describe in the next section.We have used x = x dim R(t) , t = t dim ζ to obtain the non-dimensional variables (the dimensional variables have the subscript dim).
The boundary conditions are The initial conditions are Finally, the velocity v, using the no-void condition, is obtained adding the equation of bacteria (59) to the equation of macrophages (60).That is (66) Then, integrating in x, we get

Numerical schemes and model comparison
Here, we set up a numerical scheme able to deal with the no-void equations.This is a problem with a moving boundary given by R(t), for this reason we prefer solving the non-dimensional model with a fixed boundary.Explicit and IMEX schemes gives us some stability problem, mainly due to the chemotaxis and diffusion terms.To obtain a more stable system we use an implicit numerical scheme, which is stable and monotone even using large time steps, required to solve our system.The scheme is in equations ( 68) and (69), where n = 1, ..., N and j = 1, ..., J are time and space index respectively.
where v is given by equation ( 67).This way, we have to solve an algebraic system A n+1 C n+1 = p n for C and H n+1 M n+1 = q n for M .

Parameter estimates
Also in this second model it is difficult to determine precise values for all the parameters.To have the possibility to compare the two model, we assume the same previous parameters when they coincide.The others are estimated from both modeling and biological literature.Here, Q M and Q C are estimated following the indication of [10].

Comparison between PM and no-void simulations
Here, we want to compare simulations of no-void and PM model.To this goal, we have to choose the same values of parameters whereas it is possible.A first difference between them is in the killing parameter λ, this is due to the structure of the two models: in the no-void model there are only macrophages and bacteria, then, in some way, the term BM is greater than in PM model.For this reason we choose a coefficient λ P M > λ no−void so to obtain a comparable behavior between the two models.The no-void model considers only one-half granuloma, so we have to double the number of bacteria and macrophages in PM model to compare them.
As initial condition, we choose 50 bacteria, 0 macrophages and an initial radius R(0) = 50 µm.In PM model, we assume to have an initial infection of 100 bacteria spread over a length of 100 [µm], and about 780 macrophages homogeneously spread everywhere except for the region filled by bacteria.For simplicity we assume C(x, 0) = 0.
We present simulations over a period of time of two weeks.Regarding the no-void model, in figure 11 we can see the evolution of number of macrophages (on the right) and bacteria (left) inside granuloma, while, in figure 12, is shown the behavior of the granuloma's radius during two weeks.In figure 13 are reported simulations done by PM model: the evolution of bacteria (on the left) and macrophages (on the right) concentration.Then, in figure 14 we have the evolution of the number of bacteria.Now, we can conclude that, under our hypotheses, the two models are not very different between them.In the no-void model, there is a greater initial growth of bacteria than in pressure model, but, after a while, they are contained by macrophages.This different initial growth can be due to the fact that in pressure model macrophages are present initially over all the region, then they are in touch with bacteria at t = 0.In the no-void one, macrophages arrive from outside, so they are not in touch with bacteria initially and the contact surface with bacteria is only on one side, at least in the beginning.The number of macrophages surrounding bacteria are comparable in both models, and the radius of no-void model reaches a maximum of about 0.5 [mm].In the pressure model, we have that the spread of bacteria infection has a radius of about 2 [mm].
These simulations permit us to consider the two model as comparable, even if with some difference.It is difficult to establish which is the better in absence of experimental data.We think that the PM model gives a better justification of reaction-diffusion equations used, because we have smaller value of concentrations.It is difficult think as true the diffusion of a cell in an area completely full of other cells.

Simulations in the two-dimensional case
In this section we use the pressure model in two-dimensions.In the first part of the section we simulate a case in which bacteria prevail, then, a case in which bacteria are contained.Equations of the model are given by 19-25.

Numerical simulations
We guess to have initially about 10 5 macrophages randomly distributed in an area of 1 [dm 2 ], and in the middle of this area we assume to have an initial infection of 400 bacteria.Under these assumptions we have simulated the evolution of infection in the first two months, using the parameters in table below 3

Bacteria prevailing
In the first case we use the non-dimensional value λ = 4.5 • 10 5 .It gives origin to an uncontrolled growth of bacteria.Here, we present the behavior of bacterium and macrophage concentration in the first two months.In figure 15 we can see their initial distribution, in figure 16 we can see their distribution after 1 month and in figure 17 after 2 months.In figure 17 we can observe the clustering of macrophages in the middle of region, due to the great infection located there, which attracts macrophages by chemotaxis.Finally, we show the number of bacteria evolution in the first two months in figure 18.

Bacteria contained
Here, we present simulations describing a case in which bacteria are contained by macrophages.We use the non-dimensional value λ = 9.5 • 10 5 .In figure 19 we can see their initial distribution, in figure 20 we can see their distribution after 1 month and in figure 21 after 2 months.Finally, in figure 22 we show the number of bacteria in the first two months.

Comparison with different efficiency coefficients
In figure 23 we compare the behavior of the infection (number of bacteria) varying the coefficient λ from 5 • 10 5 , case in which bacteria increases their number exponentially, and λ = 9.5 • 10 5 , case in which infection is controlled by innate immune response.It is evident that for λ ≤ 6.3 • 10 5 we have an uncontrolled growth.

Variation of immune response
In this section we want to present a couple of examples in which the immune response can change with time.As first example, we want to simulate a sort of reactivation of infection.It can be due to two different causes: the first one is a diminishing number of macrophages; the second one is a diminishing efficiency to kill bacteria.These two causes are simulated in section 8.2.1.
As second example, we want to represent a strengthening of immune system, a sort of simplified adaptive immune response.To make so, we assume that after three weeks, the efficiency of immune system (λ coefficient) increases its value of three times, and the macrophages death rate (µ) diminishes.

Case I: Reactivation
a ) Here, we assume to have a constant non-dimensional macrophage death rate µ = 10.We suppose, also, that there is no incoming macrophages from outside, then the number of macrophages is diminishing, weakening the immune system.Choosing λ = 9.5 • 10 5 we have that initially macrophages contain infection.But, after about six months, when the number of macrophages is 50% the infection of bacteria arises again.This way, we want to represent a sort of simplified reactivation, in a brief time with respect to reality (years).In figure 24 we can see the concentration of bacteria at several time points.Then, in figure 25 is shown the variation of the number of bacteria (upper) and the number of macrophages (lower).b ) Here, we guess that the efficiency of immune system is lower: λ = λ(t) is a function of time, it varies linearly from 9 • 10 5 ( at t = 0) to λ = 4 • 10 5 (for t = 8 months).The death rate lower than the previous case, it is µ = 3 • 10 5 .Under these assumptions we obtain the behavior of bacteria concentration in figure 26 and the evolution of the number of bacteria in figure 27.In this case the number of macrophages diminish about 15% in eight months.

Case II: Strengthening of immune response
Here, we guess after three weeks the non-dimensional coefficient λ increases, and the death rate of macrophages µ (non-dimensional) becomes smaller.This means that after three weeks the immune response improve (adaptive immune response) and the macrophage death rate diminishes for their improved effectiveness.Using µ = 3 in the first three weeks, after that µ = 0.5.While λ = λ 1 in the first three weeks, then λ increases linearly to the value of λ 2 in the next 9 weeks, then it is λ = λ 2 up to eighth month.We choose four cases: a) λ 1 = 3 • 10

Granuloma formation
Finally, we make a long time period simulation, assuming the same conditions of section 8.2.2, except for µ, which is, here, µ = 3 in the first three weeks, and µ = 0 in the following time.This because we suppose that, for a long period of time, the immune system can maintain its efficiency.We can see one-year simulation in figure 29, where the number of bacteria, after an initial growth, is contained, but not canceled by macrophages, its initial value is 400, and its asymptotic value is about 100.

Conclusions
In this work, we present a qualitative approach to the immune response to Mtb aggression, which is inspired by prey-predator systems.We propose a model, which works in every space dimension, by using mass balances equations and a pressure law for the velocity.We find, at least in 1-dimension, a region of linear instability, which is limited by the non-linear terms.We also observe, in long time simulations, the persistence of a small group of bacteria in equilibrium with macrophages.The reactivation mechanism has been also discussed.
Next step will be the development of new models to take into account phagocytosis, activation of cells, T-cells and all the other components which play an important role in the Mtb evolution.

Fig. 5 Fig. 6
Fig. 5 Number of bacteria with M ≈ 10 α λ on the left, and M ≈ 0.1 α λ on the right.

Fig. 7
Fig. 7 Initial perturbation of bacteria is on the left.On the right side, we have Ba (in red) and B b (in blue) after 10 years simulations.

1 .
Case I: 170 macrophages.The evolution of the number of bacteria with different λ values is in figure 8. There, macrophages are able to contain infection if λ > λ c = 3 • 10 4 .2. Case II: 430 macrophages The evolution of the number of bacteria with different λ values is in figure 9.In this second case, macrophages are able to contain infection if λ > λ c = 1.1 • 10 4 .3. Case III: 690 macrophages The evolution of the number of bacteria with different λ values is in figure 10.In this third case, macrophages are able to contain infection if λ > λ c = 6.6 • 10 3 .

Fig. 9 Fig. 10
Fig. 9 Number of bacteria with different values of λ and with 430 macrophages initially present.

Fig. 18
Fig. 18 Number of bacteria growth in 2 months.Case in which bacteria are prevailing with λ = 4.5 • 10 5 .

Fig. 23
Fig. 23 Comparison among number of bacteria with different values of λ coefficient.

Fig. 25
Fig. 25 Number of bacteria (upper) and number of macrophages (lower) in the first 8 months.

Fig. 27
Fig. 27 Number of bacteria in the first 8 months.

Fig. 29
Fig. 29 Number of bacteria evolution in 1 year.
chemoattractant C(x, t) is produced by extracellular bacteria B, at a maximal rate σ B .It comes from a Michaelis-Menten type equation, where σ B represents the saturation.Moreover, chemoattractant diffuses at a rate D C .It has a natural decay rate, Γ C C, and is used by macrophages M , at a rate Γ M .The governing equation for chemoattractant C(x, t) is

Table 2
Non-dimensional parameters table Fig. 8 Case I. Number of bacteria with different values of λ and with 170 macrophages initially present.