An efficient moments-based inference method for within-host bacterial infection dynamics

Over the last ten years, isogenic tagging (IT) has revolutionised the study of bacterial infection dynamics in laboratory animal models. However, quantitative analysis of IT data has been hindered by the piecemeal development of relevant statistical models. The most promising approach relies on stochastic Markovian models of bacterial population dynamics within and among organs. Here we present an efficient numerical method to fit such stochastic dynamic models to in vivo experimental IT data. A common approach to statistical inference with stochastic dynamic models relies on producing large numbers of simulations, but this remains a slow and inefficient method for all but simple problems, especially when tracking bacteria in multiple locations simultaneously. Instead, we derive and solve the systems of ordinary differential equations for the two lower-order moments of the stochastic variables (mean, variance and covariance). For any given model structure, and assuming linear dynamic rates, we demonstrate how the model parameters can be efficiently and accurately estimated by divergence minimisation. We then apply our method to an experimental dataset and compare the estimates and goodness-of-fit to those obtained by maximum likelihood estimation. While both sets of parameter estimates had overlapping confidence regions, the new method produced lower values for the division and death rates of bacteria: these improved the goodness-of-fit at the second time point at the expense of that of the first time point. This flexible framework can easily be applied to a range of experimental systems. Its computational efficiency paves the way for model comparison and optimal experimental design.

where matrices C, B and A are as before. As all powers of K have the same block-triangular structure, we can denote exp(tK) as, Taking the derivative of exp(tK) with respect to t gives, d dt exp(tK) = K exp(tK), and hence expressions for each of F (t), G(t), and H(t), dG(t) dt =AG(t), and dH(t) dt =CH(t) + BG(t).
Solving for H(t) gives exp(Ct)× t 0 exp(−Cs) × B × exp(As)ds , which we can use to evaluate the second moments of the stochastic system exactly. Thus, for a given set of model parameters, we can solve the first two moments exactly.

S1.2 Derivation of Hellinger and Kullback-Leibler Divergences for Multivariate Normal Distributions
For distributions P and Q of a continuous random variable, the Kullback-Leibler divergence [1,2] is defined as The Rényi divergence of order r [3] is given by and this equals the Kullback-Leibler divergence in the limit r → 1 [4]: Another measure of divergence between P and Q is the squared Hellinger distance [5] defined We derive expressions for ∆ K and ∆ 2 H with respect to P and Q as multivariate normal distributions. In both cases, the following lemma will be used in the proofs. Lemma 1: If P and Q are multivariate normal distributions with parameters (µ 1 , Σ 1 ) and (µ 2 , Σ 2 ), respectively, then

Theorem 1:
If P and Q are d-dimensional multivariate normal distributions with parameters (µ 1 , Σ 1 ) and (µ 2 , Σ 2 ), respectively, then Proof From the limit expression of equation (1), we have therefore, from Lemma 1, We can solve the limit within equation (3) as follows. Let and consider ∂a(r)/∂r. This partial derivative can be solved through the use of the fact that Given that lim r→1 a(r) = log 1 = 0, lim r→1 r − 1 = 0 and ∂a(r)/∂r exists, we can determine the limit within equation (3) via l'Hôpital's rule: therefore, upon substituting (4) into (3), we have

Theorem 2:
If P and Q are multivariate normal distributions with parameters (µ 1 , Σ 1 ) and (µ 2 , Σ 2 ), respectively, then Proof We start with the definition of ∆ 2 H given by equation (2): But, from Lemma 1,       Figures S9 and S10 demonstrate that the sum of the divergences corresponds to the smallest mean absolute relative error in the parameter estimation, on average.

S1.3 Moment computation for radial and linear networks
Note that the Kullback-Leibler divergence could be estimated in the other direction -that is, from the observed data to the model, rather than from the model to the data as we have considered thus far. We have considered the MDE tool with this alternate parameterisation and the results were very similar, and thus, are omitted.

S1.7 Inference: Simulation Study
Parameter estimate results for all scenarios, across all observation times.  Figure S11: Box plots of MDE estimates corresponding to simulated data at τ 1 = 2, where the true parameter is denoted by the red cross. S1.8 Moments of Blood, Liver, Spleen System The first two moments (M 1 and M 2 , respectively), of the blood, liver, spleen system follow a closed system of nine differential equations. With respect to the liver and spleen, denote the rate at which bacteria: move from the blood (clearance) as c L and c S ; move back into the blood (emmigration) as e L and e S ; replicate as r L and r S ; and, are killed as k L and k S .
We can write the differential equations corresponding to the first moment in matrix form as the following: where, Similarly, we can write the second moments as: where M 1 (0) are the initial moments of the system, , and, with ∆ L = r L − k L − e L , ∆ S = r S − k S − e S , and κ = c L + c S . The solutions to these equations are given by equations (3) and (4) in the main text.