Intracellular Demography and the Dynamics of Salmonella enterica Infections

An understanding of within-host dynamics of pathogen interactions with eukaryotic cells can shape the development of effective preventive measures and drug regimes. Such investigations have been hampered by the difficulty of identifying and observing directly, within live tissues, the multiple key variables that underlay infection processes. Fluorescence microscopy data on intracellular distributions of Salmonella enterica serovar Typhimurium (S. Typhimurium) show that, while the number of infected cells increases with time, the distribution of bacteria between cells is stationary (though highly skewed). Here, we report a simple model framework for the intensity of intracellular infection that links the quasi-stationary distribution of bacteria to bacterial and cellular demography. This enables us to reject the hypothesis that the skewed distribution is generated by intrinsic cellular heterogeneities, and to derive specific predictions on the within-cell dynamics of Salmonella division and host-cell lysis. For within-cell pathogens in general, we show that within-cell dynamics have implications across pathogen dynamics, evolution, and control, and we develop novel generic guidelines for the design of antibacterial combination therapies and the management of antibiotic resistance.


Burst at a fixed threshold
Consider first a single phagocyte which is infected with a bacteria at time t = 0. By assumption, the cell will not be further infected by bacteria from outside the cell, and the number of bacteria increases by fission at per capita rate α until the threshold N is reached. Rather than solving directly Eqns. [1a,b] from the Methods section of the paper, an explicit solution may more straightforwardly be obtained from an alternative (but equivalent) formalism where individual cells are considered before scaling up to the whole population. The Forward Equation for the probability P n that there are n bacteria in a cell at a time t after infection is dP n dt = α(n − 1)P n−1 − αnP n , n ≤ N − 1 This can be solved by standard techniques such as generating functions, subject to the initial condition P 1 (t = 0) = 1, to yield P n (t) = e −αt (1 − e −αt ) n−1 n ≤ N − 1 The cell lyses and releases N bacteria the instant that n = N , so the average rate of release of propagules at time t from a cell infected at time 0 is R(t) = αN (N − 1)P N −1 (t). ( Since all infections are the result of the release of propagules via cell burst, and these cells will themselves have been infected some time in the past, we can write a selfconsistent equation for the rate I(t) at which cells are being infected at time t: where t 0 is the time when the first cell was infected. This assumes that susceptible host cells are sufficiently in excess that each propagule successfully and instantaneously infects a new cell. Let us assume that the initial infection was far enough in the past so that t 0 →−∞, and that the infection rate increases exponentially, I = I 0 e γt (it is possible to justify this assumption by using Laplace transforms to solve Eqn. (3) explicitly, whereupon we need to to include an initial external infection source at t 0 ). Any reasonable initial state (i.e. one with a finite number of bacteria) will converge to this quasistationary state over time. Dividing eqn. (3) by I(t) gives where we have substituted for R from Eqn.
Since the denominator is a monotonically increasing function of γ α , there can be only one solution to this equation, which is clearly γ = α. This makes sense intuitively, because in a quasistationary state we would expect the total number of infected cells to be proportional to the total bacterial population size, which is itself proportional to e αt since there is no source of bacterial mortality.
The total number of infected cells with n bacteria will be t −∞ I(t )P n (t − t )dt , while the total number of infected cells will be proportional to I(t). The fraction of cells which have n bacteria will be .

Heterogeneous burst size
Suppose now that the pre-programmed threshold N at which each cell lyses differs between cells. Let D N be the burst size distribution, i.e. the fraction of cells that have threshold N (so N D N = 1). The probability that a cell contains n bacteria a time t after it was infected still satisfies Eqn. (1), but the mean rate of release of propagules is now The infection rate I(t) satisfies Eqn. (3), and again looking at the quasistationary state where t 0 → −∞ and I(t) = I 0 e γt we find As before, we find that this equation has only a single solution, α = γ, because B( γ α + 1, N − 1) is a monotonically decreasing function of γ α . The steady state distribution of bacteria between cells satisfies where the limits on the sum over N arise because P n = 0 for n > N (see Eqn. (1)). Substituting for P n from Eqn.
(1) and integrating we get Eqn. (5) allows us to compute the steady state distribution for a given distribution D N . The cases we studied were

Density-dependent lysis and fission
Let us now assume that cells lysis takes place at a given rate, rather than instantaneously at a threshold. Let this rate µ n depend on the number n of bacteria in the cell. Let us also assume that bacterial fission is density dependent, with α n being the percapita splitting rate for a bacteria in a cell containing n bacteria. While the approach of sections 1 and 2 can be used to give an explicit closed-form solution for some simple cases, we instead develop a method for computing the distribution numerically for general α n and µ n . The Forward equation for the total number M n of cells with n bacteria at time t is where δ n,1 is the Kroeneker delta which ensures that the final term (which represents the infection of susceptible cells) only applies to n = 1.
If we assume a quasi-stationary state with M n ∝ e γt Q n , Eqn. (6) becomes For n ≥ 2, this can be rearranged to give a recurrence relation for the distribution of bacteria between cells: Meanwhile, multiplying Eqn. (7) by n and summing gives an expression for γ: For the case where fission is density-independent, α n = α, Eqn. (9) reduces to γ = α, so it is possible to compute the Q n directly (initially in terms of Q 1 , then the whole distribution may be normalised to sum to unity) from Eqn. (8) in the form Note that the stationary distribution Q n is only sensitive to the ratio µn α , not to the absolute values of these rates.
If fission is density dependent, Eqns. (8) and (9) need to be solved self-consistently. One convergent procedure is to assume a value of γ, compute Q n from Eqn. (8), calculate an updated γ using Eqn. (9), and iterate until the initial and final values of γ are numerically indistinguishable. The distribution Q n is only sensitive to the ratios αn γ and µn αn , so our quoted values of µ n and γ are in units of α 0 .

Control of within-host proliferation
If there is inter-cellular mortality of bacteria, so that only a fraction β successfully reinfect, then Eqn. (6) is replaced by The threshold for control of within-host proliferation is obtained by setting dMn dt = 0. The stationary distribution for a given function α n is then Q n = α n−1 (n − 1) α n n + µ n Q n−1 for n ≥ 2. The threshold value of β is then obtained by setting n = 1 in Eqn. (11)

Equivalence between models
Since we have data for the stationary distribution of bacteria between cells, we can only distinguish between models via their predictions for Q n . We have had to assume specific parametric forms for either the burst size distribution D N or the lysis and fission functions µ n and α n . Does the difference in quality of fit between models reflect true differences between the classes of models, or is it merely an artefact of our chosen parametrisations? First, let us show that, for any heterogeneous burst model characterised by D N , there is a density-dependent lysis model characterised by µ n that gives exactly the same prediction for Q n . First, using Eqn. (5) we find the following expression for a heterogeneous burst model (for n ≥ 2): For a density-dependent lysis model with density-independent fission, Eqn. (10) can be written in the form Q n Q n−1 = (n − 1) 1 + n + µn α .
Combining these two expressions, we find Eqn. (12) allows us to compute the lysis function µ n for a density-dependent lysis model that is equivalent to a given heterogeneous burst model defined by D N . Since all the D n are positive, the µ n are also positive, so the lysis model is biologically meaningful. Therefore, for any heterogeneous burst model with a given number of parameters, we can construct an equivalent density-dependent lysis model with the same number of parameters that gives exactly the same fit. Figure 1 shows the computed lysis functions µ n equivalent to the maximum likelihood fit for the heterogeneous burst model with Gamma-distributed burst size distribution. It is also possible to show that there is a density-dependent fission model that is equivalent to each heterogeneous burst model, though the proof is a little harder. However, the fact that we can find a model of class 'A' (say) which is equivalent to any given model of class 'B' (say) does not logically mean that there is a model of class B which is equivalent to any model of class A -in other words, the set of models of class B might merely correspond to a subset of the models of class A. Let us now see if there is an equivalent heterogeneous burst model to any given density-dependent fission model, the class which gave us the best numerical fit. Considering Eqn. (5) for j = n and j = n − 1 where in the last line we have substituted Q j Q j−1 for a density-dependent fission model from equation (8), assuming density-independent lysis µ n = µ. While this formula allows us to compute D n from α n , it turns out that for some fission functions α n the burst distribution D n is sometimes negative-which means that the 'equivalent' heterogeneous burst model is not biologically valid. This is the case for the best-fit densitydependent fission model for the virulent strain with exponentially decreasing fission, as illustrated in Figure 2. Conversely, using the same method we can show that, for any given density-dependent lysis model with density-independent fission, there is indeed a valid equivalent heterogeneous burst model (we already proved the equivalence in the opposite direction in Eqn. (12)).
It is intuitively reasonable (though mathematically non-trivial) that heterogeneous burst models and density-dependent lysis models give rise to the same class of stationary distributions. Until the cell bursts, the dynamics are the same for both modelsthe difference being that the size at which burst takes place is pre-determined in the former case and emerges dynamically from the lysis function in the latter. However, for the density-dependent fission model the rate of bacterial splitting itself slows down over time, so there is no reason to expect the dynamics to be equivalent to the other models.
We conclude that, while there is a degree of equivalence between models, the models with density-dependent fission occupy a larger 'space' than the others in the sense that there is a density-dependent fission model equivalent to each heterogeneous burst model but not vice versa. It so happens that no heterogeneous burst model is able to give the excellent fit to the virulent strain found for our best density-dependent fission model, allowing us to reject the former class of model in favour of the latter.

Statistical procedure for model fitting
The output of the models is a stationary distribution Q n (X) of the number of bacteria per cell, which depends on a set of parameters X and is to be compared with observed : Burst size distribution for the heterogeneous burst model equivalent to the best fit density-independent lysis model with exponentially decreasing fission, computed from Eqn. (13). Note Q n is negative for the virulent strain when n < 6.
counts C n of the number of cells with n bacteria. If N = C n is the total number of counts, then assuming independence of observations the C n will be multinomially distributed with likelihood L({C n }|{Q n }) = N ! n Q Cn n C n ! .
We used a combination of conjugate gradient and simplex methods to iterate towards values of the parameters X which gave distributions Q n that maximised the likelihood L for our data C n . The numerical routines were written in Fortran 90, based on the algorithms in Press et al. (Numerical Recipes in FORTRAN: the art of computer programming, Cambridge University Press, available online at http://www.nr.com). The Akaike Information Criterion is calculated as AIC = −2 ln L + 2k, where k is the number of fitted parameters in the model.