Comparing treatment strategies to reduce antibiotic resistance in an in vitro epidemiological setting

Significance Antibiotic combination therapy is a promising strategy to combat the rising problem of resistance. However, developing such strategies is hindered by the lack of an experimental system that allows testing of strategies in a realistic epidemiological context. We present an approach to test the effects of different treatment strategies in vitro in order to narrow the gap between computational and clinical studies. Using a laboratory automation system, we approximate the epidemiological population dynamics expected in a hospital setting using experimentally evolving bacterial populations. We show that as long as there is no influx of double resistance into the system, combination therapy outperforms all other tested strategies, even at concentrations lower than those used for monotherapy.


Mathematical model
The population dynamics outlined in Fig. 1a in the main text are reflected by the system of ordinary differential equations (ODEs) given by Eqns. [1][2][3][4][5]. The variables U, S, A, B and AB reflect the populations of uninfecteds, sensitive infecteds, and resistant infecteds to drugs A, B, or both, as described in the main text. The descriptions of the parameters of the model are given in Table S1.

Model Fitting and Parameter estimation
We used two different parameter estimation procedures for model fitting. The first procedure is referred to as independent estimation, where the fitting is done independently for each dataset coming from each treatment strategy. The second procedure is referred to as simultaneous estimation, where combinations of datasets of different treatment strategies are used simultaneously for fitting. Interpreting the time series data requires two types of models: a process model which is given by Eqns. 1-5, and a statistical model for the deviations between the process model and the data. The statistical model is necessary to evaluate the goodness of the fit of the process model, and thus determine what the optimal parameter values are. Following a well established tradition in model fitting, we use the least-squares method (1), and assume that the deviations between the process model and the data are normally distributed.
The Metropolis-Hastings algorithm is used for each estimation procedure, which is a Monte Carlo Markov Chain (MCMC) method to simulate multivariate distributions (2). A uniform prior distribution on the interval (0, 1) (U(0, 1)) is assumed for all parameter values, since they represent rates with no prior information. The error is assumed to be normally distributed, leading to a likelihood function of where θ, Yi,Ŷi, and σ denote the parameter vector, data point at time i, estimation of the data point at time i, and the standard deviation of error, respectively, and SSQ is the abbreviation for "sum of squares". The natural logarithm of this likelihood is given by [10] Using the Normal likelihood to fit model parameters to data thus requires an estimate for the variance of the deviation between Y andŶ . One can either estimate this from the dataset itself, or use the maximum likelihood estimate of the error variance σ 2 , which is equal to SSQ(θ; Y )/N . Plugging this equality in Eq. 10 leads to log(L(θ|Y )) = − N 2 log SSQ(θ; Y ) + C, [11] where C absorbs all the constant terms that are independent of θ, thus can be omitted for optimization purposes. As a result, the final log-likelihood function that is used for the MCMC algorithm becomes log(L(θ|Y )) = − N 2 log SSQ(θ; Y ). [12] To calculate the posterior distributions of each parameter, we ran MCMC Metropolis-Hastings algorithm with 10000 iterations per chain for 50 randomly initiated chains. To cover all the possible regions of the parameter space, parameter vectors are randomly initialized at the beginning of each chain. The first half of each chain (first 5000 iterations) is regarded as the burn-in period and discarded. After the posterior distributions for each parameter are estimated for each estimation procedure, median values of these distributions are used as the parameters of the corresponding model, and the least squares estimation Y (θ) is calculated for each treatment strategy separately. Posterior distributions including all model parameters are presented in Fig. S2, and the estimated time series are presented in Fig. 1c-h in the main text in addition to the experimental data.
Next we tested whether the estimated parameters reproduce the experimental finding that combination therapy outperforms the other strategies also in the model. Following the analysis provided in Tepekule et al. (3), we randomly sample parameter space to identify the parameter regions where particular treatment strategies succeed in reducing the number of infecteds the most. Specifically, we sampled from a uniform distribution in the open interval (0,1) for the parameters shown in Fig. S1. Due to the high dimensionality of the parameter space we used linear discriminant analysis (LDA) and projected the results onto the two principal axes of LDA (Fig. S1). We then compared the results of this analysis with the median of the posterior distributions obtained from the simultaneous estimation, and observed that the set of estimated medians fall into the region where combination therapy outperforms other treatment strategies with a high probability (Fig. S1). The observation that all parameter combinations in the posterior fall into a narrow region that is contained in the region where combination therapy outperforms the other treatment strategies, suggest that superiority of combination therapy is highly robust. Fig. S1. LDA of random sampling results. LDA classifying parameter sets according to the strategy with the highest gain in the frequency of uninfecteds per transfer. Shaded areas represent the density of the parameter sets coloured according to which strategy wins. Each treatment strategy is represented by a different colour, and the opacity of each colour is proportional to the number of parameter sets that fall into that corresponding region. The red envelope encompasses parameter estimates obtained via simultaneous estimation. 10 6 randomly sampled parameter sets are used. LD1 and LD2 are the two principal axes of the LDA. The parameter vectors are given with their relative magnitudes and counterclockwise angles. The estimated parameters (red envelope) fall into a region where combination therapy is expected to outperform the other strategies

Supplementary Methods
Estimation of mutation probability. Mutation probabilities were estimated using a fluctuation assay (4) following a protocol adapted from Kohanski et al. (5). The fluctuation test was performed in the same medium as used in the evolution experiment (Minimal salts medium with 15µg/m chloramphenicol, see main methods). Briefly, overnight cultures of the sensitive and the two single resistant strain were diluted 1 : 10000 in 50ml medium and incubated for an additional 3.5h at 37°C. Then each culture was diluted 1:3 and ten 1ml aliquots were incubated for 24h at 37°C. All cultures were then plated on solid media without selection as well as on the antibiotics to which the strain is not already resistant to measure mutation probabilities towards those drugs. Plates were counted after 24h incubation at 37°C.
Total colony forming units as well as the distribution of the number of resistant colonies on the different drugs were then used to calculate mutation probabilities using the maximum likelihood method as implemented in the R package flan (6)