Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Stochastic dynamics of Francisella tularensis infection and replication

Abstract

We study the pathogenesis of Francisella tularensis infection with an experimental mouse model, agent-based computation and mathematical analysis. Following inhalational exposure to Francisella tularensis SCHU S4, a small initial number of bacteria enter lung host cells and proliferate inside them, eventually destroying the host cell and releasing numerous copies that infect other cells. Our analysis of disease progression is based on a stochastic model of a population of infectious agents inside one host cell, extending the birth-and-death process by the occurrence of catastrophes: cell rupture events that affect all bacteria in a cell simultaneously. Closed expressions are obtained for the survival function of an infected cell, the number of bacteria released as a function of time after infection, and the total bacterial load. We compare our mathematical analysis with the results of agent-based computation and, making use of approximate Bayesian statistical inference, with experimental measurements carried out after murine aerosol infection with the virulent SCHU S4 strain of the bacterium Francisella tularensis, that infects alveolar macrophages. The posterior distribution of the rate of replication of intracellular bacteria is consistent with the estimate that the time between rounds of bacterial division is less than 6 hours in vivo.

Author summary

Infecting a host cell is required for the replication of many types of bacteria and viruses. In some cases, infected cells release new infectious agents continuously over their lifetime. In others, such as the Francisella tularensis bacterium studied here, they are released in a single burst that coincides with the cell’s death. We show how a stochastic model, the birth-and-death process with catastrophe, can be used to characterise infection in a single cell, thereby allowing us to account for burst events and quantify the kinetics of pathogenesis in the lung, the initial site of infection, as well as in other organs that the infection spreads to. We learn about the parameters of the mathematical model of Francisella tularensis infection making use of the experimental measurements of bacterial loads, together with approximate Bayesian statistical inference methods. The most important parameter describing the pathogenesis is the rate of replication of intracellular bacteria.

Citation: Carruthers J, Lythe G, López-García M, Gillard J, Laws TR, Lukaszewski R, et al. (2020) Stochastic dynamics of Francisella tularensis infection and replication. PLoS Comput Biol 16(6): e1007752. https://doi.org/10.1371/journal.pcbi.1007752

Editor: James O’Dwyer, University of Illinois at Urbana-Champaign, UNITED STATES

Received: July 17, 2019; Accepted: February 27, 2020; Published: June 1, 2020

Copyright: © 2020 Carruthers et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Computer codes (in Python) to generate the numerical realisations of the agent-based model and to perform the cohort analysis are available in this link http://archive.researchdata.leeds.ac.uk/677/.

Funding: The experimental work shown here was funded by the UK Ministry of Defence. Mathematical research was supported by the Medical Research Council through the Skills Development Fellowship number MR/N014855/1 MLG. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Francisella tularensis, the causative agent of tularemia, is extremely infectious and considered a biothreat agent [1–3]. Treatment options are limited: a live attenuated vaccine exists but is not in mainstream use [4, 5]. Protection via antibiotics is dependent on early diagnosis and timely administration [6]. Francisella tularensis bacteria may be inhaled in an aerosol, with initial doses as low as ten colony-forming units (CFU) resulting in respiratory or pneumonic tularemia [7–10]. The bacteria enter alveolar macrophages [11–16], evading initial immune recognition and inflammatory response because of their atypical lipopolysaccharide [17]. They are able to escape from phagosomes in less than an hour and, as illustrated in Fig 1, begin multiple rounds of replication in the cytosol [18–21]. Instead of producing inflammatory cytokines, the first infected macrophages produce anti-inflammatory TGF-β cytokine. The eventual death of the host macrophage [22] returns bacteria to the extracellular environment, from where they can migrate to another organ, or again infect macrophages in the lung.

thumbnail
Download:
Fig 1. Within-host model of Francisella tularensis pathogenesis.

A single bacterium is taken up by a macrophage, inside a phagosome (top left). Without activating the macrophage, the bacterium escapes and multiplies inside the cytosol (top line) eventually causing the macrophage to rupture and release many bacteria. Free bacteria (center) may infect other macrophages, die or migrate to a different organ (bottom line). Macrophages exist in resting (green), anti-inflammatory (blue) and pro-inflammatory (red) states.

https://doi.org/10.1371/journal.pcbi.1007752.g001

Our agent-based computational model of the first 72 hours after infection is based on that of Gillard et al. [23], who considered the first 24 hours after infection. Classical mathematical models of intracellular infection consider variables describing populations of uninfected cells, infected cells and free infectious particles [24–26]. In such models, the rate of production of new infected cells is assumed to be proportional to the number of infected cells, which is true if each infected cell, independently, releases infectious particles at a constant rate. It is possible to go beyond the simplest hypothesis by considering subpopulations of infected cells: in an “eclipse” phase or productive phase [27, 28], or considering different multiplicities of infection and co-infection [29]. In this work, we seek to describe an scenario where bacteria continue to divide inside the host cell, without any being released from the cell, until the host cell ultimately ruptures, typically releasing more than a hundred bacteria at a time [23, 30–33]. Our mathematical approach is a stochastic model of the population of infectious agents inside one host cell, extending the simple birth-and-death process by the occurrence of catastrophes.

A consequence of the extreme virulence of Francisella tularensis is that initial doses of bacteria used in experiments are small enough for it to be reasonable to assume that host macrophages are infected by only one bacterium each [18]. Thus, the ensemble of realisations of the stochastic process describing the dynamics of a population of bacteria inside one cell can be thought of as describing the dynamics inside a set of host cells, that behave independently until they rupture.

In this paper we model the number of bacteria inside an infected macrophage, making use of a birth-and-death process with catastrophe [34], and the bacterial populations in multiple organs in the first few days after infection. Expressions for these variables are first computed for the lung before considering the mesenteric lymph nodes (MLNs), liver, kidney and spleen. We describe host cell rupture by a load-dependent “hazard rate”: an infected macrophage’s rupture probability per unit time is proportional to its bacterial load. Thus, we assume that cells with high bacterial load at a given time are more likely to rupture than those with a lower bacterial load, but there is no fixed maximum or minimum load [35]. Our assumption is consistent with observations of infection, and apoptosis, of murine macrophage-like J774.A1 cells [36]. Alternatively, a mathematical modeller may take the rupture time and the number of bacteria released per rupture event to be fixed parameters that can be inferred from experimental data [30], or assume a distribution of rupture times that is independent of the intracellular dynamics [33]. Similar questions arise in the modelling of Salmonella enterica infections [37–39]. In this work, we calculate the distributions of rupture times and number of bacteria released per rupture event as a consequence of the stochastic description of an infected macrophage and its bacterial contents.

We make use of the Sobol method of global sensitivity analysis to identify which parameters, in the mathematical model that describes the early days of infection, have the greatest effect on bacterial counts in each organ [40, 41]. Using a decomposition of the variance, this approach allows us to see how the variance in bacterial counts changes when combinations of parameters are fixed. If fixing a parameter results in a large reduction in the variance, this parameter will be of greater importance. With an approximate Bayesian computation algorithm [42], we learn about the most important parameters by comparing predictions of the model with the bacterial counts measured in the experiments.

Results

Birth-death-catastrophe process

In order to describe the course of Francisella tularensis infection for a given host, we need to consider the dynamics inside one infected cell, beginning at the time of entry of one infectious agent into the cytosol and ending either with the rupture of the cell and release of its bacterial content, or with the elimination of the infection from the cell. To describe the dynamics inside a host cell, we assume the most simple hypothesis; namely, that each bacterium has a constant division and death rate. The corresponding rates for the bacterial population inside a single host cell are proportional to its size. Since the probability per unit time that the host cell ruptures is a function of its intracellular bacterial population size, the independence property required for a branching process description does not hold in this case. Rupture of the host cell is a “catastrophe” event that affects every bacterium in the cell at the same time. Catastrophes have been considered mathematically as an extension of birth-and-death processes [43, 44], including scenarios where a subset of the population is removed [45–51], or where a catastrophic event kills the entire population [52–57]. Our interest here is in the process depicted in Fig 2, where two distinct absorbing states exist, both of which represent the loss of all intracellular bacteria. The first, state 0, represents the recovery of the macrophage due to successful elimination of its bacterial load. The second, state H, represents the rupture of the macrophage and the release of its entire content of bacteria [34]. Thus, the number of bacteria inside the cytosol of a single infected cell is denoted Xt where In this Section, we take X0 = 1. That is, we assume that the cell only phagocytoses one bacterium, and t = 0 is the time at which the bacterium escapes from the phagosome. We also ignore, for now, the possibility of reinfection. Macrophages containing any number of bacteria may rupture; those with high bacterial load are more likely to do so than those with low bacterial loads [58]. Thus, we assume that the rate associated with the catastrophe event is proportional to the instantaneous number of bacteria: if a macrophage contains Xt bacteria at time t, the probability that it ruptures before time t + Δt is δXtΔt, in which case all of the Xt bacteria are released from the cell.

thumbnail
Download:
Fig 2. A birth-and-death process with catastrophe representing division, death and rupture.

The state n represents a macrophage with n cytosolic bacteria. There are three types of events: transition to state n + 1 (division of a bacterium, rate βn), transition to state n-1 (death of a bacterium, rate μn), and transition to state H (rupture of the macrophage with release of n bacteria, rate δn). In this work we assume that βn = βn, μn = μn and δn = δn. The states 0 and H are absorbing. The infected macrophage survives for as long as it does not reach the state H.

https://doi.org/10.1371/journal.pcbi.1007752.g002

We assume that Francisella tularensis bacteria replicate in the cytosol of their host macrophages with rate β per bacterium, and are susceptible to intracellular death through misfortune or cellular defence mechanisms with rate μ. Estimates based on observations of Francisella tularensis proliferation suggest that β ≃ 0.16 h−1 and μ ≪ β [23, 59]. As with cell rupture, the probability that an intracellular bacterium divides or dies in a short interval (t, t + Δt) is proportional to Xt. Thus, if Xt is the bacterial load at time t, then the bacterial load at time t + Δt is either

Dynamics of bacterial load in a single infected cell

The first quantity of interest from the perspective of an infected macrophage is its survival function, S(t): the probability that a single infected macrophage has not ruptured by time t. S(t) is also the fraction of macrophages with a single bacterium in their cytosol at time t = 0 that survive to time t. We can think of this as the surviving fraction in a large group of macrophages, each initially infected with a single bacterium. We can write If a macrophage is carrying Xt bacteria in its cytosol at time t, its probability of rupture between t and t + Δt is equal to δXtΔt. The function S(t) is the average over realisations, so . In other words, S(t) satisfies the following differential equation (1) One way to calculate the survival function is to find the distribution of Xt and make use of (1). Before proceeding to the explicit calculation of Xt, however, it is useful to consider a direct method for calculating S(t), using an extended definition of the survival function. Let S(k)(t) be the survival function of a single macrophage with k cytosolic bacteria at t = 0: If X0 = 1 then, as Δt → 0, either

  • XΔt = 1 with probability 1 − (β + μ + δ)Δt,
  • XΔt = 0 with probability μΔt,
  • XΔt = 2 with probability βΔt,
  • XΔt = H with probability δΔt.

Thus, we have and (2) If k = 2, then the number of bacteria at time t can be written as the sum of two families: the first initial bacterium with its progeny and the second initial bacterium with its progeny (see Fig 3). The probability of rupture between t and t + Δt is given by , with S(1)(t) = S(t). That is We therefore conclude that S(2)(t) = [S(t)]2. Similarly, the probability that a macrophage, initially infected by k bacteria, is alive at time t is equal to the probability that k independent macrophages, each infected by a single bacterium, are all alive at time t: (3) Thus, from (2), the survival function S(t) obeys (4)

thumbnail
Download:
Fig 3. Realisations of the birth-death process with catastrophe.

On the left, the initial number of bacteria, k = 1. On the right, the initial number of bacteria, k = 2. On the right, the two families follow, independently, the same stochastic process as in the case k = 1. However, the catastrophe affects both families at the same instant. The parameter values are β = 0.15 h−1, μ = 0.01 h−1 and δ = 0.01 h−1 [23].

https://doi.org/10.1371/journal.pcbi.1007752.g003

The solution of (4) agrees with that obtained by Karlin and Tavare [34] (5) where a and b are the zeros of the function F(S) = βS2 − (β + μ + δ)S + μ; that is, and (see Fig 4). The survival function itself is shown in Fig 5. Note that S(t)→a as t → + ∞, so that a is equal to the probability that the infected macrophage eliminates the infection, as opposed to rupturing. If μ = 0 then a = 0, and the survival function takes the simpler form,

thumbnail
Download:
Fig 4. Polynomial governing the survival function.

F(S) = βS2 − (β + μ + δ)S + μ is the RHS of (4) with parameter values β = 0.15 h−1, μ = 0.01 h−1 and δ = 0.01 h−1. The constants a and b satisfy and . The value of a is equal to the probability that the infected macrophage eliminates the infection, rather than ruptures.

https://doi.org/10.1371/journal.pcbi.1007752.g004

thumbnail
Download:
Fig 5. Survival function, probability density function and bacterial release on rupture.

Left: the macrophage survival probability function, S(t). Centre: the bacterial load is proportional to . The vertical line is (6). Right: the function that gives the mean number of bacteria released per macrophage. The parameter values, taken from Ref. [23], are β = 0.15 h−1, μ = 0.01 h−1 and δ = 0.001 h−1. The dashed lines show the corresponding functions when μ = 0.

https://doi.org/10.1371/journal.pcbi.1007752.g005

The second quantity of interest is f(t), the probability density function of the time until rupture of a macrophage initially infected with a single bacterium, given by

The maximum value of f(t) is found at (see Fig 5) (6) If μ = 0 then (7) The third function of interest provides more detail about the kinetics of pathogenesis: the rate of release of bacteria, per infected macrophage, as a function of time. The probability that a macrophage, containing Xt bacteria at time t, ruptures before t + Δt is δXtΔt, in which case the number of bacteria released is simply Xt. The mean number of bacteria released between t and t + Δt is therefore r(t)Δt where (8) We note that The most elegant way to evaluate moments of Xt is making use of the probability generating function: (9) In Materials and methods, we show that (10) where Note that G(k)(1, t) = S(k)(t) [34]. When μ = 0, we have In the case k = 1, we make use of (10) to obtain and the probability that there are n bacteria at time t is which is a geometric distribution. The function r(t), defined in (8), can be written as the product : (11) Thus, , illustrated in Figs 5 and 6, is interpreted as the mean number of bacteria released by a macrophage, given that it ruptures at time t. Note that .

thumbnail
Download:
Fig 6. Agent-based realisation compared to predicted means: First cohort of macrophages.

In a numerical realisation, N = 30 macrophages are infected, by one bacterium each, at t = 0. The red line shows the number of those macrophages surviving up to time t. The dotted red curve is NS(t), using the survival function (5). Each blue dot in the lower panel coincides with a downward step in the red line, corresponding to a macrophage rupture event. The dotted blue curve is , given by (11). Parameter values have been taken from Ref. [23]: β = 0.15 h−1, μ = 0 and δ = 0.001 h−1. Agent-based realisations are simulated making use of tau-leaping time-stepping with Δt = 0.01 h−1, M = 104, ρ = 0.01 h−1, ϕ = 2 h−1, μE = 0.01 h−1 and γ = 1 h−1.

https://doi.org/10.1371/journal.pcbi.1007752.g006

Let K be the random variable denoting the number of bacteria released from a single infected macrophage when it ruptures. If μ = 0 then [23]. If μ ≠ 0 but Xt is ultimately absorbed into the catastrophe state, the mean number of bacteria released is b/(b − 1) and [34]. Karlin and Tavare analysed the processes obtained by conditioning on the ultimate fates: elimination or catastrophe [34]. Taking both possibilities into account, the mean number of bacteria released is

Cohort analysis

Bacterial loads can be measured in different organs of a given infected mouse, but only at one time point. In an agent-based simulation, on the other hand, the entire history of every macrophage and bacterium is available. We classify bacteria into cohorts, according to their lineage, assigning a “cohort number” attribute to each bacterium as follows. At the start of the realisation, the cohort number of every bacterium is equal to zero. The cohort number increases by one whenever a bacterium enters a macrophage. When bacteria divide, the daughters inherit their cohort number from their mother. Thus, for the initial dose of bacteria, each bacterium has a cohort number equal to zero. Following their uptake by macrophages, each of their cohort numbers increases to one. A realisation of first cohort macrophage rupture events is shown in Fig 6.

In order to calculate the total number of bacteria in a particular organ, we consider cohorts of bacteria contained within macrophage phagosomes and cytosols. We define the quantities

  • Pn(t), the mean number of cohort n bacteria in macrophage phagosomes at time t ≥ 0, and
  • Cn(t), the mean number of cohort n bacteria in macrophage cytosols at time t ≥ 0.

The initial condition is P1(0) = N. That is, we assume phagocytosis of the initial dose of bacteria occurs instantaneously. Bacteria enter the phagosome from a previous cohort rupture event and escape to the cytosol with rate ϕ. The bacteria inside the cytosol then replicate with rate β before being released in a cohort rupture event. In the calculations of this section, we assume that all bacteria released in macrophage rupture events are immediately absorbed by uninfected macrophages. This assumption is consistent with the dynamics of the agent-based model, as long as the supply of uninfected macrophages is much larger than the number of extracellular (or free) bacteria.

Given the per bacterium rate of phagosomal escape, ϕ, the mean number of first cohort bacteria in phagosomes is simply (12) In the agent-based computation, the bacteria escape from the phagosome to the cytosol at a different time in each macrophage, drawn from an exponential distribution with mean 1/ϕ. Accounting for this delay, the mean number of first cohort bacteria in macrophage cytosols is, thus, equal to Here, f(t)/δ is the mean of the birth-death-catastrophe process considered in the previous section, , where t = 0 was taken to be the time of phagosomal escape of the initial bacterium to the cytosol. The mean rate of release of first cohort bacteria from rupturing macrophages at time t after the start of the experiment, is r1(t) where (13) For the second cohort, the functions P2 and C2 satisfy (14) and (15) If we define then, in general, higher order cohorts of phagosomal and cytosolic bacteria satisfy (16a) (16b) For the initial 48 hours of Francisella tularensis infection, it is sufficient to only consider the first three cohorts of bacteria, or equivalently, first and second order cohort rupture events. A comparison between the mean of 102 simulations of the agent-based model and these approximations is provided in Fig 7. With the cohorts of bacteria in macrophage phagosomes and cytosols, the total number of intracellular bacteria in the lung at time t is found by summing the total number in each cohort.

thumbnail
Download:
Fig 7. Cohorts of bacteria in the lung.

The formulæ for the number of bacteria in phagosomes (left) and cytosols (right), obtained from (16a) and (16b) are shown as dashed lines. Averages over 102 realisations of the agent-based computational model are shown as solid lines. The parameter values, taken from Ref. [23], are β = 0.15 h−1, μ = 0, δ = 0.001 h−1, M = 104, ρ = 0.01 h−1, ϕ = 2 h−1, γ = 0.1 h−1, μE = 0.01 h−1, N = 102 and Δt = 0.01h−1.

https://doi.org/10.1371/journal.pcbi.1007752.g007

In order to determine the mean number of bacteria in the MLNs, liver, kidney and spleen, the mean number of extracellular bacteria in the lung must first be calculated. Bacteria released through rupture events either reinfect macrophages in the lung with rate Mρ(bacteria h)−1, are killed extracellularly with rate μE, or migrate to a different organ with rate γ. Here, γ is the per bacterium rate of exiting the lung. The destination of a migrating bacterium is then determined by weights assigned to each organ. These weights are given by wj for , and satisfy ∑j wj = 1 [60]. The rate at which a bacterium migrates from the lung to the liver, say, is then γwliver. If E(t) denotes the mean number of extracellular bacteria in the lung at time t ≥ 0, then this variable satisfies where M is the number of macrophages within the lung capable of taking up bacteria.

In the agent-based model, the dynamics in the other organs is equivalent to that in the lung. However, since all bacteria are intracellular until first cohort rupture events occur, and are likely to again infect macrophages in the lung following these events, the mean number of extracellular bacteria is small during the early stages of infection. Assuming that bacteria migrating away from the lung are quickly phagocytosed upon reaching their destination, the mean number of bacteria contained within macrophage phagosomes and cytosols in organs, aside from the lung, then satisfy (17a) (17b) for . Together, (16a)–(17b) provide an elegant and accurate approach to describe the mean bacterial loads for the first 48 hours in the lungs, as well as in other organs.

Parameter inference

Total-order Sobol sensitivity indices, presented in Fig 8, quantify the overall effect of a single parameter on model output [61] with respect to total bacterial counts in the lung and MLNs. The parameters were varied over the ranges indicated by the prior distributions in Table 1, with ϕ ∈ [0.5, 5] and log10 μE ∈ [−4, −1]. The parameter β is initially the most important, with δ having an increasing effect at later times. If infected macrophages rupture quickly, bacteria are not able to replicate as effectively in the cytosol and are more frequently found in extracellular environments and in phagosomes, where bacterial replication does not take place. Thus, larger values of δ are associated with slower bacterial population growth. In addition to β and δ, the per bacterium phagocytosis rate, Mρ, and the total migration rate, γ, are important for describing the dynamics outside the lungs. Together with μE, they determine with what probability an extracellular bacterium in the lung is killed, migrates to a different organ, or infects a cell in the lung. For the two remaining parameters, μE and ϕ, the total-order Sobol indices remain low during the initial 48 hours. This allows us to fix their values when performing the Bayesian inference, with the confidence that any uncertainty in these estimates will have little effect on the total bacterial counts. We therefore only expect to learn about β, δ, Mρ and γ, and thus, set μE = 0.01 h−1 and ϕ = 2 h−1 [62]. In this section, the decision has also been made to fix the rate of intracellular bacterial death to zero; that is, μ = 0, given the belief that macrophages will likely rupture rather than clear their bacterial load. Including μ in the inference would also affect our ability to learn about β, with these two parameters not identifiable from measurements of bacterial counts alone. Finally, the values of the weights that dictate which organs bacteria exiting the lung migrate to are selected, based on the data (summarised in Table 2). displayed in Fig 9. In order to do so we consider the fraction of bacteria that are present in each organ after 48 hours. The mean fraction of bacteria in each organ was used to obtain the following weights: wMLN = 0.8, wliver = 0.11, wspleen = 0.05 and wkidney = 0.04. The larger weight assigned to the MLNsis reasonable, since bacteria may be drained rapidly through the lymphatic system to the MLN [63, 64].

thumbnail
Download:
Fig 8. Total-order Sobol indices.

An indication of the most important parameters that describe the dynamics of total bacterial counts in the lung (left) and MLNs(right) during the first 48 hours of infection. The ranges over which each parameter is varied are: (Mρ) ∈ [10−2 h−1, 105 h−1], ϕ ∈ [0.5h−1, 5h−1], β ∈ [10−2 h−1, 1h−1], δ ∈ [10−5 h−1, 10−1 h−1], μE ∈ [10−4 h−1, 10−1 h−1] and γ ∈ [1h−1, 103 h−1].

https://doi.org/10.1371/journal.pcbi.1007752.g008

thumbnail
Download:
Fig 9. Bacterial counts in the lung, MLN, liver, kidney and spleen.

Mice are exposed to either 160.33 CFUs (high), 13.7 CFUs (medium) or 2 CFUs (low) of Francisella tularensis SCHU S4 bacteria. The observed data are denoted by shaded points, whilst the geometric mean and standard deviation are represented by solid points and bars, respectively. Zero counts have been replaced by one in order to calculate the geometric mean and standard deviation.

https://doi.org/10.1371/journal.pcbi.1007752.g009

thumbnail
Download:
Table 1. Agent-based model parameters.

A description and value of the parameters and prior distributions used to determine the total bacterial load in each organ with approximate Bayesian inference. Migration probabilities are calculated using the proportion of bacteria in each organ after 48 hours, starred parameters are inferred using ABC from the observed bacterial counts in each organ.

https://doi.org/10.1371/journal.pcbi.1007752.t001

thumbnail
Download:
Table 2. Table with experimental data sets. Bacterial counts in the lung, MLN, liver, kidney and spleen used for the parameter inference.

Mice are exposed to either 160.33 CFU (top) or 13.7 CFU (bottom) of Francisella tularensis SCHU S4 bacteria. Geometric means and standard deviations (SD) are also given.

https://doi.org/10.1371/journal.pcbi.1007752.t002

The aim of this section is to make use of the experimental data (see Fig 9) and the mathematical cohort model described in the previous section, to learn about the selected model parameters with an approximate Bayesian computation (ABC) rejection sampling algorithm [42]. Initial uncertainty regarding each parameter is encoded in the prior distributions described in Table 1. For every set of parameters sampled from the prior distributions, the total number of bacteria in each organ is calculated (mod), and compared to the experimental observations (exp) using the distance function where is the set of initial doses and is the set of times at which measurements are provided for a given dose and organ . For the lung, model predictions of bacterial burdens, , are found by summing the initial three cohorts of phagosomal and cytosolic bacteria, (12)–(16b). For each of the remaining organs, where it is assumed that a bacterium is only able to infect one macrophage during the initial 48 hours, is found by summing (17a) and (17b). For the geometric mean, , and geometric standard deviation, , experimental observations yielding no bacteria are set to one bacterium.

In total, 106 iterations of the rejection sampling algorithm were performed. The acceptance rate is 0.5%, which leads to an accepted posterior sample of 5×103 parameter sets. Pointwise median predictions, provided in Fig 10, confirm that the posterior samples can reproduce the behaviour experimentally observed during the first 48 hours of infection.

thumbnail
Download:
Fig 10. Pointwise median predictions.

A comparison between model predictions of total bacterial counts and observed bacteria counts for medium (left) and high (right) initial doses. Solid curves and shaded regions, respectively, denote pointwise median predictions and 95% credible regions. These have been constructed using all parameter sets from the posterior sample.

https://doi.org/10.1371/journal.pcbi.1007752.g010

The posterior distribution of β, shown in Fig 11, is narrow. This indicates that the experimental data together with the mathematical model allows us to learn about this parameter. We find β = 0.154 ± 0.014 h−1, a range that includes the value considered in Ref. [23] based on reports that the doubling time of a single Francisella tularensis bacterium is approximately five hours. The wide posterior distribution of δ, in Fig 11, suggests that it is only possible to identify an upper bound of δ = 10−2 h−1. However, if τrupture is a random variable for the time a Yule-catastrophe process [34] takes to reach the rupture state H, the mean time until rupture is where f(t) is the density of the time until rupture, given in (7). As this is a function of β and δ only, and β is confined to a small range of values (as described above), this independent estimate of the mean rupture time allows us to improve our learning about δ.

thumbnail
Download:
Fig 11. Posterior histograms for β and δ.

From the posterior sample, the histogram for β (left) shows the significant learning that has been achieved making use of the experimental data. When comparing to bacterial counts, only an upper bound for δ can be identified (centre). However, this distribution can be refined when also considering model predictions for the mean time until rupture (right).

https://doi.org/10.1371/journal.pcbi.1007752.g011

Wood et al., by comparing to data from an in vitro study involving the infection of macrophages with Francisella tularensis bacteria, estimate a mean rupture of time of 44.4 h [30]. By choosing pairs (β, δ) such that the mean rupture time is 44.4±0.5 h, this additional knowledge allows us to refine the posterior distribution of δ (see Fig 11 (right)), now yielding a posterior median estimate of δmedian ≈ 1.5 × 10−4 h−1 and a significantly narrower range of 6.3 × 10−5 h−1 ≤ δ ≤ 3.8 × 10−4 h−1.

For the parameters Mρ and γ, refined individual learning is not possible. However, the bivariate posterior histogram in Fig 12 shows that a strong correlation exists between these parameters: increasing Mρ increases the likelihood that an extracellular bacterium again infects a macrophage in the lung, which can be balanced by also increasing the rate at which bacteria leave the lung. Summary statistics for each of the posterior samples are reported in Table 3. A useful quantity is , which is the probability that a bacterium migrates to a different organ, rather than dying or infecting another macrophage in the lung. The corresponding posterior distribution, constructed using the posterior samples of Mρ and γ, along with μE = 0.01 h−1, is provided in Fig 12. Here, a posterior median value suggests that approximately 4% of extracellular bacteria in the lung are directly involved in the early dissemination (first 48 hours) of Francisella tularensis infection to other organs.

thumbnail
Download:
Fig 12. Posterior histogram for Mρ and γ.

Bivariate histogram (left) depicting the strong correlation between Mρ and γ and a histogram of the migration probability (right), constructed using posterior samples of Mρ, γ and μE = 0.01 h−1.

https://doi.org/10.1371/journal.pcbi.1007752.g012

thumbnail
Download:
Table 3. Summary statistics of the posterior sample for each parameter included in the approximate Bayesian inference.

Posterior samples contain 5 × 103 values.

https://doi.org/10.1371/journal.pcbi.1007752.t003

It is often difficult to predict the course of infection when infecting mice at low initial doses of bacteria. Here, of the 24 mice infected with 2 CFUs of Francisella tularensis bacteria and culled between 12 and 48 hours post-infection, only seven had detectable levels of bacteria present in their lungs. Only one mouse had detectable levels in any of the remaining organs measured. These bacterial counts are presented in Fig 13, alongside model predictions obtained using the posterior distributions that were previously inferred from the medium and high infectious dose data. The predictions for the lung agree well with the observed bacterial counts, whilst the predictions for other organs are informative for understanding expected disease progression. More importantly at low initial doses, the stochastic nature of the agent based model described here would allow us to estimate the probability that all bacteria are cleared and the mouse recovers.

thumbnail
Download:
Fig 13. Low infectious dose predictions.

Predictions of total bacterial counts in each organ following infection at a low initial dose. Posterior distributions inferred by performing ABC with the medium and high doses have been used to create pointwise median predictions and 95% credible regions.

https://doi.org/10.1371/journal.pcbi.1007752.g013

Discussion

There is a tradition of mathematical models that consider populations of infectious agents, target cells and infected cells [24, 25, 66–70]. The usual assumption that new infectious agents are produced at a rate proportional to the number of infected cells, perhaps after an “eclipse” phase [27, 28], may be appropriate in situations where infected cells, independently, release new infectious agents, one or a few at a time, on multiple occasions during their lifetime, a process known as “budding” [32]. It is more problematic when infectious particles are released in a single “burst” as the infected cell dies. The burst scenario is found in many types of infection, including the pathogen of interest in this work, Francisella tularensis.

Our approach to understanding the first 48 hours of infection is to focus on dynamics within a single cell. Our within-macrophage model describes the intracellular replication of bacteria and the rupturing and death of the macrophage. Because a rupture event immediately releases all the bacteria from a cell, the stochastic process is a birth-death-catastrophe process (or Yule-catastrophe process if the death rate is zero) [34]. We assume that the rupture rate of a macrophage containing Xt cytosolic bacteria at time t is δXt. In the first 48 hours post-infection, the bacteria released in catastrophe events rapidly enter new macrophages. Therefore, the growth rate of the total bacterial load is not strongly dependent on this assumption or on the value of δ. We show that Xt has a geometric distribution, the probability density of macrophage rupture times is and the mean rate of release of bacteria, as a function of time, is .

Our modelling approach is applicable to other intracellular pathogens, such as Salmonella enterica, and Bacillus anthracis, where models must also consider germination of spores [71]. Levofloxacin and ciprofloxacin are antibiotics commonly used to treat tularemia, although their success relies on early administration, which is often made difficult given the non-specific symptoms [72]. Pharmacokinetic and pharmacodynamic models can be used to describe the concentration of antibiotic in each organ and the effect it has on the bacterial load [73]. In epidemic models, birth and death events describe the infection and recovery of animals within a disease reservoir; an analogue of a catastrophe event is the “spillover” of disease into a human population [74]. In models of receptor-mediated signalling events [75], a sequence of states is used to represent reversible phosphorylation events initiated by the binding of a ligand, with dissociation of the ligand leading to the termination of the signal.

In addition to the mathematical analysis focusing on a single cell, we have also described Francisella tularensis pathogenesis with an agent-based computational model. Because the history and family tree of every bacterium is available in an agent-based model, it is useful to classify bacteria by cohorts, according to how many different macrophages the bacterium and its ancestors have entered. Knowing the distributions of rupture time and numbers of bacteria released allows us to provide a cohort description of the total bacterial load, in phagosomes and cytosols of infected macrophages, in the lung and in other organs.

With experimental data (see Fig 9) and approximate Bayesian computation, we learn about the parameters of the cohort model. The rate of growth of the experimentally measured bacterial load is primarily determined by the rate of intracellular bacterial division, a parameter that does not appear in a model that only counts numbers of free bacteria and infected cells. Our experimental data, from in vivo measurements in lungs, lymph nodes, spleen, kidney and liver, have allowed us to tightly confine the posterior distribution of β to a range that is consistent with published estimates based on in vitro data. On the other hand, the experimental data do not allow us to determine the mode of migration from the lung to other organs, nor to place tight constraints on the associated timescales.

Data from human or primate infection is even more rare than murine data [76, 77]. In silico models serve as a bridge between animal and human research, with the advantage that human pharmacokinetic and pharmacodynamic parameters can be directly applied. Mathematical models, suitably developed and validated, can provide a suite of tools to estimate the result of experiments, inform their design and extrapolate to humans.

Materials and methods

Experimental procedures

Six-to-eight week old female BALB/c mice were challenged (inhalational exposure) with Francisella tularensis SCHU S4. In these experiments, mice were infected with either 160 (high), 14 (medium) or 2 (low) colony forming units. At each challenge dose, six mice were culled at 1, 18, 24, 48, 72 and 96 hours and the bacterial burden was measured in the lung, liver, spleen, kidney, MLNsand blood. An additional measurement was taken at 12 hours from mice receiving the low or medium dose. All manipulations were carried out under Advisory Committee for Dangerous Pathogens Level 3 containment conditions in a (Level 3 containment) safety cabinet complying with BS 5726. Francisella tularensis SCHU S4 was cultured from frozen stock for two days on blood cysteine glucose agar (BCGA) with cysteine at 37°C. Subsequently, bacteria were harvested to inoculate 50 ml of modified cysteine partial hydrolysate broth with cysteine and glucose and incubated overnight at 37°C on a rotary shaker (150 rpm). The suspension was then adjusted using phosphate-buffered saline (PBS) until the optical density at 590 nm was 0.10, where the estimated bacterial density would be 5 × 108 CFU per ml. Bacterial numbers for challenge were determined on agar following serial dilution (1:10) of samples. The work was conducted under the terms of a licence granted in accordance with the UK Animal (Scientific Procedures) Act, 1986. Female BALB/c mice (Charles River Laboratories Ltd, Margate, Kent, UK) were habituated to the experimental animal unit for one week prior to challenge. Environmental conditions were maintained at 21°C ± 2°C and 55% ± 10% humidity with lighting set to mimic a 12/12 (hour) dawn to dusk cycle. The mice were then transferred to a Level 3 containment rigid isolator for a further 5 and 7 days. They were housed in polycarbonate cages (six animals per cage) with steel cage tops and corncob bedding (International Product Supplies, Wellingborough, UK). The mice were fed a Teklad TRM 19% protein irradiated diet ad libitum (Harlan Teklad, Bicester, UK) and given fresh water daily.

Mice were challenged by aerosol, using a Henderson-type apparatus [78] and a Collison nebuliser [79]. Briefly, 10 ml of Francisella tularensis SCHU S4 culture was aerosolised using a Henderson Apparatus over an exposure time of 10 minutes [80]. The aerosol was delivered at a flow rate of 12 L/min with impinger samples from the exposure apparatus plated on BCGA to calculate retained dose. Using the known flow rate of the Henderson exposure apparatus (66 L/min), bacterial counts from these samples were then converted to bacterial counts per litre of air. The breathing rate of the animals in the apparatus (approximately 20 ml of air per minute) was then added to the calculation along with the length of exposure (10 minutes) to yield an estimated delivered dose expressed in CFU per animal. Previous studies have determined that aerosol uptake in obligate nasal breathers, such as the mouse, is approximately 40% [81]. Using this conversion factor, the estimated retained dose was calculated for each exposure. Mice were culled for analysis of tissues at different time points and exsanguinated using cardiac puncture following terminal anaesthesia. Blood was placed into 1.5 ml heparin tubes. Lungs, spleen and liver were removed and placed into bijoux tubes filled with 2 ml of PBS. Duplicate experiments were performed. All procedures and housing were in accordance with the Animal (Scientific Procedures) Act (1986). Organs were processed at less than 1 h post-mortem. Blood was diluted 1:10 in PBS. Collected organs were placed into 6-well trays containing 40 μm cell sieves with 1,800 μl of PBS, then disrupted through the cell sieve using the plunger of a 2 ml syringe. Cell suspensions were collected. 100 μl aliquots of the cell suspension or blood were used for enumeration of bacteria on agar plates following serial dilution in PBS.

Agent-based model

In our computational model, each macrophage and each free bacterium has a unique identity and set of mutable attributes. The attributes of a macrophage are: spatial location, state of activation, cohort counter (see Cohort analysis Section), and list (and number) of phagosomal and cytosolic bacteria. The spatial location is either lung, liver, spleen, MLNsor kidney. Although bacterial counts have also been measured in the blood, the numbers were small enough (<10 CFUs) that this compartment can been neglected in the model, as a first approximation. Once phagocytosed, a bacterium remains in a macrophage’s phagosome for an exponentially-distributed time with mean 1/ϕ, then escapes to the macrophage’s cytosol. There, it becomes the founder of a population of intracellular bacteria, governed by the birth-death-catastrophe process, that lasts until either the bacterial population is eliminated from the macrophage, or the macrophage ruptures and releases its contents. Newly-released bacteria suffer one of three fates: phagocytosis by a macrophage in the same organ, death or migration to a different organ. Given that the initial number of resting macrophages is much larger than the initial number of bacteria, events in which an infected macrophage is reinfected by another bacterium are rare in the first 72 hours of pathogenesis.

Macrophages may exhibit a variety of activation states in different tissues [71, 82–85]. At any time in our computational model, each macrophage in a computational realisation is in one of three states: resting, suppressed (anti-inflammatory) or activated (pro-inflammatory). Every one of the initial M macrophages begins in the resting state. On phagocytosis, resting macrophages enter a suppressed state in which they are unable to kill bacteria and secrete the anti-inflammatory cytokine TGF-β that contributes to the suppression of other macrophages. Resting macrophages can become activated through the detection of host damage molecules released by rupturing macrophages [22]: each rupture event can result in the activation of a resting macrophage in the same organ. Activated macrophages kill the bacteria they phagocytose. They also secrete IFN-γ that provokes activation of neighbouring macrophages.

A numerical realisation starts with the arrival of a number, chosen from a Poisson distribution with mean N, of free bacteria in the lung (see Table 1 for values of N). Individual rupture times are recorded, and we explicitly track cohorts of bacteria by assigning a “cohort number” attribute to each bacterium (see Cohort analysis Section). Similarly, the set of bacteria inside each macrophage is subdivided by cohort number and each rupture event is classified according to the minimum cohort number of the bacteria released. Computer codes (in Python) to generate the numerical realisations of the agent-based model and to perform the cohort analysis are available in this link http://review.researchdata.leeds.ac.uk/id/eprint/1399/.

There are nine types of events in the computational agent-based model: phagocytosis, escape, division, intracellular or extracellular death of a bacterium, rupture, migration, cytokine-mediated suppression or activation of a macrophage. The rates are as follows:

  • ρ is the rate of phagocytosis per macrophage,
  • ϕ is the rate at which bacteria escape the phagosome,
  • β and μ are the birth and death rates of bacteria in the cytosol,
  • μE is the death rate of free (or extracellular) bacteria,
  • δXt is the rupture rate of a macrophage containing Xt cytosolic bacteria at time t, and
  • γ is the rate at which bacteria migrate to other organs.

Cytokine-mediated activation and suppression of macrophages are included in the computational model by means of two dimensionless functions of time, G(t) and T(t), in each organ. The first summarises the levels of inflammatory cytokines, such as IFNγ; the second summarises the levels of anti-inflammatory cytokines, such as TGF-β. The functions are updated according to the following differential equations where MA(t) and MS(t) are the numbers of activated and suppressed macrophages at time t, respectively. The parameters have been chosen as follows: αG = 10−3 h−1 and αT = 10−1 h−1 are the per macrophage production rates of TGF-β and IFNγ, respectively; the decay rates are μG = 9 × 10−2 h−1 and μT = 10−1 h−1, respectively [86]. At any time when T(t) exceeds the threshold level 102 in any organ, each resting macrophage in that organ has a rate νγ = 0.04h−1 of transition to the suppressed state. At any time when G(t) exceeds the threshold level 102 in any organ, each resting macrophage in that organ has a rate νβ = 0.01h−1 of transition to the activated state. Although cytokine-mediated events have been included in the computational agent-based model, their effect in Francisella tularensis infected mice is minimal during the initial 72 hours [87], and thus, they are not included in the subsequent approximations or parameter inference. Despite this, future experimental measurements of the concentration of these cytokines would enable us to use this agent-based model, along with any learning about the remaining model parameters achieved here, to obtain more accurate estimates of parameters, such as the level of IFN-γ required for macrophage activation.

Two types of time-stepping are available:

  • The Gillespie algorithm, where the time increments are inter-event times which are drawn from exponential distributions; the probability of each type of event is proportional to its rate, and all rates are updated after each event [23].
  • Tau-leaping, where the time increment (or step size), Δt, is fixed; the number of occurrences of each type of event per step is a Poisson random variable with mean proportional to its rate [88].

When results from agent-based simulations are reported here, the tau-leaping procedure has been applied with a step size (or time increment) of Δt = 10−2 hours. Due to the large number of macrophages present at the start of the simulation, along with the rapid growth of the bacterial population, the number of agents is too large for the Gillespie algorithm to be implemented efficiently.

Solution of the probability generating function of a birth-and-death process with catastrophe

Let G(z, t) denote the probability generating function of a birth-and-death process with catastrophe, We then have Using the method of characteristics, we may write Given the initial condition G(z, 0) = z and the substitution ξ = (z − a)/(z − b), we find Setting ξ = e−β(b−a)t(z − a)/(z − b), one obtains If the process instead starts with Xt = k, then G(k)(z, 0) = zk and the solution is G(k)(z, t) = [G(z, t)]k.

Acknowledgments

We acknowledge and are grateful to a number of technical staff and scientists who assisted in performing these experimental procedures. Numerical work was undertaken on ARC3, which is part of the High Performance Computing facilities at the University of Leeds, UK. We acknowledge and are grateful to the International Centre for Mathematical Sciences (ICMS), Edinburgh, where this manuscript was completed during a Research in Groups programme (JC, GL, MLG, TRL and CMP). JC, MLG and CMP acknowledge the support and hospitality of ICTS, Bangalore, India, where the final version of this manuscript was discussed and completed. Content includes material subject to Dstl Crown copyright (2019).

References

  1. 1. Kortepeter MG, Parker GW. Potential biological weapons threats. Emerging infectious diseases. 1999;5(4):523. pmid:10458957
  2. 2. Dennis DT, Inglesby TV, Henderson DA, Bartlett JG, Ascher MS, Eitzen E, et al. Tularemia as a biological weapon: medical and public health management. Jama. 2001;285(21):2763–2773. pmid:11386933
  3. 3. Ellis J, Oyston P, Green M, Titball R. Tularemia. Tularemia Clinical microbiology reviews. 2002;15(4). pmid:12364373
  4. 4. Fortier AH, Slayter M, Ziemba R, Meltzer M, Nacy C. Live vaccine strain of Francisella tularensis: infection and immunity in mice. Infection and immunity. 1991;59(9):2922–2928. pmid:1879918
  5. 5. Oyston PC. Francisella tularensis vaccines. Vaccine. 2009;27:D48–D51. pmid:19837286
  6. 6. Piercy T, Steward J, Lever M, Brooks T. In vivo efficacy of fluoroquinolones against systemic tularaemia infection in mice. Journal of Antimicrobial Chemotherapy. 2005;56(6):1069–1073. pmid:16223941
  7. 7. Oyston PC, Sjöstedt A, Titball RW. Tularaemia: bioterrorism defence renews interest in Francisella tularensis. Nature Reviews Microbiology. 2004;2(12):967–978. pmid:15550942
  8. 8. McCaffrey RL, Allen LAH. Francisella tularensis LVS evades killing by human neutrophils via inhibition of the respiratory burst and phagosome escape. Journal of Leukocyte Biology. 2006;80(6):1224–1230. pmid:16908516
  9. 9. Oyston PC. Francisella tularensis: unravelling the secrets of an intracellular pathogen. Journal of Medical Microbiology. 2008;57(8):921–930. pmid:18628490
  10. 10. Larsson P, Oyston PC, Chain P, Chu MC, Duffield M, Fuxelius HH, et al. The complete genome sequence of Francisella tularensis, the causative agent of tularemia. Nature Genetics. 2005;37(2):153–159. pmid:15640799
  11. 11. Mosser DM. The many faces of macrophage activation. Journal of Leukocyte Biology. 2003;73(2):209–212. pmid:12554797
  12. 12. Clemens DL, Lee BY, Horwitz MA. Francisella tularensis enters macrophages via a novel process involving pseudopod loops. Infection and Immunity. 2005;73(9):5892–5902. pmid:16113308
  13. 13. Bosio C. The subversion of the immune system by Francisella tularensis. Frontiers in Microbiology. 2011;2. pmid:21687406
  14. 14. Jones CL, Napier BA, Sampson TR, Llewellyn AC, Schroeder MR, Weiss DS. Subversion of host recognition and defense systems by Francisella spp. Microbiology and Molecular Biology Reviews. 2012;76(2):383–404. pmid:22688817
  15. 15. Dai S, Rajaram MV, Curry HM, Leander R, Schlesinger LS. Fine tuning inflammation at the front door: macrophage complement receptor 3-mediates phagocytosis and immune suppression for Francisella tularensis. PLoS Pathogens. 2013;9(1):e1003114. pmid:23359218
  16. 16. del Barrio L, Sahoo M, Lantier L, Reynolds JM, Ceballos-Olvera I, Re F. Production of anti-LPS IgM by B1a B cells depends on IL-1β and is protective against lung infection with Francisella tularensis LVS. PLoS pathogens. 2015;11(3):e1004706. pmid:25768794
  17. 17. Okan NA, Kasper DL. The atypical lipopolysaccharide of Francisella. Carbohydrate research. 2013;378:79–83. pmid:23916469
  18. 18. Golovliov I, Baranov V, Krocova Z, Kovarova H, Sjöstedt A. An attenuated strain of the facultative intracellular bacterium Francisella tularensis can escape the phagosome of monocytic cells. Infection and immunity. 2003;71(10):5940–5950. pmid:14500514
  19. 19. Jones JW, Broz P, Monack DM. Innate immune recognition of Francisella tularensis: activation of type-I interferons and the inflammasome. Frontiers in Microbiology. 2011;2. pmid:21687410
  20. 20. Leander R, Dai S, Schlesinger LS, Friedman A. A Mathematical Model of CR3/TLR2 Crosstalk in the Context of Francisella tularensis Infection. PLoS Computational Biology. 2012;8(11):e1002757. pmid:23133361
  21. 21. Celli J, Zahrt TC. Mechanisms of Francisella tularensis intracellular pathogenesis. Cold Spring Harbor Perspectives in Medicine. 2013;3(4):a010314. pmid:23545572
  22. 22. D’Elia RV, Harrison K, Oyston PC, Lukaszewski RA, Clark GC. Targeting the “cytokine storm” for therapeutic benefit. Clin Vaccine Immunol. 2013;20(3):319–327. pmid:23283640
  23. 23. Gillard JJ, Laws TR, Lythe G, Molina-París C. Modeling early events in Francisella tularensis pathogenesis. Frontiers in cellular and infection microbiology. 2014;4:169. pmid:25566509
  24. 24. Perelson AS, Nelson PW. Mathematical analysis of HIV-I: dynamics in vivo. SIAM Review. 1999;41(1):3–44.
  25. 25. Ciupe SM, Heffernan JM. In-host modeling. Infectious Disease Modelling. 2017;2(2):188–202. pmid:29928736
  26. 26. Conway JM, Perelson AS. Early HIV Infection Predictions: Role of Viral Replication Errors. SIAM Journal on Applied Mathematics. 2018;78(4):1863–1890. pmid:31231142
  27. 27. Handel A, Liao LE, Beauchemin CA. Progress and trends in mathematical modelling of influenza A virus infections. Current Opinion in Systems Biology. 2018;.
  28. 28. Conway JM, Ribeiro RM. Modeling the immune response to HIV infection. Current Opinion in Systems Biology. 2018;.
  29. 29. Koelle K, Farrell AP, Brooke CB, Ke R. Within-host infectious disease models accommodating cellular coinfection, with an application to influenza. Virus evolution. 2019;5(2):vez018.
  30. 30. Wood R, Egan J, Hall I. A dose and time response Markov model for the in-host dynamics of infection with intracellular bacteria following inhalation: with application to Francisella tularensis. Journal of The Royal Society Interface. 2014;11(95):20140119.
  31. 31. Pantha B, Cross A, Lenhart S, Day J. Modeling the macrophage-anthrax spore interaction: Implications for early host-pathogen interactions. Mathematical biosciences. 2018;305:18–28.
  32. 32. Pearson JE, Krapivsky P, Perelson AS. Stochastic theory of early viral infection: continuous versus burst production of virions. PLoS Computational Biology. 2011;7(2):e1001058.
  33. 33. Carruthers J, López-García M, Gillard JJ, Laws TR, Lythe G, Molina-París C. A novel stochastic multi-scale model of Francisella tularensis infection to predict risk of infection in a laboratory. Frontiers in Microbiology. 2018;9.
  34. 34. Karlin S, Tavaré S. Linear birth and death processes with killing. Journal of Applied Probability. 1982;19(3):477–487.
  35. 35. Brock SR, Parmely MJ. Complement C3 as a prompt for human macrophage death during infection with Francisella tularensis strain SCHU S4. Infection and Immunity. 2017;85(10):e00424–17.
  36. 36. Lai XH, Golovliov I, Sjöstedt A. Francisella tularensis induces cytopathogenicity and apoptosis in murine macrophages via a mechanism that requires intracellular bacterial multiplication. Infection and Immunity. 2001;69(7):4691–4694.
  37. 37. Brown SP, Cornell SJ, Sheppard M, Grant AJ, Maskell DJ, Grenfell BT, et al. Intracellular demography and the dynamics of Salmonella enterica infections. PLoS biology. 2006;4(11):e349.
  38. 38. Mastroeni P, Grant A, Restif O, Maskell D. A dynamic view of the spread and intracellular distribution of Salmonella enterica. Nature Reviews Microbiology. 2009;7(1):73.
  39. 39. Gog JR, Murcia A, Osterman N, Restif O, McKinley TJ, Sheppard M, et al. Dynamics of Salmonella infection of macrophages at the single cell level. Journal of The Royal Society Interface. 2012;9(75):2696–2707.
  40. 40. Sobol IM. Sensitivity estimates for nonlinear mathematical models. Mathematical modelling and computational experiments. 1993;1(4):407–414.
  41. 41. Saltelli A, Tarantola S, Campolongo F, Ratto M. Sensitivity analysis in practice: a guide to assessing scientific models. Chichester, England. 2004;.
    • 42. Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MP. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface. 2008;6(31):187–202.
    • 43. Karlin S, McGregor J. The classification of birth and death processes. Transactions of the American Mathematical Society. 1957;86(2):366–400.
    • 44. Novozhilov AS, Karev GP, Koonin EV. Biological applications of the theory of birth-and-death processes. Briefings in Bioinformatics. 2006;7(1):70–85.
    • 45. Brockwell PJ. The extinction time of a birth, death and catastrophe process and of a related diffusion model. Advances in Applied Probability. 1985;17(1):42–52.
    • 46. Brockwell PJ. The extinction time of a general birth and death process with catastrophes. Journal of Applied Probability. 1986;23(4):851–858.
    • 47. Brockwell PJ, Gani J, Resnick SI. Birth, immigration and catastrophe processes. Advances in Applied Probability. 1982;14(4):709–731.
    • 48. Cairns B, Pollett P. Extinction times for a general birth, death and catastrophe process. Journal of applied probability. 2004;41(4):1211–1218.
    • 49. Pakes AG. Limit theorems for the population size of a birth and death process allowing catastrophes. Journal of Mathematical Biology. 1987;25(3):307–325.
    • 50. Pakes AG, Pollett P. The supercritical birth, death and catastrophe process: limit theorems on the set of extinction. Stochastic Processes and their Applications. 1989;32(1):161–170.
    • 51. Pollett P, Zhang H, Cairns BJ. A note on extinction times for the general birth, death and catastrophe process. Journal of applied probability. 2007;44(2):566–569.
    • 52. Di Crescenzo A, Giorno V, Nobile AG, Ricciardi LM. A note on birth–death processes with catastrophes. Statistics & Probability Letters. 2008;78(14):2248–2257.
    • 53. Krinik A, Mortensen C. Transient probability functions of finite birth–death processes with catastrophes. Journal of statistical planning and inference. 2007;137(5):1530–1543.
    • 54. Zeifman A, Satin Y, Panfilova T. Limiting characteristics for finite birth–death-catastrophe processes. Mathematical Biosciences. 2013;245(1):96–102.
    • 55. van Doorn EA. Conditions for the existence of quasi-stationary distributions for birth–death processes with killing. Stochastic processes and their applications. 2012;122(6):2400–2410.
    • 56. Van Doorn EA, Zeifman AI. Extinction probability in a birth-death process with killing. Journal of applied probability. 2005;42(1):185–198.
    • 57. Van Doorn EA, Zeifman AI. Birth-death processes with killing. Statistics and Probability Letters. 2005;72(1):33–42.
    • 58. Lai XH, Golovliov I, Sjöstedt A. Expression of IglC is necessary for intracellular growth and induction of apoptosis in murine macrophages by Francisella tularensis. Microbial Pathogenesis. 2004;37(5):225–230.
    • 59. Attie O, Daefler S. An agent-based model of tularemia. J Data Mining Genomics Proteomics. 2013;4(1):2153–0602.
    • 60. Gallagher M, Brooke C, Ke R, Koelle K. Causes and Consequences of Spatial Within-Host Viral Spread. Viruses. 2018;10(11):627.
    • 61. Zhang XY, Trame M, Lesko L, Schmidt S. Sobol sensitivity analysis: a tool to guide the development and evaluation of systems pharmacology models. CPT: pharmacometrics & systems pharmacology. 2015;4(2):69–79.
    • 62. Ray K, Marteyn B, Sansonetti PJ, Tang CM. Life on the inside: the intracellular lifestyle of cytosolic bacteria. Nature Reviews Microbiology. 2009;7(5):333.
    • 63. Bosio CM, Bielefeldt-Ohmann H, Belisle JT. Active suppression of the pulmonary immune response by Francisella tularensis Schu4. The Journal of Immunology. 2007;178(7):4538–4547.
    • 64. Hall CA, Flick-Smith HC, Harding SV, Atkins HS, Titball RW. A Bioluminescent Francisella tularensis SCHU S4 Strain Enables Noninvasive Tracking of Bacterial Dissemination and the Evaluation of Antibiotics in an Inhalational Mouse Model of Tularemia. Antimicrobial agents and chemotherapy. 2016;60(12):7206–7215.
    • 65. Lowrie D, Aber V, Carrol M. Division and death rates of Salmonella typhimurium inside macrophages: use of penicillin as a probe. Microbiology. 1979;110(2):409–419.
    • 66. Perelson AS, Kirschner DE, De Boer R. Dynamics of HIV infection of CD4+ T cells. Mathematical Biosciences. 1993;114(1):81–125.
    • 67. Nowak M, May RM. Virus dynamics: mathematical principles of immunology and virology: mathematical principles of immunology and virology. Oxford University Press, UK; 2000.
      • 68. Wodarz D, Nowak MA. Mathematical models of HIV pathogenesis and treatment. BioEssays. 2002;24(12):1178–1187.
      • 69. De Boer RJ. Which of our modeling predictions are robust? PLoS computational biology. 2012;8(7):e1002593.
      • 70. Perelson AS, Ribeiro RM. Introduction to modeling viral infections and immunity. Immunological Reviews. 2018;285(1):5–8.
      • 71. Day J, Friedman A, Schlesinger LS. Modeling the immune rheostat of macrophages in the lung in response to infection. PNAS. 2009;106(27):11246–11251.
      • 72. Caspar Y, Maurin M. Francisella tularensis susceptibility to antibiotics: a comprehensive review of the data obtained in vitro and in animal models. Frontiers in cellular and infection microbiology. 2017;7:122.
      • 73. Kirschner D, Pienaar E, Marino S, Linderman JJ. A review of computational and mathematical modeling contributions to our understanding of Mycobacterium tuberculosis within-host infection and treatment. Current opinion in systems biology. 2017;3:170–185.
      • 74. Singh S, Schneider DJ, Myers CR. Using multitype branching processes to quantify statistics of disease outbreaks in zoonotic epidemics. Physical Review E. 2014;89(3):032702.
      • 75. Castro M, López-García M, Lythe G, Molina-París C. First passage events in biological systems with non-exponential inter-event times. Scientific reports. 2018;8(1):15054.
      • 76. Nelson M, Lever MS, Dean RE, Savage VL, Salguero FJ, Pearce PC, et al. Characterization of lethal inhalational infection with Francisella tularensis in the common marmoset (Callithrix jacchus). Journal of Medical Microbiology. 2010;59(Pt 9):1107.
      • 77. Hamblin KA, Armstrong SJ, Barnes KB, Davies C, Laws TR, Blanchard JD, et al. Inhaled liposomal ciprofloxacin protects against a lethal infection in a murine model of pneumonic plague. Frontiers in Microbiology. 2017;8:91.
      • 78. Druett H. A mobile form of the Henderson apparatus. Epidemiology and Infection. 1969;67(3):437–448.
      • 79. May K, Harper G. The efficiency of various liquid impinger samplers in bacterial aerosols. British Journal of Industrial Medicine. 1957;14(4):287.
      • 80. Eyles J, Hartley M, Laws T, Oyston P, Griffin K, Titball R. Protection afforded against aerosol challenge by systemic immunisation with inactivated Francisella tularensis live vaccine strain (LVS). Microbial pathogenesis. 2008;44(2):164–168.
      • 81. Harper G, Morton J. A method for measuring the retained dose in experiments on airborne infection. Epidemiology & Infection. 1962;60(2):249–257.
      • 82. Gordon S, Taylor PR. Monocyte and macrophage heterogeneity. Nature Reviews Immunology. 2005;5(12):953–964.
      • 83. Mosser DM, Edwards JP. Exploring the full spectrum of macrophage activation. Nature Reviews Immunology. 2008;8(12):958–969.
      • 84. Murray PJ, Wynn TA. Protective and pathogenic functions of macrophage subsets. Nature Reviews Immunology. 2011;11(11):723.
      • 85. Gordon S, Plüddemann A. Tissue macrophages: heterogeneity and functions. BMC biology. 2017;15(1):53.
      • 86. Marino S, Kirschner DE. The human immune response to Mycobacterium tuberculosis in lung and lymph node. Journal of Theoretical Biology. 2004;227(4):463–486.
      • 87. D’Elia RV, Laws TR, Carter A, Lukaszewski R, Clark GC. Targeting the “Rising DAMP” during a Francisella tularensis infection. Antimicrobial agents and chemotherapy. 2013;57(9):4222–4228.
      • 88. Gillespie DT. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007;58:35–55.