An Agent Based Model of Tularemia

Francisella tularensis is a formidable intracellular pathogen. Upon inhalation it leads to systemic disease (tularemia) with a high mortality rate. We sought here to develop a computational tool for infections with this class. A bio-threat microbe for which extensive datasets are not available to infer parameters that might determine the outcome of the infection. We present a two-compartment agent based model that simulates inhalational tularemia with subsequent dissemination to the liver and incorporates experimental data and validated general parameters of host defense mechanisms. This systems approach suggests that the initial number of macrophages, the probability of dissemination, and the initial clearance rate of bacteria correlate with the outcome of infections with Francisella . These findings underline the importance of early innate immune defense mechanisms in the prevention of tularemia.


Introduction
A wide variety of individual host cell defense mechanisms and bacterial virulence factors have been described and their singular importance during infection has been experimentally established in many cases. Computer models of infectious diseases can provide a systems biology approach for the inference of critical parameters of host-pathogen interaction. However, few such models exist and they are particularly scarce for infectious processes of the lung. Disease models for infections with Mycobacterium tuberculosis were able to identify chemokine secretion and the efficiency with which T cells activated macrophages as important parameters for the outcome of latent versus active infection [1]. Modeling of influenza infection demonstrated the importance of the spatial structure of the initially infected cell [2]. In another model the growth rate of the parasite was found to be the most important parameter for infection with Leishmania [3]. While model-building is frequently done to integrate large and complex existing datasets, we sought here to apply computer modeling to the infection with the class. A bio-threat agent Francisella tularensis, which results in the systemic disease tularemia and for which limited experimental data are available. Francisella tularensis is a category of bio-threat because of its high infectivity, its high case-fatality rate when untreated, and because of previous attempts to 'weaponize' this microbe [4]. The ability of Francisella tularensis to cause disease is closely correlated with its ability to enter into and survive inside of what is currently thought to be its permissive host cell, the macrophage [5][6][7]. Inhalational tularemia progresses from necrotic pulmonary foci to systemic disease with the formation of ill-defined granulomas. During the initial phase of infection with Francisella via inhalation bacteria accumulate in macrophages or leukocytes in pulmonary alveolar ducts. From the second day on multiple foci of necrosis can be detected with significant macrophage infiltration. After the second day of infection lesions can also be detected in liver and spleen [8]. Francisella can be found in hepatocytes, where they rapidly recruit macrophages, which lyse hepatocytes with subsequent release into the extracellular space [9].
The goal was to investigate if a systems biology approach for the case of a rare disease such as tularemia could infer critical parameters for the outcome of this infection, which could then become the focus of further targeted investigations.

Mathematical model
In our model, the immune system is modeled by cells which are treated as C++ objects, and move according to a biased random walk in the direction of the chemokine gradient. This random walk is the solution of the chemotaxis PDEs [10].
where X C (t) is the position of the cell, and φ (X,t) is the concentration of chemoattractant, and ξ is a white noise variable, α, λ, β and D1 are constants.
The implementation of our model is based on Segovia-Juarez et al. [1]. Chemokine diffusion is implemented with the following rule: in each micro-compartment (i,j) the chemokine concentration C i,j diffuses to and from the four micro-compartments in its immediate (von Neumann) neighborhood: where λ is a diffusion constant. λ measures the proportion of C i,j that diffuses out micro-compartment (i,j) during each time-step. We calculate a value for λ from λ = 4 c ∆t/∆x where ∂ c is the diffusion constant for chemokine molecules in the diffusion PDE, x  is the scale of spatial discretization, and ∆t is the scale of time discretization. Using ∆t = 0.1 min and ∆x = 10 −5 m, respectively. Using a value of ∂ c = 10 -7 we obtain a value of λ=0.6.
time-step, a certain proportion δ decays as in: Extracellular bacteria B E replicate with an estimated doubling time of 3-5 h. We assume an upper bound of B E , K BE =400 within each microcompartment, which is twice the number for M. tuberculosis, since the Francisella tularensis LVS bacterium is half the size of M. tuberculosis. Extracellular bacterial replication follows logistic growth with respect to this upper bound: In addition we have each bacterium secrete 1 unit of chemokine per unit time step.
Macrophages are represented as discrete agents that reside on the lattice. At most one macrophage can occupy a given microcompartment. Macrophages have the following attributes: position, age, the number of intracellular bacteria (B I ), and the state of the macrophage. B I is a non-negative integer such that 0 ≤ B I < K BI where K BI is a parameter representing the average intracellular bacterial carrying capacity of macrophages. The state is defined from one of the following: resting (M R ), infected (M I ), chronically infected (M C ), or activated (M A ). Macrophages move according to a random walk biased towards neighboring micro compartments with higher concentrations of chemokine.
The rules for resting macrophages are as follows: We assume that if there are a small number N RK = 2 of bacteria in the same micro-compartment as a resting macrophage, then the macrophage phagocytoses and kills those bacteria and remains in the resting state. If there are more than N RK bacteria present, then there is some small probability P K that the macrophage still succeeds in killing N RK intracellular bacteria. If it cannot it becomes infected (MI) with an intracellular level of N RK .
The rules for infected macrophages are as follows. At every timestep, each infected macrophage releases a given amount C I = 5000 units of chemokine into its micro-compartment. The model is invariant to changes in C I , for 1 Intracellular bacteria replicate within infected macrophages. The discrete intracellular growth rate α BI 2 is estimated to be between 0.002 and 0.006/min: If the intracellular bacterial load exceeds a threshold N C = 10, the macrophage becomes chronically infected.
The rules for chronically infected macrophage rules are as follows. A chronically infected macrophage secretes chemokine into its micro-compartment. We assume that within chronically infected macrophages, intracellular bacteria replicate logistically with respect to the bacterial carrying capacity of macrophages, K BI .
The growth rate α BI 2 is estimated to be between 0.006/min and 0.008/min. We use a value of K BI = 50. If the intracellular bacterial load of a chronically infected macrophage exceeds the bacterial carrying capacity of macrophages (K BI ), the macrophage bursts and its intracellular bacteria are released.
We model necrosis of lung tissue by counting how many times a chronically infected macrophage bursts. A micro-compartment (i,j) is declared to be necrotic if the number for bursting exceeds N necr , we use a value of N necr = 8.
We model dissemination of macrophages from the lung to the liver as follows. If macrophages move to a sink compartment they will move to the liver with a probability M diss between 1% and 20%. Because macrophages are disseminating to the liver, we have them replenished for the liver and lung according to the following rule. Macrophages are recruited with a probability of ( ) , max 1,0.001 i j recr C M * + , where C i,j is the concentration of chemokine at the source location with coordinates i and j at which the macrophage is recruited. We estimate M recr for both the lung and liver to be between 0.2 to 0.7.
Because the macrophages become infected in 1 minute (10 time steps) after the inoculation, and the movement of macrophages is updated every 100 time steps, we assume that there are initially three infected macrophages which are placed 1 μm away from a sink compartment, as infected macrophages can enter the bloodstream within 1 minute.

Infection assays
Francisella tularensis subspecies holarctica vaccine strain (F. tularensis LVS, army lot 11) was generously provided to us by Dr. Karen Elkins (FDA). Francisella was grown on chocolate II agar enriched with IsoVitaleX (BD Biosciences, San Jose, CA) for 40-48 hrs at 37°C. As liquid medium we used Mueller-Hinton broth supplemented with IsoVitaleX RAW 264.7 murine macrophages and A549 type II alveolar epithelial cells were obtained from ATCC. Dulbecco's Modification of Eagle's Medium (DMEM; Cellgro) was supplemented with 10% fetal bovine serum (Hyclone, not heat-inactivated) and penicillin (100 IU/ ml) and streptomycin (100 μg/ml). When cells were used for Francisella infection assay no antibiotics were added 24 h prior to infection. Cells were grown at 37°C and 5% CO 2 . Infection assays were performed as described in Tamilselvam and Daefler [11].

Agent-based model for infection with Francisella
We base our model on a tested agent-based model of mycobacterial infection of the lung [1,12]. This model, which approximates the solution of a system of Partial Differential Equations (PDE), comprises the basic innate immune defense mechanisms during infection and goes beyond the purely deterministic attempt of modeling with ordinary differential equations (ODE) [1]. We modified this agent-based model as described below and implemented it with experimental data from our laboratory and from the literature for infection with Francisella. Our model uses parameters of the human immune system and assumes infection with Francisella tularensis subsp. tularensis. When appropriate, data from infection of mice with Francisella tularensis subspecies holarctica vaccine strain (F. tularensis LVS) were incorporated.
Since infection with Francisella represents a systemic disease with dissemination to the liver, we designed a novel two-compartment agent-based model of Francisella infecting lung and liver. Recruitment of macrophages occurs at "source" microcompartments located at microcompartments 400 μm away from the four corners of the lattice. Macrophages enter the lattice with a probability M necr at these locations. These source compartments represent locations where blood vessels enter the lung tissue through which new macrophages arrive at the infection site. In order to model dissemination into the bloodstream, and then into the liver, we designate certain compartments in the lung,  μm away from the center of the lattice as "sink" compartments. These sink compartments represent locations where blood vessels enter the lung tissue through which macrophages leave the infection site. Macrophages leave the lung with probability M diss . We also introduced replenishment of macrophages in the lung and liver in response to chemokines and 'chemoattractants' elaborated by bacteria.
An agent-based model is a generalization of a "lattice-gas" cellular automation, in which particles move about and interact in a prescribed fashion [13]. In an agent-based model, elements of the system are represented as discrete agents, which interact with each other and with the environment according to sets of rules. In our model, the immune system is modeled by cells that are treated as C++ objects, and move according to a biased random walk in the direction of the chemokine gradient. This random walk is the solution of the chemotaxis PDEs (equations 1 and 2 in Materials and Methods) [10,14]. Macrophages and T-cells are modeled individually as cells in a square region of organ tissue with source and sink compartment where they enter and leave the lung. Extracellular and intracellular bacteria are modeled as continuous variables. Rules describe events such as the recruitment, infection and bursting of macrophages, the killing of bacteria by macrophages, and activation of macrophages by T-cells. Our model is a two-component model, with one component being the lung, and the second component being the liver. Marino and Kirschner [15] have built a two component model of M. tuberculosis infection of the lung and lymph nodes. We model this by making two identical copies of the same model with different parameters, corresponding to the parameters of the lung and liver, denoting the lung with the index 0 and the liver with the index 1. The extracellular growth rate of the bacteria in the liver and the recruitment probability of macrophages in the liver are different from those in the lung, the rest remain the same. We also model the same square area of liver and lung tissue. Dissemination is modeled by changing the component from 0 to 1 with a probability of M diss , when a macrophage arrives at a ``sink" compartment. Dissemination to the liver takes one time-step or 6 seconds. We assume there are no macrophages in the liver initially.
We model 2 mm square of organ tissue (liver or lung) as a 100×100 lattice of micro-compartments. Each micro-compartment can simultaneously hold at most one macrophage, one T-cell, chemokine, and bacteria. The entities of the model consist of objects, which are macrophages and T-cells and continuous variables representing chemokine and bacteria. The chemotaxis of macrophages and T-cells is modeled using a biased random walk for each cell, biased by the chemokine gradient. Infected macrophages and bacteria secrete chemokines, which is a composite of all the chemo-attractants from the bacteria and cytokines from the macrophages. Chemokines diffuse and degrade according to equation 3.
Time is discrete in this model; each time-step corresponds to 6 seconds. Macrophages are modeled as objects with the attributes of position, age, the number of intracellular bacteria and the state of the macrophage. The state of the macrophage is either resting (MR), infected (MI), chronically infected (MC), or activated (MA). Infected macrophages release chemokines at each time-step. If the intracellular bacterial load exceeds a certain threshold, the macrophage becomes chronically infected. Chronically infected macrophages continue to secrete chemokines. If the intracellular load of a chronically infected macrophage exceeds the bacterial carrying capacity, the macrophage bursts and its intracellular bacteria are released. Necrosis of lung tissue results from this bursting of macrophages. A micro-compartment (i, j) is declared to be necrotic if the number of burstings exceeds N necr =8 as in Segovia-Juarez et al. [1]. Recruitment of macrophages occurs at "source" microcompartments located at microcompartments 400 μm away from the four corners of the lattice. Macrophages enter the lattice with a probability M necr at these locations. These source compartments represent locations where blood vessels enter the lung tissue through which new macrophages arrive at the infection site. In order to model dissemination into the bloodstream, and then into the liver, we designate certain compartments in the lung, 200 μm away from the center of the lattice as "sink" compartments. These sink compartments represent locations where blood vessels enter the lung tissue through which macrophages leave the infection site. Intracellular bacteria also follow logistic growth according to the equation 7. Here, intracellular growth, B I , is the continuous variable for intracellular bacteria, and K BI is the upper bound on the number of intracellular bacteria, which is also determined by the size of the bacterium. For Mycobacterium tuberculosis, K BI is 20, for Francisella tularensis K BI is 50.
Some parameters of the immune system are fixed. The speed of resting and infected macrophages has been found by Webb et al. [16] to be 1 μm/min for resting macrophages and 0.0007 μm/min for infected macrophages. The lifespan of macrophages are fixed at 100 days for a resting macrophage [17]. We also assume that the number of bacteria killed by a resting macrophage is 2. Other parameters of the immune system need to be estimated. The chemokine diffusion coefficient and the chemokine degradation coefficient were estimated between 0.5 and 0.8 per time step and 0.000288 and 0.0011 per time-step respectively [18]. The initial number of macrophages was estimated from Mercer et al. [19] and Stone et al. [20] to be between 40 and 400. The dissemination probability of a macrophage from the lung to another organ is estimated to be between 1% and 20%. The probability of a resting macrophage killing 2 bacteria is estimated between 1% and 10%.

Determination of Francisella's growth parameters
Growth rates for Francisella tularensis subspecies holarctica vaccine strain (F. tularensis LVS, army lot 11) were determined in a model infection of murine macrophage cell lines (Raw 264.7) and human type II alveolar epithelial cells (A549). We determined the intracellular growth rates for these models as between 0.002/min and 0.006/min for infected, and between 0.006/min and 0.008/min for chronically infected cells. While the initial rate of infection varied, no differences in growth rates were found between Raw 264.7 and A549 cell lines.
These findings are in agreement with results described in the literature, where Francisella follows logistic growth with doubling time of 3-5 hours [9,21] for the lung and liver, which gives calculated growth constants of 0.00231/min and 0.00785/min respectively.

Parameters that determine the outcome of infection
The model as described above was used for simulation of infection with Francisella for a period of 5 days. Nine parameters were allowed to vary in order to determine how they would correlate with outcome of infection (measured here as an increase or decrease in the number of extracellular bacteria): the diffusion coefficient λ, the degradation coefficient δ, the intracellular growth rate for infected macrophages α BI1 , the intracellular growth rate for infected macrophages α BI2 , the initial number of macrophages, M init , the recruitment probability of macrophages to the lung, M recr , the recruitment probability of macrophages to the liver, M recrL , the dissemination probability M diss , and the probability of a resting macrophage killing bacteria P K .
We performed Latin Hypercube Sampling [1,22,23] on these nine parameters of the model. The Partial Rank Correlation Coefficient (PRCC) with the number of extracellular bacteria in the liver and lung was calculated as shown in table 1. From this we see the parameter M diss , the initial number of macrophages, is most strongly correlated with high numbers of bacteria in lung and liver. M diss , the dissemination probability is negatively correlated with the number of bacteria in the liver, but this seems largely due to the stochastic nature of the model, where a dissemination probability higher than 10% can result in no or little liver dissemination. With a higher dissemination probability, macrophages migrate to the liver rather than stay in the lung and hence do not become infected.
The intracellular growth rate for infected macrophages also has a significant correlation with the number of extracellular bacteria in the liver and the lung, but the intracellular growth rate for chronically infected macrophages is not as significantly correlated. This suggests that slowing the growth of Francisella in the early phases of infection could have a significant impact on the outcome of infection.

Recruitment of macrophages correlates with dissemination of infection
We further investigated the effect of the dissemination probability on the number of bacteria in liver and lung when the initial number of macrophages was kept constant. Simulations were run and each individual outcome was plotted as the number of bacteria in lung ( Figure  1A) or liver ( Figure 1B) after 5 days versus dissemination probability. The number of bacteria in the lung decreases monotonically with increasing dissemination probability. This levels off at a probability of 10%.
The number of extracellular bacteria in the liver is within a range between dissemination probability 1% and 10%, and levels off roughly at dissemination probability of 5%. Above dissemination probability 10%, there is some noise, where there is little or no dissemination to the liver ( Figure 1B). The reason for this is that at higher rates of dissemination, more macrophages leave the lung and cannot become infected. To reduce this noise, we increased the chemoattractant secreted by the bacteria to be 5000 units per bacterium. The following plot shows that increased chemoattractant leads to more liver dissemination at M diss above 10% ( Figure 2C).
These findings suggest that the likelihood with which macrophages are recruited to the lung and leave this compartment are important factors in the pathogenesis of systemic disease. Very little quantitative data are available about the true dynamics of macrophage recruitment  and release in the lung, which are very difficult, if not impossible, to measure. Our simulations suggest, however, that a balance of recruitment and release plays an important role and should be further experimentally investigated.

Initial number of macrophages is significant
Since the outcome of infection in terms of bacteria in lung and liver appears to level off at a dissemination probability of 5%, we investigated the influence of the initial number of macrophages on bacterial proliferation. At a constant dissemination probability of 5% the simulation was run 10 times for each set of parameters with the initial number of macrophages between 40 and 400.
Bacterial growth in the lung (Figure 2A) or in the liver ( Figure 2B) was plotted against the initial number of bacteria in the lung for each event. The number of extracellular bacteria in liver and lung increases monotonically with the initial number of macrophages. This number, due to the probabilistic nature of the model falls within a range of values, with the error bars of size between 10109.32 and 46139 in either direction. There is no clearance of bacteria either in the lung and the liver.
We further simulated the time-course of infection over a five days period for 40, 100, and 400 initial macrophages with a dissemination probability of 5% (Figure 3). One can observe that in spite of the stochastic nature of the model, there is not a large amount of variability. With the initial number of 40 macrophages, there are more bacteria in the lung than in the liver ( Figure 3A). As the initial number of macrophages increases to 100, the number of extracellular bacteria in the liver eventually increases above the number of extracellular bacteria in the lung ( Figure 3B). With 400 initial macrophages we can see that there are more bacteria in the liver than in the lung ( Figure 3C). There is also some stochastic variability in this simulation, but this variability is consistent with the fact that the number of extracellular bacteria in the liver and lung fall within a range of values, as described above.

Probability of clearing infection
With the initial number of macrophages and their dissemination dynamics being the most important factors to determine the proliferation of Francisella, we reasoned that the probability with which a macrophage might clear its intracellular infection might be pivotal for the overall outcome. We introduced a probability of activation for resting macrophages (P act ) to see if the infection clears. P act was varied between 1% and 100%. For each time the movement of a resting macrophage was updated, the macrophage could be activated with this probability. We then simulated the infection with 100 initial macrophages. We found that a probability greater than 4% leads to clearance within 10 days. We then set the probability of a resting macrophage killing 2 bacteria at 10% and placed 5 bacteria in the center of the grid instead of infected macrophages. The macrophages then could either kill the bacteria or become infected. We ran 10 runs for each value of the initial number of macrophages. We found that between 250 and 400 initial macrophages, there were runs in which all the bacteria were phagocytosed, and that the infection cleared immediately. The probability of clearance this way was 30% clearance with 250 initial macrophages, 10% with 300 initial macrophages, 30% with 350 initial macrophages and 30% with 400 initial macrophages. A representative simulation is demonstrated in figure 4.

Comparison with M. tuberculosis infection
The strong positive correlation of M init with the number of extracellular bacteria in the liver and in the lung contrasts with infections with M. tuberculosis where M init is not significantly correlated with this number. To further investigate this difference, we ran our simulation with the parameters for Mycobacterium tuberculosis with a dissemination probability of 5% and found that there was no liver dissemination, and that the infection eventually cleared out of the lung ( Figure 5A). To further test our model, we eliminated the replenishment of macrophages. In this case, there is also no dissemination to the liver, as the infection cannot spread to sink microcompartments within 10 days ( Figure 5B).
T-cells then arrive and activate most of the infected macrophages, which leads to a "wall" forms around the infection, where all of the macrophages surrounding the infection are either resting or activated. As the graph shows, the infection clears eventually and no infected macrophages can enter "sink" compartments to disseminate to the liver. Only activated macrophages able to disseminate to the liver, since the macrophages surrounding the extracellular bacteria all become activated before they can get to the "sink" compartments. In the case of F. tularensis, all of the macrophages in the granuloma surrounding the extracellular bacteria would be infected.

Discussion
When Francisella is inhaled it rapidly leads to systemic disease with a high mortality rate. We present here an agent-based model that simulates inhalational tularemia with subsequent dissemination to the liver as a two-compartment system. We used statistical analysis from multiple simulations to infer factors that influence the outcome of the infection. The current model suggests that the initial number of macrophages, the probability of dissemination, and the initial clearance rate of bacteria by innate immune defense mechanisms correlate with the outcome of infection with Francisella.
Agent-based modeling of infectious diseases has been employed previously for tuberculosis and leishmaniasis [1,3]. The immunological parameters and host defense mechanisms that are the basis for our model are the same as those which have been used in these validated agent-based models. However, when we introduce Francisella-specific attributes we identify different parameters that are crucial for the outcome of the disease process. The model presented here points to a preeminent role for early defense mechanisms (number of macrophages recruited and rate of initial clearance) to thwart a successful infection. This contrasts to infections with M. tuberculosis, where mounting a delayed specific immune response is the more important host defense and where variations in the number of macrophages recruited do not result in a clearance of the infection. A novel aspect of our model is the two-compartment system, which demonstrates the importance of dissemination as a measure to evade clearance mechanisms.
The power of this modeling approach is illustrated that by integrating pathogen-specific parameters such as growth rate, intracellular replication rates, and immuno-stimulatory effects as measured by the release of chemoattractants, disease-specific outcomes are modeled that correlate with experimental evidence. This also demonstrates how relatively simple parameters are key determinants in host-pathogen interactions and can be accurately assessed in an agent-based model. Even though parameters such as growth rate might be simple, they are the outcome of complex interactions of host and pathogen metabolism, elaboration of virulence factors, and host defense pathways. However, the power of our model is to allow variations of such summary and individual parameters, some of which are difficult or impossible to determine, and then correlate them with the outcome of infection in order to infer the dominating parameters for pathogenesis. This is of particular importance for infections with a less well studied microbe such as Francisella. Our modeling approach might thus delineate suitable parameters in lieu of extensive experimentation that can then be further investigated.