A Novel Stochastic Multi-Scale Model of Francisella tularensis Infection to Predict Risk of Infection in a Laboratory

We present a multi-scale model of the within-phagocyte, within-host and population-level infection dynamics of Francisella tularensis, which extends the mechanistic one proposed by Wood et al. (2014). Our multi-scale model incorporates key aspects of the interaction between host phagocytes and extracellular bacteria, accounts for inter-phagocyte variability in the number of bacteria released upon phagocyte rupture, and allows one to compute the probability of response, and mean time until response, of an infected individual as a function of the initial infection dose. A Bayesian approach is applied to parameterize both the within-phagocyte and within-host models using infection data. Finally, we show how dose response probabilities at the individual level can be used to estimate the airborne propagation of Francisella tularensis in indoor settings (such as a microbiology laboratory) at the population level, by means of a deterministic zonal ventilation model.


INTRODUCTION
Francisella tularensis is a gram-negative, facultative bacteria and the causative agent of tularemia (Oyston et al., 2004;Oyston, 2008). Due to its high infectivity and ability to cause a debilitating disease with the inhalation of as few as 10 organisms, F. tularensis has been classified as a category A bioterrorism agent by the Centers for Disease Control and Prevention (CDC). Following inhalation, bacteria are deposited in the lungs where, to begin with, they are primarily taken up by alveolar phagocytes through phagocytosis, as described by Hall et al. (2008). By escaping the Francisellacontaining phagosome (FCP), bacteria enter into the cytosol of the phagocyte. F. tularensis can resist killing in the cytosol from reactive oxygen species (ROS) and can subsequently undergo multiple rounds of division within the host cell. Following this intracellular bacterial replication, the host phagocyte ruptures and dies, releasing its bacterial load back into the extracellular environment (Cowley and Elkins, 2011). For up to 72 h post-infection, F. tularensis is capable of preventing immune recognition. Therefore, it is important to understand how an individual may react to the infection, and when they develop tularemia.
Dose response models have been developed in an attempt to quantify the risk to a population associated with chemical and biological agents. However, unlike with chemical agents where the initial dose is the total amount available to cause a response, the ability of biological agents to reproduce post-infection means that the extent of replication within the host must be taken into account (Huang and Haas, 2009). Furthermore, since this timescale of infection is in the order of days, and since the window of opportunity for effective medical treatment is often limited, a better understanding of the infection timescale could provide valuable information to guide optimal treatment strategies. Attempts have therefore been made to incorporate time into such dose response models. Many of these original approaches involved adjusting existing dose response models, such as the classical exponential and beta-Poisson models, or probit analyses to allow for time dependency of the model parameters (Chen, 2007;Huang and Haas, 2009). However, by choosing convenient statistical distributions, the link between the dose response model and the underlying within-host biological mechanisms that govern the level of bacterial replication is tenuous. A stochastic mechanistic model is proposed by Pujol et al. (2009) for the within-host interaction dynamics between immune effector cells and pathogens, which takes into account both the total dose inhaled by the host and the total exposure period during which this dose is inhaled. It is also worth mentioning the work by Gillard et al. (2014), where a stochastic within-host computational model is proposed for the infection process, in the BALB/c mouse, following inhalational exposure to Francisella tularensis SCHU S4. By focusing on a compartmental agent based model, Gillard et al. (2014) consider the intracellular dynamics of a single infected phagocyte, and model the stages of bacterial replication and phagocyte rupture as a birth process with catastrophe, where the number of bacteria released in a single rupture event follows a geometric distribution. The average number of bacteria released is then estimated using the mean of this geometric distribution.
Another recent example is the Markov chain model described by Wood et al. (2014), which addresses these issues by considering the key interactions between F. tularensis bacteria and host (human) phagocytes within the lung space. Using the Markovian nature of the process, the probability and time for the total number of bacteria to reach some threshold can be computed, this threshold being identified as the necessary amount of bacteria needed for host illness onset. Despite this, fitting procedures are still used to obtain quantities, such as the time until a single infected phagocyte ruptures, which are required to parameterize the model. A particular limitation suggested by Wood et al. (2014) is the consideration of a deterministic (constant) time for the time to rupture of each infected phagocyte. This does not allow for modeling the experimentally observed variability in this time among different phagocytes, where in fact a log-normally distributed rupture time is predicted by Wood et al. (2014), but not explicitly incorporated into the model. Also, by using a deterministic approach to modeling the intracellular growth of F. tularensis bacteria, Wood et al. (2014) assume a constant number of bacteria released on rupture of any infected phagocyte, not accounting for the existing variability in the number of bacteria released by different phagocytes.
In this paper, an extension to the model described by Wood et al. (2014) is proposed. By incorporating a second, within-phagocyte, model into the existing within-host model, the stochastic intracellular dynamics of F. tularensis can be replicated. This can account for the log-normally distributed rupture time, leading to a rupture size probability distribution (i.e., number of bacteria released upon phagocyte rupture) which enables us to account for inter-phagocyte variability at the withinhost level. Thus, the within-phagocyte model can be linked with the within-host model for the interaction between extracellular bacteria and susceptible phagocytes by means of the distribution of the number of bacteria released by a single infected phagocyte, obtained from analyzing the within-phagocyte model, which allows for varying phagocyte rupture sizes in the within-host model. In summary, this multi-scale model allows us to relax the assumption made by Wood et al. (2014) that a fixed number of bacteria is released from every single infected phagocyte on rupture. For both the within-host and within-phagocyte models, analytical approaches to calculate the summary statistics (dose response probability and mean time until response) defined by Wood et al. (2014) are outlined. However, by exploiting the structure of the resulting Markov processes, more efficient approaches than the methods proposed by Wood et al. (2014) are described here. Finally, a zonal ventilation model for the indoor airborne spread of F. tularensis is presented in order to illustrate how dose response probabilities at the individual level, computed from the within-host model, can be used in order to make predictions at the population level.

MATERIALS AND METHODS
In this section, our aim is to develop a multi-scale model for the infection dynamics of F. tularensis bacterium, by linking a within-phagocyte, a within-host and a population-level model. In section 2.1 we develop a stochastic within-phagocyte model for the infection dynamics of a single phagocyte by F. tularensis. We show how the log-normally distributed rupture time estimated by Wood et al. (2014) from experimental data (Lindemann et al., 2011), can be incorporated into this model, while maintaining the Markovian nature of the underlying stochastic process, and how first-step arguments allow one to compute the probability distribution of the total number of bacteria released by an infected phagocyte upon rupture. This distribution is used in section 2.2 to link the within-phagocyte model to the withinhost model for the interaction between extracellular bacteria and phagocytes within the host. This within-host model accounts for inter-phagocyte variability in the amount of bacteria released upon rupture. The aim of the within-host model is to compute the probability of host response (in terms of the onset of symptoms), as well as the time to this response. Finally, we illustrate in section 2.3 how these dose response probabilities at the individual level might be used for predicting, at the population level, the number of individuals showing symptoms upon indoor release and airborne spread of F. tularensis, by means of a zonal ventilation model and under different ventilation settings in an hypothetical microbiology laboratory.

Within-Phagocyte Model
The first level of the multi-scale model is a within-phagocyte model, used to replicate the dynamics of an F. tularensis bacterium after it has been ingested by a host phagocyte, assuming that the bacterium escapes the FCP, entering into the cytosol and starting replication. Phagocytosis leading to successful bacterial killing will be considered in the within-host model, and is not analyzed here. This includes the replication of bacteria within the cytosol, and the subsequent rupturing and death of the phagocyte (Cowley and Elkins, 2011). These stages of the intracellular lifecycle can be modeled using a continuous-time stochastic process X = {X(t) : t ≥ 0} that follows the structure of a birthand-death process with catastrophe (Karlin and Tavaré, 1982;Di Crescenzo et al., 2008), where X(t) is the number of bacteria within the phagocyte at time t ≥ 0. In particular, replication and death of bacteria within the phagocyte can be modeled as a stochastic logistic growth process over states in N = {1, 2, . . . }, representing the number of bacteria contained within the cytosol (see Figure 1A). Birth and death rates for state n ∈ N are obtained by following arguments by Allen (2003, section 6.8), where we assume that each bacterium replicates independently of all others at rate λ, so that: (1) We denote by λ [hours −1 ] the per bacterium birth rate, and by C [bacteria] the carrying capacity of the population of intracellular bacteria within a single phagocyte, which represents limitation of nutrients necessary for replication, such as iron or tryptophan (Jones et al., 2012). The decision to assume logistic growth for the intracellular bacteria reflects the competition for resources within the phagocyte. The rate γ 1 is set to zero since only phagocytes experiencing an effective long-term bacterial infection (and within-phagocyte replication) are later considered in the within-host model. The initial state of the process X corresponds to the number of bacteria taken up by a phagocyte. Experimental evidence by Golovliov et al. (2003) suggests that the uptake of F. tularensis is relatively ineffective in monocytic cells so that, during the initial phase of the infection, on average only one or two intracellular bacteria per cell were observed. Thus, we assume that a single phagocyte will take up a single bacterium, hence the process X will always begin in state X(0) = 1. Once infected, the possibility of the phagocyte taking up more bacteria is neglected (Wood et al., 2014), and an increase in its bacterial load is solely due to the replication of the bacterium initially ingested. We refer the reader to the Supplementary Material where the impact of this assumption is further explored. The number of bacteria released upon rupture of an infected phagocyte will depend on the stochastic dynamics of this logistic growth process, as well as on the actual time when this rupture takes place. In order to describe this rupture event, we consider additional transitions to an absorbing (rupture) state, B, from any of the transient states in N, as shown in Figure 1B. The rate at which this rupture event occurs is assumed to be independent of the number of bacteria within the phagocyte. This is based on the fact that bacterial escape into the cytosol has been shown to be FIGURE 1 | Within-phagocyte model. (A) Logistic growth process for the within-phagocyte replication of bacteria; (B) logistic growth process with log-normally distributed phagocyte rupture, moving the process to absorbing state B; (C) approximation of the process in (B) by using a PH(η, T) distribution for the rupture time.
Frontiers in Microbiology | www.frontiersin.org both essential and sufficient for triggering caspase-3 activation, which is the mechanism thought to induce cell death (Santic et al., 2010). In fact, a recent experimental study by Brock and Parmely (2017) shows that cell death does not require high bacterial burden, nor does a large number of intracellular bacteria ensure that phagocyte rupture would result soon. This implies that X can be thought of as a stochastic birth-and-death process where t = 0 marks the start of a "clock" that counts up toward the time of rupture of the phagocyte. At this moment, X immediately transitions into state B, from whichever of the transient states this may be, this state accounting for the number of bacteria released upon rupture (i.e., the rupture size). By fitting a deterministic model to experimental data, Wood et al. (2014) found that the time T rupture taken for an infected phagocyte to rupture is log-normally distributed, T rupture ∼ logN(3.72, 0.385), so that the average rupture time is E[T rupture ] = 44.4 h. Instead of incorporating this log-normally distributed time in the within-phagocyte model, Wood et al. (2014) consider a deterministic logistic growth process for the amount of bacteria within the phagocyte. Finally, Wood et al. (2014) set the number of bacteria released to be equal to the amount of bacteria in this logistic growth process at time Median[T rupture ] hours (i.e., by considering Median[T rupture ] and neglecting the actual distribution of the random variable T rupture ), which leads to a constant and deterministic number of bacteria released for any infected phagocyte.
If a log-normal distribution of T rupture is used in our model to compute the probability distribution of the number of bacteria released upon phagocyte rupture, this leads to the process described in Figure 1B. However, by considering a log-normally distributed inter-event time in the stochastic process, the resulting process X in Figure 1B is no longer Markovian. In order to address this difficulty, we propose to approximate this log-normally distributed rupture time T rupture ∼ logN(3.72, 0.385) by a phase-type (PH) distribution, T rupture approx.
∼ PH(η, T), since the family of phase-type distributions is dense within the family of continuous nonnegative distributions (He, 2014). This leads to the process shown in Figure 1C. In the Supplementary Material, we explain in detail how one can approximate this log-normal distribution by an approximate phase-type distribution, which depends on parameters η (a vector) and T (a matrix). The resulting estimated parameters η and T obtained for a PH distribution which approximates the logN(3.72, 0.385) distribution, as well as a graphical representation of this approximation, are reported in Figure 2.
Once the log-normal distribution for the rupture time has been approximated by a PH distribution, the resulting withinphagocyte stochastic process X in Figure 1C is Markovian, and the probability distribution of the number of bacteria R released upon rupture can be analytically computed (see Supplementary Material). The probability distribution of R, defined in terms of the following probabilities R k = P(R = k) = probability that the infected phagocyte releases k bacteria upon rupture, is used in section 2.2 to incorporate inter-phagocyte variability (in the amount of bacteria released upon phagocyte rupture) in the within-host infection dynamics.

Within-Host Model
The within-host model proposed here is a birth-death-rupture model that replicates the dynamics of F. tularensis within the lung, following inhalation of some initial quantity of bacteria (initial dose), and is largely based on the original model by Wood et al. (2014). Within the lung, bacteria can be killed by host immune cells or ingested by host phagocytes. In the latter case, the phagocyte might kill the corresponding bacterium (e.g., if the phagocyte is activated), or this bacterium can escape the FCP and enter into the cytosol, resulting in rapid proliferation of the bacteria and the subsequent rupture and death of the phagocyte, as described by the within-phagocyte model. Three events are therefore included in the within-host model, as well as their effect on the total population of bacteria and the number of infected phagocytes, and are detailed as follows: • Phagocytosis and bacterial survival (rate α > 0 [hours −1 ]): this phagocytosis event refers to the phagocytosis and intracellular survival of a bacterium; that is, to phagocytosis resulting in bacterial escape from the FCP, and in the subsequent events represented by the within-phagocyte model.
• Extracellular bacterial death (rate µ > 0 [hours −1 ]): host defense mechanisms such as the complement system, antibodies, natural killer cells, activated phagocytes and antimicrobial peptides directly contribute to the killing of extracellular bacteria (Jones et al., 2012). These methods of killing, including phagocytosis with successful intracellular bacterial killing, are collectively represented in the within-host model as this single event, with rate µ.
• Rupture of infected phagocytes (rate δ = Median[T rupture ] −1 [hours −1 ]): following phagocytosis of bacteria that results in their survival and intracellular proliferation, infected phagocytes rupture and die. The distribution of the number of bacteria released, computed by means of the within-phagocyte model, is incorporated here in terms of probabilities R k . This then accounts for an addition to the original model by Wood et al. (2014), allowing for inter-phagocyte variability in the rupture size.
In this way, the within-phagocyte model in section 2.1 allows one to represent the intracellular bacterial dynamics for bacteria surviving the phagocytosis event, escaping the FCP and entering into the cytosol, eventually leading to phagocyte rupture and bacterial release. Phagocytosis leading to successful bacterial killing is one of the several mechanisms described above leading to bacterial death at the withinhost level. Furthermore, intracellular bacterial replication is not explicitly considered in the within-host model, since one bacterium is considered per infected phagocyte. Once rupture of an infected phagocyte occurs, the number of bacteria released to the extracellular environment is given by the rupture size distribution computed from the within-phagocyte model. Given that R k is the probability that an infected phagocyte, initially infected by a single bacterium, releases k bacteria on rupture, the rate at which an infected phagocyte ruptures releasing k bacteria in the within-host model is then given by δR k . We note that since ∞ k=1 R k = 1, δ can be interpreted as the total rate of rupture of a single phagocyte.
The within-host model can be described using a continuoustime two-dimensional Markov process Y = {Y(t) = (B(t), P(t)) : t ≥ 0}, where B(t) denotes the total number of extracellular bacteria and bacteria-containing phagocytes at time t ≥ 0, and P(t) represents the number of infected phagocytes at time t ≥ 0, B(t) ≥ P(t) for any time instant t ≥ 0. An initial state of Y given by Y(0) = (k, 0) represents that k is the number of bacteria initially inhaled by the individual (initial dose), and there are 0 infected phagocytes. When the total population of bacteria reaches a threshold M ∈ N, a response is assumed to occur and reflects the onset of symptoms in the infected individual (Wood et al., 2014). This state, M, referred to as the response state, is one of two absorbing states of Y; the other is state 0 and represents the clearance of infection without reaching this response threshold. A depiction of the model is provided in Figure 3.
Two summary statistics of interest in the within-host model are the probability of response and the mean response time. For each of these, an efficient analytic approach for their exact computation can be found in the Supplementary Material. In particular, we define π (i,j) as the probability of response given the initial state Y(0) = (i, j) (3) By applying first-step arguments, the following recursive formula for π (i,j) may be obtained with the boundary condition π (0,0) = 0 representing that the probability of response is equal to zero if the recovery state is reached. A detailed derivation of this expression, as well as an algorithmic solution to the previous equations, are provided in the Supplementary Material.
One may define the mean time until the infected host responds in terms of the onset of symptoms. This can be done by choosing a threshold in the total number of extracellular bacteria equal to M, and considering the time to get to M, Since absorption into the response state M is not certain, there is no guarantee that the time to response, T (i,j) , will be finite (i.e., T (i,j) = +∞ if the individual recovers without reaching the threshold state M, while T (i,j) < +∞ if this threshold is reached, leading to the onset of symptoms). Thus, one may compute the restricted mean response time, after which the conditioned mean response time can be obtained. That is, if T (i,j) denotes the time to reach state M provided that Y(0) = (i, j), then the restricted and conditioned mean response times are given respectively by where 1 A is equal to 1 if A is satisfied and 0 otherwise. Following a first-step analysis, a recursive formula for the scalar quantities r (i,j) is given by for 0 ≤ j ≤ i ≤ M − 1, with the boundary condition r (0,0) = 0 representing the restricted time to a response if the recovery state is reached. Similar arguments to those used for computing the dose response probabilities, and described in the Supplementary Material, may be used for solving Equation (6) in an algorithmic and matrix-oriented way.

Population-Level Model
With a multi-scale model of F. tularensis infection that captures both the intracellular and within-host dynamics, we can now formulate a population scale model. At the population level, outbreaks of tularemia, as a result of infection by F. tularensis, have been declared in microbiology laboratories (Shapiro and Schwartz, 2002). This is directly related to the fact that diagnosis of tularemia requires a high level of suspicion, since the disease often presents with non-specific symptoms and nonspecific results of routine laboratory tests (Report, 2008). Because F. tularensis is a risk to laboratory personnel, clinicians are required to notify the laboratory when tularemia is suspected in a given patient, and if this is not notified, an outbreak can occur within the laboratory when manipulating the contaminated samples, as in the outbreak reported by Shapiro and Schwartz (2002). In particular, this notification allows for manipulation of the corresponding samples to be carried out under strict control measures, such as Biosafety Level 2 (BSL-2) or BSL-3 conditions (Report, 2008). If proper control measures are not taken when manipulating these samples, or if an accident occurs, F. tularensis can be released to the air, triggering its airborne dispersal and spread. Specific high-risk sample manipulation activities that have been identified in the literature are centrifuging, grinding or vigorous shaking (Report, 2008).
Recent work has been carried out for the indoor airborne spread of pathogens while taking into account the ventilation regime in place at the facility under analysis, such as the airborne spread of bacteria in health care facilities (Liao et al., 2005). For these scenarios, zonal ventilation models that are able to link the airflow dynamics within the facility and the epidemic dynamics at the population level are extremely useful (Noakes and Sleigh, 2009;López-García et al., under review). We consider here the scenario of a laboratory consisting of two rooms joined by a corridor, and where a given ventilation setting, in terms of the airflow dynamics, is in place (see Figure 4). We consider that a fixed amount of bacteria is released at time t = 0 in a given room. This could represent some high-risk manipulation of a contaminated sample or any accident causing airborne release of F. tularensis bacteria, which we assume here passes unnoticed for the staff (Shapiro and Schwartz, 2002). Our aim  is to estimate, for any given ventilation setting and any possible spatial position of the release location within the laboratory, the total number of individuals who would develop symptoms in the near future.
We propose here to follow the approach introduced by Noakes and Sleigh (2009), recently extended by López-García et al., (under review), where a system of ordinary differential equations (ODEs) is used to model the concentration of F. tularensis in the air in the different spatial compartments of the laboratory. In particular, a ventilation regime is defined by splitting this laboratory into ventilation zones, where the main assumption is that the air is well-mixed in each zone, but that airflow imbalances across the different zones can lead to different pathogen concentrations in the air at each zone (Noakes and Sleigh, 2009). Therefore, individuals in the same ventilation zone are assumed to have equal probability of inhaling the F. tularensis bacteria. Airflow dynamics could be further refined by considering a larger amount of ventilation zones. If C i (t) [bacteria · m −3 ] denotes the concentration of bacteria in the air in zone i at time t, and p i (t) [bacteria] is the total amount of bacteria inhaled by each individual in this zone up to time t, then C i (t) and p i (t) satisfy the system of ODEs given in Figure 4. Here, V i [m 3 ] denotes the volume of zone i, Q i [m 3 · min −1 ] is the rate at which air is extracted from zone i by the ventilation system, β ij [m 3 · min −1 ] is the rate at which air flows from zone i to zone j, n i is the number of individuals in zone i, and ρ [m 3 · min −1 ] is the pulmonary rate, or the rate at which individuals inhale air (Noakes and Sleigh, 2009). We set n i = 2 for i ∈ {1, 2, 4, 5} to represent two individuals working in each of these zones during the bacterial release, where the propagation occurs in the timescale of minutes (see section 3), and n i = 0 for i ∈ {3, 6} (i.e., corridor areas).
We propose to link the dose response probabilities obtained from the within-host model with this zonal ventilation model, by considering that the steady state value of p i (t) is equal to the total dose that an individual in zone i inhales. Thus, we implicitly assume that the timescale at which p i (t) reaches equilibrium (minutes, see section 3), is significantly shorter than the timescale of the within-host infection dynamics (days, see section 3), so that lim t→+∞ p i (t) can be considered as the initial dose for individuals in zone i. We note that the differential equations in Figure 4 depend on the rates of the ventilation setting under analysis, and on the initial conditions C i (0), 1 ≤ i ≤ 6 (related to where the bacterial release occurs in the first place). In section 3, we consider different ventilation settings and potential initial locations of the bacterial release.

PARAMETER VALUES
In this section, we discuss how to calibrate our within-phagocyte and within-host models from data. We also consider different ventilation settings for the population model, according to values reported by Noakes and Sleigh (2009) for the airborne spread of bacteria within a health care facility.

Within-Phagocyte Model
In order to use the within-phagocyte model described in section 2 to compute the rupture size distribution of any given infected phagocyte, parameters λ and C must first be estimated for the logistic growth process in Figure 1 for the within-phagocyte bacterial replication. We do so, making use of experimental data of the number of intracellular bacteria within a phagocyte (Lindemann et al., 2011). In this experiment, measurements of the number of intracellular bacteria were only considered for phagocytes that were still alive and had not ruptured (Lindemann et al., 2011). Therefore, we obtain estimations for λ and C by calibrating the logistic growth process in Figure 1A, where rupture events are neglected.
A sequential Approximate Bayesian Computation (ABC) method is used to get estimations for these parameters. When implementing the ABC method, unknown parameter values are sampled from a prior distribution, and model predictions (e.g., number of intracellular bacteria at different time instants) are obtained for these parameter values. Once these predictions are in hand, one can compare these model predictions with experimental data by using a particular distance measure, and accept or reject these sampled parameter values depending on this distance being below or above a given threshold ǫ. Accepted sampled parameter values lead to a posterior distribution for the corresponding parameters (Kypraios et al., 2017).
We consider as prior distributions for each parameter λ ∼ U(0.01, 1) and C ∼ U(100, 1500), which have been set according to values previously estimated by Wood et al. (2014). We sequentially implement the ABC algorithm by considering successively smaller tolerances, ǫ, to refine the parameter space. For each pair (λ, C) of parameters sampled from the priors, the birth-and-death process is simulated using the Gillespie algorithm to obtain the number of intracellular bacteria as predicted by the model (Gillespie, 2007). Once this number is predicted from our model, these values are compared with data by Lindemann et al. (2011). In particular, if X(t) is the amount of bacteria predicted by our within-phagocyte model at time t, and Data(t) is the number of bacteria observed at that time instant according to data by Lindemann et al. (2011), which are available for a set of time instants T, we make use of the Euclidean distance d(Model Prediction, Data) = t∈T (X(t) − Data(t)) 2 1 2 , (7) so that each corresponding parameter pair (λ, C) is accepted only if d(Model Prediction, Data) < ǫ. At first the tolerance is set to ǫ = 100, so that an estimate of where the true parameters lie can be found. After this, the prior distributions are adjusted accordingly and the ABC algorithm is repeated for tolerances ǫ = 50, 25, 15, to converge around the posterior distribution (Toni et al., 2009). We note that threshold values ǫ = 100, 50, 25, 15 were chosen after a preliminary exploration of the parameter space and the corresponding distance measures between the model predictions and experimental measurements, so that a posterior sample of size 10 5 could be obtained in around 48 h, by using the high performance computing ARC3 facilities at the University of Leeds. A bivariate histogram of the sample posterior distribution obtained in this way is provided in Figure 5, with the median of the sample indicated. Univariate histograms for each parameter are given on the corresponding axes.

Within-Host Model
Estimated parameter values α and µ for the within-host model proposed by Wood et al. (2014) were obtained using non-linear least squares to fit their within-host model to experimental data for the number of extracellular bacteria within the host during the initial 48 h post infection. Since our within-host model is part of a multi-scale model which incorporates a variable number of bacteria released on rupture of any infected phagocyte, new estimations for these parameter values are now required. Thus, ABC is used to calibrate the parameters (α, µ) of the within-host model depicted in Figure 3 by using within-host infection data (Eigelsbach et al., 1962;White et al., 1964). We note that this requires the distribution of the number of bacteria released on rupture. This has been described in the Supplementary Material, using the posterior median values of λ and C of Figure 5. This same rupture distribution is used in each iteration of the ABC algorithm at the within-host level. In keeping with Wood et al. (2014), and to represent the heterogeneities at the population level in individual susceptibility, M is not fixed and is considered instead a random value M ∼ logN(26.2, 6.05), according to data by Saslaw et al. (1961) and Sawyer et al. (1966). These data report the amount of bacteria found within infected individuals at the time of symptoms onset. For small to moderate values of M, the exact analysis carried out in the Supplementary Material can be applied to compute the probability of response and the mean response time in the withinhost model. On the other hand, stochastic simulation approaches need to be implemented for large values of M. We note that given the potential extremely large values of M, the Gillespie algorithm is not a viable choice to simulate the within-host infection dynamics for these values, and an approximate τleaping procedure is used instead, with adaptive step size (Cao et al., 2006).
Prior distributions assumed for each parameter are α ∼ U(0, 1) and µ ∼ U(0, 25). Because of the shorter intervals considered in the priors of these parameters compared to those in section 3.1, we carry out here a standard rejection ABC (i.e., not sequential) where 2 × 10 5 iterations of the ABC algorithm were performed. Tolerance is set so that an acceptance rate of 1% is obtained, and a sample of size 2 × 10 3 is obtained for the posterior distributions. Due to the large orders of magnitude for the number of extracellular bacteria within the host observed in the data by Eigelsbach et al. (1962) and White et al. (1964), we propose here to use the Euclidean distance as for (λ, C) but over the logarithm of the predicted values and the observed data by Eigelsbach et al. (1962) and White et al. (1964). That is, we consider the distance The results of the ABC lead to the posterior bivariate histogram of Figure 6, which clearly indicates a positive correlation between parameters α and µ, where most of the learning occurs about the ratio µ/α. We note that this positive correlation is directly related to the fact that, intuitively, α and µ rates correspond  to within-host events which can be considered as opposite events in this system (one representing bacterial escape from the extracellular environment, facilitating disease, and the other representing bacterial death, preventing disease). Thus, our within-host model dynamics can replicate the experimental data by either considering that both events occur, simultaneously, at a slower or faster pace. However, we point out that since the (α, µ) joint distribution in Figure 6 (left) does not have the accepted sampled values homogeneously located all around the elliptic shape, where more accepted values can be found around the center of the ellipse than in the corners, one should consider that these parameter values (near the corresponding medians, given by the red circle) have larger posterior probability than the estimated values obtained by Wood et al. (2014) (red triangle). Final parameter values for the within-phagocyte and the within-host models are reported in Table 1.
We can compare our within-host model predictions, in terms of the number of bacteria throughout time, with the data by Eigelsbach et al. (1962) and White et al. (1964). In Figure 7 we plot the predictions made by our within-host model and compare them with the bacterial load data by Eigelsbach et al. (1962) and White et al. (1964), where the initial conditions are given as the corresponding data values at time t = 0. Similarly to results by Wood et al. (2014), our within-host model does better in predicting the data by White et al. (1964), where larger amounts of bacteria were measured within the host.   Eigelsbach et al. (1962) and White et al. (1964), vs. data points by Eigelsbach et al. (1962) and White et al. (1964). Median values of α and µ considered as computed in Figure 6.

Population-Level Model
Four different scenarios A1, A2, C1, and C2 are considered depending on two potential bacterial release locations (see Figure 4). Two potential ventilation regimes (A and C) within the microbiology laboratory have been chosen, as described in Table 2: ventilation regime A (scenarios A1 and A2) and ventilation regime C (scenarios C1 and C2) considered by Noakes and Sleigh (2009) and López-García et al. (under review). Regardless of the particular location where it occurs, it is assumed that 10 5 bacterial counts are released at time t = 0. In each scenario it is assumed that V i = 36m 3 for i ∈ {1, 2, 4, 5} and V i = 12m 3 for i ∈ {3, 6}. The pulmonary rate is set to ρ = 0.01m 3 · min −1 (Noakes and Sleigh, 2009), while the remaining parameters in Figure 4 are provided in Table 2, along with the steady state values p (k) = lim t→∞ p (k) 1 (t), p (k) 2 (t), p (k) 4 (t), p (k) 5 (t) , k ∈ {A1, A2, C1, C2}. A graphical representation of scenarios A1, A2, C1 and C2 is given in Figure 8, and the time course of the variables C i (t), 1 ≤ i ≤ 6, and p j (t), j ∈ {1, 2, 4, 5}, are plotted for scenario A1 in Figure 9 for illustrative purposes.

RESULTS
The distribution of the number of bacteria released by an infected phagocyte, for posterior median values of λ and C from Figure 5, is provided in Figure 10. In order to compare with results by Wood et al. (2014), let us note that the approach they use involves evaluating a deterministic logistic growth process at the median (log-normally distributed) time taken for an infected phagocyte to rupture. The method here may instead be interpreted as computing the distribution of the number of bacteria generated by means of the analogous stochastic logistic growth process, but when the actual log-normally distributed rupture time is incorporated into the model (see Figures 1B,C).
Since the deterministic and stochastic processes have both been parameterized using the same data set, they are comparable, and the median number of bacteria released from our predicted distribution in Figure 10 is approximately equal to the fixed value of 358 bacteria released upon rupture estimated by Wood et al. (2014), supporting the fact that the median number of bacteria released had previously been estimated correctly. Despite this, the method outlined here is more general, since it allows to incorporate the log-normal distribution of rupture times, and thus, a more comprehensive analysis of the number of bacteria released can be conducted, and incorporated into the within-host dynamics, by considering inter-phagocyte rupture size variability. Moreover, we note that the mean number of bacteria released on rupture is predicted to be 288, significantly lower than the fixed value 358 considered by Wood et al. (2014). This is directly related to the bimodal shape of our predicted rupture size distribution, which suggests that some phagocytes will likely rupture with just a few bacteria, and that the total number of bacteria released by each single infected phagocyte was slightly over-estimated by Wood et al. (2014) on average. We note that our model is able to predict that a significant amount of phagocytes might rupture releasing few bacteria, which is something that the deterministic approach followed by Wood et al. (2014) does not reflect. We also note that the actual rupture size distribution, to the best of our knowledge, has not been experimentally measured in vitro yet, which would allow us to do model selection based on predictions in Figure 10. However, it has been recently experimentally observed by single-cell analysis (Brock and Parmely, 2017) that a significant amount of phagocytes can die releasing very few bacteria. While the deterministic amount of bacteria proposed by Wood et al. (2014) cannot account for this, our model predicts indeed a significant amount of phagocytes releasing very few bacteria, which is represented by the first mode in Figure 10. This suggests that this mode is not an artifact caused by the stochastic within-host model, but that phagocytes rupturing soon (according to the estimated log-normally distributed rupture time) would not have enough  = (145, 82, 13, 17) = β 56 = β 45 = β 21 = β 32 = β 65 = β 54 = 9 A2 β 12 = β 23 = β 36 = β 63 Q i = 3, i = 1, .., 6 5 p (A2) = (17, 23, 82, 110) = β 56 = β 45 = β 21 = β 32 102,46,9,9) = β 56 = β 45 = 9 Q 2 = Q 3 18,18,92,92) = β 56 = β 45 = 9 Q 2 = Q 3 Airflow parameters for the four scenarios considered, and steady state bacterial intake values representing initial dose for individuals at each zone. Airflow parameters have been chosen according to those in the ventilation regimes A and C considered by Noakes and Sleigh (2009) and López-García et al. (under review).
FIGURE 8 | Ventilation scenarios considered in the microbiology laboratory. Four scenarios A1, A2, C1, and C2 corresponding to two potential release locations (zone 1, scenarios A1 and C1; zone 5, scenarios A2 and C2). Ventilation regime in scenarios A1 and A2 represents a well-mixed ventilation, where airflow (arrows, with β ik rates given as red numbers) is well balanced across zones and same extract ventilation (circled values) is considered in all zones. Ventilation regime in scenarios C1 and C2 represents airflow occurring from the corridor areas to the opposed side of the rooms, where extract ventilation is in place.
time for substantial bacterial proliferation, leading to small rupture sizes predicted by the model and being experimentally observed.
By using this rupture size distribution, and the within-host model in section 2.2, the probability of response and mean response times can be computed for varying initial doses. In Figure 11 (left), we plot the cumulative probability of response (i.e., cumulative probability of the process in Figure 3 reaching state M), as predicted from our model for different initial doses. We note that the asymptotic values in Figure 11 (left) represent the probabilities of response for each initial dose. We plot in Figure 11 (right) the (conditioned) mean time until response predicted for different initial doses, and compare this with the predictions by Wood et al. (2014) and with data of the time until symptoms onset observed in infected individuals (Saslaw et al., 1961;Sawyer et al., 1966). Our predictions are obtained by using the posterior median parameter values in Figures 5, 6. We note that, once parameters (α, µ) are estimated as explained in section 3.1, results obtained here for the probability of response and the (conditioned) mean response time are very similar to those previously found by Wood et al. (2014), indicating that the multi-scale model is not only capable of reproducing their results, but also corresponds well with the two experimental data sets by Saslaw et al. (1961) and Sawyer et al. (1966). However, we note that our multi-scale model only replicates well these results for posterior distribution of (α, µ) in Figure 6, where our predicted median values are far away from those parameters estimated by Wood et al. (2014). In particular, although these parameters are highly correlated and determining their individual true values is difficult, the histogram in Figure 6 suggests that the ratio of α and µ ranges from 24.69 to 27.54, which is lower than the ratio of 31.95 found by Wood et al. (2014). Moreover, our results in Figure 6 suggest that both α and µ were underestimated by Wood et al. (2014) (see the red circle and triangle in Figure 6).
At the population level, one can use the probability of response for each individual computed from the within-host model, where their initial dose is given by the steady state values in Table 2, in FIGURE 9 | Predicted airborne spread and inhalation of bacteria in the laboratory. Time course of the variables C i (t), 1 ≤ i ≤ 6, and p j (t), j ∈ {1, 2, 4, 5}, for scenario A1.
FIGURE 10 | Predicted rupture size distribution. The distribution of the predicted number of bacteria released by a single phagocyte on rupture, as computed from the within-phagocyte model, compared to the fixed value assumed by Wood et al. (2014). The posterior median values for λ and C in Figure 5 are used to compute this distribution.
order to compute the distribution of the number Z of individuals within the laboratory showing symptoms after the bacterial release, for each of the four scenarios considered in Table 2. These distributions are plotted in Figure 12, together with the corresponding expected values E [Z]. From this, it can be seen that scenarios associated with smaller number of responses are A1 and C1, that is, when the bacteria are released from zone 1 as opposed to zone 5. This might be expected since air can flow from zone 5 into other areas more easily, whereas it only flows into one other zone from zone 1. However, an interplay between the ventilation regime (i.e., airflow dynamics) and the bacterial release location can be observed, where the ventilation regime in scenario C1 helps to decrease pathogen concentration in the release zone (zone 1), due to significant extract ventilation in place in this zone, while this same ventilation implies in scenario C2 the airborne spread of pathogen from zone 5 toward zone 4, causing more infections at the population level.

DISCUSSION
In this work, we propose a multi-scale model for the infection dynamics of F. tularensis which covers the within-phagocyte, within-host and population scales. The within-host model should be considered an extension of the model originally proposed by Wood et al. (2014), where inter-phagocyte rupture size variability is incorporated in the distribution of the number of bacteria released upon rupture by any infected phagocyte. This distribution is computed by means of a stochastic logistic growth process for the replication of bacteria at the within-phagocyte level, but where the log-normally distributed rupture time predicted by Wood et al. (2014) is explicitly incorporated by means of a PH-type approximation. This approximation allows us to consider a Markovian stochastic process for the within-phagocyte infection dynamics. Once the extended within-host model is set up, we provide analytical approaches for computing the probability of response (in terms of the number of extracellular bacteria within the host to reach some response threshold M), and the mean time until this response takes place (conditioned on this response actually occurring). By calibrating the within-phagocyte and within-host model parameters using experimental infection data, our multi-scale model predictions are in agreement with experimental data both at the within-phagocyte and within-host level.
The main advantages of our multi-scale model are: • The within-phagocyte model allows us to incorporate the estimated log-normally distributed rupture time into the bacterial proliferation dynamics, while keeping the Markovian nature of the original process. This allows the exact distribution of the rupture size to be computed.   Figure 8, for scenarios A1, A2, C1, and C2. That is, probabilities P(Z = z), 0 ≤ z ≤ 8.
We believe that our methodology, using phase-type approximations for incorporating non-Markovian events in these intracellular processes, as well as the first-step arguments considered here for computing the rupture size distribution, is applicable to other intracellular bacterial replication systems.
• The rupture size distribution computed by our model and plotted in Figure 10 is able to capture the fact that a significant amount of phagocytes might die releasing very few bacteria, which has been recently experimentally observed (Brock and Parmely, 2017). • The stochastic nature of the within-phagocyte model incorporates inter-phagocyte variability in the rupture size in the within-host model, relaxing the assumption made by Wood et al. (2014) that every phagocyte releases a fixed amount of bacteria. Relaxing this assumption leads to different predictions in the posterior estimated values of within-host parameters (α, µ), as shown in Figure 6, with respect to previous predictions made by Wood et al. (2014). This is directly related to the fact that different behaviors can be expected when the within-host model is simulated with the actual rupture distribution (so that each phagocyte, upon rupture, can release different numbers of bacteria with different probabilities) instead of considering that every phagocyte releases a fixed number of bacteria, even if this fixed release is set equal to the median value of the distribution in Figure 10. • The zonal ventilation model is a simple but flexible way of representing airborne spread of bacteria, and of linking this spread with the initial doses infecting each of the individuals in the laboratory under study. Our results suggest that there is a clear interplay between the potential release location and the ventilation in place within the laboratory, where an appropriate ventilation regime can decrease the number of individuals developing symptoms.
The original model by Wood et al. (2014), as well as the extended model proposed here, should be considered as one of the few and recent attempts to propose mechanistic models for the computation of dose-response probabilities and the mean time until individuals showing symptoms following bacterial infection. Many of the original approaches in the literature to this aim usually involve adjusting exponential and beta-Poisson models to data (Chen, 2007;Huang and Haas, 2009). These models are limited since the real within-host biological mechanisms at play are not explicitly considered, and the distributions are selected only due to their ability to approximate the experimental or clinical data. Moreover, timescales for the different within-host processes are usually not explicitly considered in these models, where the final output of the model is usually limited to the dose-response probability curve. Thus, recent attempts are being made in order to explicitly consider the biological mechanisms following bacterial infection, leading to computational models which can analyse the timescales of these intracellular and within-host processes, not only for F. tularensis but also for other pathogens such as anthrax (Day et al., 2011).
Developing new mathematical and computational models that can explicitly account for biological mechanisms requires a significant amount of quantitative experimental data, and a balance between model complexity and experimental information must always be struck. For example, in our within-host model, all the mechanisms leading to extracellular bacterial death, such as the complement system, antibodies, natural killer cells, antimicrobial peptides or phagocytosis leading to bacterial killing are represented as a single event occurring at rate µ. If one were to distinguish all of these events in the model, experimental measurements of the specific contribution of each mechanism would be required, and a new version of our multi-scale model could be proposed. An additional limitation of our model, at the within-phagocyte level, is the fact that the rupture time is modeled as a lognormally distributed time which is independent of the bacterial proliferation dynamics simultaneously occurring within the phagocyte. Ideally, if we had enough experimental knowledge about the effect that the bacterial load has on the rupture of the phagocyte, one could consider that the rate of rupture from any state n in Figure 1 (i.e., n bacteria within the phagocyte at a given time) is a function δ n of this bacterial load. Thus, using the independent log-normally distributed time estimated by Wood et al. (2014) should be seen as a compromise between current experimental knowledge and model complexity, and is based on the fact that bacterial escape into the cytosol has been shown to be both essential and sufficient for triggering caspase-3 activation, which is the mechanism thought to induce cell death (Santic et al., 2010). This also agrees well with recent experimental evidence (Brock and Parmely, 2017) showing that cell death does not require high bacterial burden, nor does a large number of intracellular bacteria ensure immediate phagocyte rupture. Finally, at the population-level, we note that more elaborated fluid dynamics simulations could be considered for the airborne spread of F. tularensis in the microbiology laboratory. We propose here a zonal ventilation model as a simple but flexible way of linking the indoor airflow dynamics with the initial dose of each individual after a bacterial release. We note however that the imprecisions inherently caused by the spatial discretisation in this zonal ventilation approach, where the indoor setting is split in a number of zones and the air is assumed to be well-mixed within each zone, can be reduced by increasing the number of zones under consideration.
The development of a mathematical model of infection dynamics at different scales is a challenging problem for which few successful attempts have been made in the literature so far (Bauer et al., 2009). To the best of our knowledge, this is the first multi-scale model for F. tularensis trying to account for the infection dynamics from the intracellular to the population level. It is conceivable that the future of in silico modeling will consist of a large number of interconnected models at different scales, and where one of the main aims will be to predict the effects that perturbations of model parameters along the different scales can have in the global infection dynamics. Finally, the approach presented in this article could also be readily applied to investigate the potential casualty impacts resulting from a deliberate bioterrorism or biological warfare attack in civilian and military scenarios. For instance, our multi-scale model may be used in conjunction with a larger-scale outdoor dispersion model that produces F. tularensis concentration estimates over large areas of terrain.