Combination Therapy Outperforms Other Strategies.

We compare alternative treatment regimens based on their effect on 1) maximizing the frequency of uninfected individuals and 2) minimizing the frequency of individuals infected by a resistant strain. The first criterion can be applied to all six treatment arms. The second one, however, does not allow comparison with control using no antibiotics, as the control trivially prevents the evolution of resistance. The resistance phenotype of the bacterial population in a well was measured by transferring small aliquots from the well onto four different agar plates containing either no antibiotic, only antibiotic A or B, or both. An A-resistant phenotype would thus grow on the no-drug and the A-containing plate, but not on the B- and AB-containing plates. Additionally, we assessed resistance phenotype by growth (i.e., optical density at 595 nm [OD₅₉₅] > 0.1) in the experimental cultures.

For the first scenario (scenario $$ \varnothing $$), we assumed that there was no influx of antibiotic-resistant infections from outside (i.e., the wells representing newly admitted patients contained either no bacteria or sensitive bacteria only). Hence, all resistance observed evolved in response to treatment in the hospital. Patients were discharged and admitted with a rate of 0.2 per day. The rate of infection was 0.3. The time course of 40 successive days of treatment for all six treatment arms is shown in Fig. 1 C–H. With regard to the treatment effect on the frequency of uninfecteds, we found that combination therapy was most effective, as the entire infected population was cleared throughout the course of the experiment. Cycling and mixing lead to a lower level of uninfecteds, but still substantially higher than both monotherapies. Because of the strong bottleneck at every transfer (a dilution of $$ 1/166$$), low-frequency, nongrowing phenotypes, such as persisters, are not expected to be maintained in the experiment. With regard to the treatment effect on resistance, again, combination therapy outperformed the other strategies. The resistance frequency under cycling and mixing was similar, and lower than under the monotherapies. Surprisingly, resistance to both antibiotics did not emerge and spread, despite the fact that there was at least one combination of two single-point mutations (i.e., gyrA S83 L and rpsL K43 R), each giving high-level resistance to one antibiotic, that conferred high-level resistance to both antibiotics at low fitness cost (18). Sequencing the resistance-determinant regions of gyrA and rpsL revealed mostly different mutations at those positions (SI Appendix, Table S2). Trindade et al. (18) report both positive and negative epistatic interactions between mutations at these positions. This raises the possibility that these interactions prevent the evolution of double resistance in our experiment.

The findings, that combination therapy outperforms the other strategies, and that cycling and mixing have similar effects, are in good agreement with the computer simulations of Tepekule et al. (12). Moreover, the observation that cycling and mixing have similar effects is also in agreement with the outcomes of a recent multicenter RCT (17).

Parameter Estimates Vary across Different Treatment Arms.

Fig. 1 C–H also shows the population dynamics of a model that is fitted simultaneously to the time course of all treatment arms. The model is an extension of the model used by Tepekule et al. (12) (Materials and Methods and SI Appendix). Fitting the model independently to each treatment improved the fit by permitting the same parameter to have different posterior distributions for different treatment arms. Ideally, we would expect similar posterior distributions for a given parameter across those treatment arms for which the parameter can be estimated, since building an epidemiological model fundamentally assumes that the dynamics of different treatment arms can be explained by using the same set of processes. The identifiability of the parameters for the different treatment arms is shown in a Venn diagram (SI Appendix, Fig. S4). However, fitting the model separately to each treatment revealed substantial differences for some key parameters between the alternative treatment regimens (Fig. 2). Although the same processes (e.g., de novo evolution of resistance in response to drugs) occurred in the monotherapies and cycling or mixing, parameter estimates obtained from the different treatment arms differed substantially. This implies that the established approach of using a general model succeeds in explaining the dynamics of different treatment outcomes up to a certain level (Fig. 1 A–F), but does represent an oversimplification by excluding treatment-specific processes.

Fig. 2.

Posterior distributions (estimated probability distribution for a given parameter [Materials and Methods ]) obtained from fitting the model to different parts of the experimental data. Mono A, mono B, combination, cycling, and mixing refer to fitting the model to data only from these treatment arms, whereas simultaneous refers to the overall fit to the time course of all treatment arms simultaneously. A, B, and C show, respectively, the posterior distributions of the parameters $$ {\tau }_A$$, $$ {\tau }_B$$, and $$ {\nu }_B$$ describing the clearance rate after treatment with antibiotic A, clearance rate after treatment with antibiotic B, and the rate of de novo emergence of B resistance for the sensitive strain. Posterior distributions of the parameters that cannot be identified from a given treatment arm are shown as thin box plots and equal to their prior distribution, which is uniform between zero and one. Posterior distributions, including all model parameters, are presented in SI Appendix, Fig. S5.

Cycling and Mixing Perform Better than Predicted Based on Monotherapies.

Although the model gives an overall good fit to all treatment arms, it systematically overestimates the frequency of uninfecteds in monotherapy A and underestimates it for cycling. The differences in key parameters (Fig. 2) indicate that the population dynamics of multidrug strategies, such as cycling and mixing, cannot be predicted accurately from measurements based on the monotherapies only. Importantly, using the parameter estimates from the monotherapies to model the outcome of cycling and mixing leads to an overestimation of the frequency of resistants and an underestimation of the frequency of infecteds, implying that both cycling and mixing perform better than would be anticipated on the basis of the effect of the monotherapies (SI Appendix, Fig. S3).

Preexisting Double Resistance Negates Benefit of Combination Therapy.

For scenario $$ \varnothing $$ (admission of uninfected and sensitive infected patients only), combination therapy outperforms all other treatment strategies. But how does combination therapy perform when resistant patients are also admitted? To address this, we studied two further scenarios with progressive levels of resistant inflow into the hospital. In scenario I, also, patients were admitted who carry bacteria resistant to a single antibiotic. In scenario II, additionally, patients carrying double resistance were admitted. For the three scenarios, the different treatment regimens were compared on the basis of the frequencies of uninfecteds and resistants, respectively, over transfers 9 to 12 (Fig. 3; full time series are shown in SI Appendix, Figs. S7–S8). Under scenarios $$ \varnothing $$ and I, combination therapy significantly outperformed the other treatment regimens, with regard to both maximizing the frequency of uninfecteds and minimizing the frequency of resistants. In scenario II, combination therapy was no longer significantly different from cycling and mixing, and all strategies selected for high levels of double resistance. In summary, combination therapy offers a benefit, as long as there is no influx of double resistance, and all strategies fail otherwise.

Fig. 3.

Frequencies of uninfected and resistant populations for different treatment strategies in different resistance inflow scenarios. Points show population phenotype frequency averaged over transfers 9 to 12 for four biological replicates for scenarios $$ \varnothing $$ (no inflow of resistance), I (inflow of both single resistances), and II (inflow of double resistance). Bars show the median. Samples not sharing a letter are significantly different (generalized linear hypothesis test across treatment strategies, P < 0.05; ANOVA tables and all P values can be found in SI Appendix, Tables S4–S9). R, resistant populations; U, uninfected (sum of phenotypes A, B, A/B, and AB; Fig. 1). The full time series are shown in SI Appendix, Figs. S6–S8.

Benefits of Combination Therapy Do Not Arise only from Higher Dosage.

Does the advantage of combination therapy stem from increased drug dosage? More specifically, does combination therapy outperform the other strategies because it is based on using each antibiotic at the respective concentration used in the single drug therapies? To address this question, we performed another series of analogous experiments for combination therapy, lowering the drug concentration from 2$$ \times $$ MIC, as used in scenarios $$ \varnothing $$, I, and II, to 1$$ \times $$ and 0.5$$ \times $$ MIC for both drugs (plus a control of no drugs). We measured the effect on the emergence of double resistance in two ways: 1) growth on solid medium containing 10$$ \times $$ MIC, indicating high-level resistance; and 2) growth in liquid culture at the respective drug concentration used for treatment indicating low-level resistance. The results show that combination therapy prevents the emergence of both high- and low-level double resistance, as long as the concentration of both antibiotics are 1$$ \times $$ MIC or higher (upper left triangles, Fig. 4). This finding argues that it is not the effect of drug dosage, but the probability of simultaneous acquisition of multiple mutations that underlies the success of combination therapy. For lower concentrations, low-level, but not high-level, resistance evolves (lower right triangles, Fig. 4). Note, however, that the concentrations used in liquid are below MIC for at least one antibiotic, where even bacteria sensitive to these antibiotics are expected to grow. Given that at 1$$ \times $$ MIC for both antibiotics, combination therapy prevents the emergence of resistance; this indicates that the benefit of combination therapy does not arise only from the higher combined concentration compared to the other treatment strategies. This is another possible benefit of combination therapy: Combining antibiotics might enable the reduction of the concentration of potentially toxic antibiotics, while still minimizing the likelihood of resistance evolution.

Fig. 4.

Mean frequencies of double-resistant populations after five transfers at different concentrations of nalidixic acid and streptomycin. The upper left triangle indicates frequency of double resistance to a high concentration of drugs used on plate (10$$ \times $$ MIC), while the lower right triangle indicates double resistance to the drug concentrations used in liquid culture. Full time series of population phenotype frequencies can be found in SI Appendix, Fig. S9.

Strains, Drugs, and Media.

All strains were grown in minimal salts (MS) medium [MS contains 1 g/L (NH₄)₂SO₄, 3 g/L KH₂PO₄, and 7 g/L K₂HPO₄ supplemented with 0.4 mM MgSO₄, 3.33 nM FeSO₄, 1.2 mM Na₃C₆H₅O₇, and 0.8 ng/L thiamine], and 0.2% glucose was added as a carbon source. Furthermore, 15 $$ \mathrm{\mu }$$g/mL chloramphenicol was added to maintain the used mutants and to prevent external contaminations. Solid medium was prepared the same way, but supplemented with 1.5% (weight/volume) agar.

Streptomycin was used at a concentration of 12.5 or 100 $$ \mathrm{\mu }$$g/mL in liquid and solid medium, respectively. Nalidixic acid was used at a concentration of 20 or 40 $$ \mathrm{\mu }$$g/mL in liquid and solid medium, respectively. These concentrations correspond to 2$$ \times $$ MIC in liquid and 10$$ \times $$ MIC (as measured by using a standard broth microdilution assay) on solid medium for both antibiotics.

All strains used were derived from E. coli K-12 MG1655. A gyrA S83 L (C248 T) or a rpsL K43 R (A129 C) mutation was added by using single-stranded DNA recombineering (26) and subsequently transduced alone or in combination into a fresh background by using P1 transduction (27, 28). The mutations conferred resistance to nalidixic acid or resistance to streptomycin, respectively, at relatively low fitness cost (18) and were selected on plates containing 40 $$ \mathrm{\mu }$$g/mL nalidixic acid and/or 100 $$ \mathrm{\mu }$$g/mL streptomycin. Additionally, and following the same methodology, a fluorescent marker linked to a chloramphenicol resistance gene (cat-YFP) from E. coli MDS42(YFP) (29) was transduced into the single- and double-resistant strains and selected on plates containing 25 $$ \mathrm{\mu }$$g/mL chloramphenicol. In order to increase mutation supply, the mutS knockout from E. coli JW2703 (30) was also transduced into the strains and selected on plates containing 25 $$ \mathrm{\mu }$$g/mL kanamycin, resulting in the strains used in the experiment. All used genotypes were verified by resequencing of the modified genes. A list of all strains and their genotype can be found in SI Appendix, Table S3.

Evolution Experiments.

In order to test the effect of different treatment strategies on the evolution of antibiotic resistance, we set up large-scale, serial passage evolution experiments using a Tecan Evo 200 automated liquid handling system (Tecan) with an integrated, automated incubator (Liconic STX100, Liconic) and a Tecan Infinite F200 spectrophotometer (Tecan).

Our experiments were designed to simulate a hospital or hospital ward. We modeled beds as wells in a 384-well plate (Greiner catalog no. 781186) and patients as growth media in those wells. All beds were always occupied, i.e., filled with a total of 50 $$ \mathrm{\mu }$$L of growth medium, allowing for a maximum population size of about $$ 10^7$$ CFU. One 384-well plate contained four replicate units with 94 patients each. The remaining eight wells were used as sterile controls. We serially passaged every well every day into new plates, incorporating the following processes: 1) treatment; 2) admission and discharge; 3) transfer; 4) infection and superinfection; and 5) phenotyping. Before each transfer, a custom R script computed the necessary pipetting steps using the parameters outlined below and prepared pipetting worklists. For scenario $$ \varnothing $$, we initially planned the experimental evolution experiment for 40 transfers. However, we decided to extend the experiment and switch the two monotherapies and to continue the experiment for another 12 transfers. We finished the experiment after 52 transfers; however, we here focus on the 40 transfers before the drug switch. The complete time series for scenario $$ \varnothing $$ can be found in SI Appendix, Fig. S6. Scenario I and scenario II were run for 13 and 22 transfers, respectively. Because these experiments were largely dominated by preexisting resistance, and because of the trajectories observed in the scenario $$ \varnothing $$, we decided to stop these scenarios earlier, as the dynamics had already run its course (SI Appendix, Figs. S7 and S8). Every 10 transfers, and at the end of each experiment, 35 $$ \mathrm{\mu }$$L of 40% glycerol was added to the assay plates, and the plates were stored at −80 °C.

Treatment.

First, we prepared new 384-well plates by filling all wells with 40 $$ \mathrm{\mu }$$L of growth medium. We then added 5 $$ \mathrm{\mu }$$L of diluted antibiotic stock solution (20× final concentration, or no antibiotic for control wells and wells without treatment) to all wells, according to the treatment strategy. We used a separate 384-well plate for every treatment strategy. All wells, except the sterile controls, were treated, irrespective of their infection status. We considered the following treatment strategies: 1) no treatment; 2) monotherapy with nalidixic acid; 3) monotherapy with streptomycin; 4) combination therapy with both drugs; 5) cycling, where the used antibiotic was changed every other transfer; and 6) mixing, where one of the antibiotics was randomly assigned to each patient. For the initial setup of the first 384-well plate, no antibiotics were added. While all plates were prepared and treated at the same time, the following steps were conducted in series to minimize the difference in incubation time between the plates.

Admission and Discharge.

The experimental design also reflected the influx of patients from a community outside the hospital. This community was assumed sufficiently large, such that it did not change during the course of the experiment. Experimentally, this was implemented as follows. A certain fraction of wells (0.2 for all experiments reported here) were randomly chosen for admission and discharge. These wells were exempt from the transfer to the new plate (see Transfer section). Instead, these wells were inoculated with 5 $$ \mathrm{\mu }$$L from replicate overnight cultures, or sterile medium for input of uninfecteds, that was prepared from frozen stocks every day for the following day. This represents patient influx from the outside community. The fraction of the different phenotypes (uninfecteds, sensitive infecteds, single-resistant A, single-resistant B, and double-resistant AB) constituting the overall influx differed for the different scenarios considered. These proportions were $$ \left(0.15,0.85,\mathrm{0,0,0}\right)$$ for scenario $$ \varnothing $$, $$ \left(0.21,0.57,0.11,0.11,0\right)$$ for scenario I, and $$ \left(0.21,0.52,0.11,0.11,0.05\right)$$ for scenario II, respectively. All wells that were not chosen received 5 $$ \mathrm{\mu }$$L of sterile medium to equalize culture volumes.

Transfer.

Next, all wells that had not been chosen for admission and discharge were diluted into the new plate at the same position. This was done by using a custom pintool developed in-house that allows excluding certain positions by removing certain pins (i.e., the pins corresponding to the wells that are chosen for admission and discharge). The pintool transferred a fixed volume of about 0.3 $$ \mathrm{\mu }$$L of the culture for a dilution of approximately 1/166 at each transfer (about 5 × 10⁴ CFU). For the initial setup, this step was skipped.

Infection and Superinfection.

To implement infection, we randomly chose a fraction of wells (0.3 in all experiments) as the sources of infection. The targets of infection were chosen randomly with replacement from the remaining wells. Approximately 0.3 $$ \mathrm{\mu }$$L of the chosen source wells were then transferred from the previous plate to the chosen target well on the new plate by using the custom pintool for a dilution of approximately 1/166 (approx. 5 × 10⁴ CFU). The newly prepared plates were then moved to the incubator and incubated at 37° C and 95% relative humidity for approximately 24 h.

Phenotyping.

After the newly prepared plates were placed in the incubator, an aliquot of the cultures from the previous transfer was spotted onto agar plates containing no antibiotic, nalidixic acid, streptomycin, or both antibiotics. The agar plates were then also moved into the incubator. At the beginning of the next transfer, agar plates were fetched from the incubator, and a picture of each plate was taken. Pictures were analyzed for growth (growth/no growth) by using the Pickolo Software package (SciRobotics), as well as manual inspection of each picture by using a custom R script.

To measure growth of cultures in the 384-well plates, OD₅₉₅ of each well was measured at the beginning of each transfer. We defined growth in 384-well plates as reaching an $$ O{D}_{595}>0.1$$ after incubation.

From these two measurements, we defined six possible phenotypes: U (uninfected), for populations that showed no growth on any of the agar plates and in liquid culture; S, (sensitive) for populations that grew only in the absence of antibiotics on agar or in liquid culture; A or B, for populations that grew on agar plates or in liquid culture in the presence of nalidixic acid or streptomycin, respectively; A/B, for populations that grew in the presence of nalidixic acid and streptomycin (on agar plates or in liquid culture), but not in the presence of both drugs at the same time, indicating a mixed population; and AB, for populations that grew both in the presence of both single drugs and the combination. All other populations were classified as E, for erroneous phenotypes.

Sequencing.

To identify common resistance mutations, we picked three random populations from time points between 20 and 40 from all four replicates of the monotherapy treatments, as well as cycling and mixing. Populations were inoculated from frozen samples and grown for 24 to 48 h, and DNA was isolated by using the Wizard Genomic DNA purification kit (Promega). We sequenced the main resistance-determinant regions for streptomycin and nalidixic acid, which are in gyrA and rpsL, respectively (18). The regions were amplified and sequenced by using the same primers as used in Trindade et al. (18).

Statistical Analysis.

To compare the three different scenarios studied, we averaged the population phenotype frequency over transfers 9 to 12 for all four biological replicates for each scenario. Transfers 9 to 12 were chosen because they are the latest timepoints for which we have data for all three scenarios. The effect of different treatment strategies on the frequencies of uninfected and resistant populations in the three different scenarios was tested by using an ANOVA, followed by a post hoc generalized linear hypothesis test for a model that included both treatment strategies and phenotype as a main effect, as well as their interaction. ANOVA tables, as well as all tested linear hypotheses and the corresponding test statistics, can be found in SI Appendix, Tables S4–S9. All statistical analyses were performed in R 4.0.2 (31) using the packages tidyverse (32), multcomp (33), and multcompView (34).

Mathematical Modeling.

Adapting the mathematical model by Tepekule et al. (12), we used a compartmental model to describe the population dynamics observed in the evolution experiment. The mathematical model is described in greater detail in SI Appendix. The model considered each well as a patient and ignored the population dynamics of bacteria within the wells. Hence, at each time point, each well belonged to one of the five compartments: U, S, A, B, or AB. The phenotype A/B, which represents a mixed population of single resistants, was not considered in the model, since the frequency of A/B in the experiment was negligible compared to the other phenotypes. The epidemiological processes included in our mathematical model are analogous to those implemented during the experiments, which were admission and discharge of the patients, infection, superinfection, clearance due to successful treatment, and de novo emergence of resistance during incubation (Fig. 1A). Each process was associated with a rate, and these rates were translated into model parameters of the corresponding system of ordinary differential equations (ODEs). The system of ODEs and the model parameters with their corresponding descriptions and units are provided in SI Appendix.

Parameter Estimation and Model Fitting.

We estimated the parameters of the ODE-based model by fitting it to the time course of the frequency of the different phenotypes observed for scenario $$ \varnothing $$, as this experiment had the largest number of successive days of treatment. Since our experimental results are highly replicable (Fig. 1 C–H), we used the mean frequency of phenotypes averaged over the four replicates. To fit the mathematical model to the experimental data, we implemented the Metropolis–Hastings algorithm (35) and used the least-squares method (36), assuming that the error between the model outcomes and the experimental data were normally distributed (see also SI Appendix). By minimizing this error, we obtained posterior distributions for each parameter, which were then used to simulate the mathematical model. A uniform prior distribution on the interval $$ \left(\mathrm{0,1}\right)$$ ($$ \mathcal{U}\left(\mathrm{0,1}\right)$$) was assumed for all parameter values. We implemented two different fitting procedures, where 1) the parameters were estimated separately for each treatment by minimizing the error independently for different treatment arms (independent estimation); and 2) the parameters were estimated simultaneously by minimizing the error over the concatenated time course of all treatments (simultaneous estimation). Posterior distributions of the model parameters for independent and simultaneous estimations are given in SI Appendix, Fig. S5.