Experiment design and platform

This study aims to gain insights into bacterial evolution under fluctuating antibiotic stress by the design of a series of long-term (24 days) laboratory evolution experiments. Experiments were designed to expose bacterial cultures to single antibiotic stress and a range of configurable alternating antibiotic stresses. The drug type and the concentration of the drug applied were either switched (on/off) or oscillated over a range of variable time periods (Fig. 1a). To achieve this type of programmable environmental drug condition we designed and constructed an automated morbidostat platform7, application of which enabled assessment of the development of antibiotic resistance (measured as changes in susceptibility) as a function of multi-drug composition, relative concentrations and time (Supplementary Fig. 1). In the morbidostat, the antibiotic concentration was automatically controlled by a custom algorithm7; when a bacterial culture reaches a certain density, 1 ml of a growth medium containing an antibiotic was administered into 12 ml of a bacterial culture, while a medium containing no antibiotic was added at low cell density (Methods). The antibiotic concentration in the administered growth medium was at least 10 times higher than the minimum inhibitory concentration (MIC) of a wild-type strain. Typically, antibiotics were added at least a few times once the cell density exceeds a threshold level, and thus the antibiotic concentration was increased to the range of drug concentrations that select for drug-resistant mutants (so-called ‘mutant-selection window’33). The use of the automated cell culture system was important not only in delivering programmable dosing, but the morbidostat system enabled maintenance of bacterial populations in exponential growth phase. Thus the platform was expected to provide good phenotypic reproducibility between parallel experiments7 whilst enabling genetically diverse populations to co-exist34,35. Note that we here defined a wild-type Escherichia coli MG1655 strain as a susceptible strain, in contrast to the clinical definition in which a strain can be treated with an antibiotic at the recommended dosage. Resistant strains were defined as any (evolved) bacterial samples whose MICs exceed that of the wild-type strain.

Figure 1: Evolution of antibiotic resistance under antibiotic cycling condition. (a) Experimental design for bacterial evolution under cyclic exposure to two antibiotics. Wild-type E. coli was cultured in the morbidostat system with dynamic antibiotic concentrations. The Every t days (where t=1, 3 and 6), the drugs were switched in an alternating manner and the cell cultures were sampled daily and stored at −80 °C. After a total of 12 or 24 days the MICs of the collected samples were determined. (b) Evolutionary trajectories of KAN and POL resistance (green and yellow curves, respectively) under single antibiotic condition. Trajectories of parallel experiments were shown as pale lines and the average of parallel experiments was in solid line with empty circles. The resistance levels were measured as MIC fold change relative to the wild-type, calculated by log 2 (MIC i /MIC WT ), where MIC i and MIC WT are the MICs of the i-th day sample and wild-type E. coli, respectively. (c) Trajectories of KAN and POL resistance when cycled with 3 day interval. Pale lines indicate trajectories of parallel experiments and the average was shown as solid line with filled circles. The pale coloured background indicates the antibiotics used. (d,e) evolutionary trajectories of CHL and POL resistance (blue and yellow) under single antibiotic condition and cycling with 3 day interval, respectively. Similarly, for (f,g) CHL and NIT resistance (blue and purple), and (h) and (i) KAN and NAL resistance (green and orange). Samples sizes are n=3 for single antibiotic condition and n=2 for antibiotic cycling (biological duplicates). Full size image

Bacterial evolution under alternating antibiotic stresses

Starting with an isogenic E. coli strain, a series of 24-day experiments in the morbidostat system was performed using two alternating antibiotic stresses (Fig. 1a), as well as 12-day experiments under a single antibiotic stress. Antibiotics were chosen from different classes so that they have different targets (Table 1), which minimize the risk of developing potential cross resistance between antibiotics12,13,29. To assess the reproducibility of the evolutionary trajectories, experiments were duplicated for each antibiotic cycling pair and triplicated for each single antibiotic. Experiments with reversed order of antibiotic cycling were also performed in duplicate to evaluate the effect of order of cycling. These experiments revealed that the bacterial populations exhibited two evolutionary patterns in the development of drug resistance under the cycled conditions as shown in Fig. 1. The first group developed the resistance to only one of the cycled antibiotics, whereas the second group became resistant to both antibiotics. We refer to these antibiotic-resistant conditions of the bacterial populations as a ‘single drug’ and ‘multi-drug’ resistant state, respectively. Within the second category, two different patterns were observed: An oscillation30 with a gradual increase of baseline (Fig. 1g) and a lowered maximum resistance (Fig. 1i). In the former case, the rates of adaptation to chloramphenicol (CHL) and nitrofurantoin (NIT) were both delayed (the steepness of fitted logistic curves k=0.4±0.03 and 0.19±0.01 day−1, mean±s.e.m., respectively) compared to the case of single antibiotic stress (Fig. 1f, k=0.79±0.13 and 0.55±0.05 day−1, respectively). In the latter case, the drug resistance was developed rapidly and at a similar rate to the single antibiotic case k=1.7±0.12 and 31.9±3.12 day−1 for kanamycin (KAN) and nalidixic acid (NAL), respectively, for cycling case and k=7.11±6.01 and 5.06±1.65 day−1 for single antibiotic case. However, the maximum level of NAL resistance was significantly lowered (P<0.001, Welch’s two sample t-test. Supplementary Fig. 4a), while that of KAN resistance was comparable to the single antibiotic case. Evolutionary trajectories were reproducible between parallel experiments and the evolution experiments with reversed cycling order of antibiotics showed similar evolutionary patterns (Supplementary Fig. 2a).

Table 1 List of antibiotics used in this study. Full size table

While the development of multi-drug resistance was expected, the results of cycling with polymyxin B (POL) shown in Fig. 1c,e were surprising. This is because, depending on the antibiotic combinations, either the development of POL resistance or the counteracting antibiotic was suppressed. As a result, the bacterial population was resistant to only one antibiotic. In contrast, the resistance was developed in all cases when cultured under single antibiotic stress (Fig. 1b,d,f,h). To gain a better understanding of this behaviour, the bacterial population was challenged with various antibiotic combinations (Supplementary Fig. 3). Here, the antibiotics were cycled over an interval of 1 day because shorter intervals tended to delay the development of the antibiotic resistance (Supplementary Fig. 4b–i). The results showed that cycling with POL often drove the bacterial population to a single drug-resistant state, or kept the population in that state for a long period (Supplementary Fig. 3a). In contrast, cycling with the other antibiotics ended up in a multi-drug-resistant state in all the cases (Supplementary Fig. 3b).

Bacterial evolution under single antibiotic stress

We hypothesized that the single drug-resistant state with POL might be related to the evolutionary patterns of resistance under single antibiotic condition. Indeed, the development of POL resistance under single antibiotic stress exhibited a unique evolutionary pattern with a relatively long ‘silent phase’ where the population does not develop drug resistance (up to 6 days) followed by sharp increase as shown in Fig. 1b and Supplementary Fig. 5. As such, the evolutionary patterns of the other antibiotic resistance profiles can be categorized into two types: (1) The rapid development of NAL, KAN and rifampicin resistance increasing to the maximum level with little or no silent phase, and (2) a gradual increase in resistance for CHL, NIT, ampicillin and tetracycline over time7,8. Cycling POL with the antibiotics in the former category tends to suppress the resistance to POL, while populations developed POL resistance at a relatively early stage when cycled with the antibiotics in the latter category.

Collateral sensitivity profile

We also investigated the cross resistance and collateral sensitivity profiles12,13,29 using the final day samples from the evolution experiments under a single antibiotic condition (Fig. 2). This showed that POL resistant strains showed collateral sensitivity to many antibiotics, indicating that having POL resistance tends to make bacterial cells more susceptible to other antibiotics than wild-type strains. On the other hand, strains that have already evolved resistance to the other antibiotics often developed cross resistance to the other antibiotics or remained neutral. While these results were in general consistent with previous studies on collateral sensitivity (Supplementary Fig. 12)29,36, some collateral sensitivity profiles (for example, KAN) were less evident compared to previous cases. This may be due to different selection strengths during evolution between previous and our cases because the collateral sensitivity/cross resistance profile is known to be dependent on the selection strength36.

Figure 2: The collateral sensitivity and cross resistance profile of resistant strains. MICs of the final day samples from the evolution experiments under single antibiotic conditions were tested with all eight antibiotics. Blue and red grids indicate that the antibiotic pair is collateral sensitive and cross resistance, respectively. The numbers in the grids are log 2 -transformed MIC fold change. The values shown in the matrix were rounded. See Supplementary Table 2 for the original data. Full size image

Possible mechanisms of drug resistance suppression

The multi-drug and single drug evolution experiments demonstrated various evolutionary patterns depending on the types and combinations of antibiotics. While antibiotic cycling strategies that eliminated bacteria with sublethal antibiotic dosages have been reported10, our case maintained exponential growth yet suppressed the emergence of antibiotic resistance (Fig. 1c,e). These various patterns of evolutionary trajectories make a stark contrast with a previous laboratory evolution study on antibiotic cycling where monotonic increase in drug resistance was consistently observed9. To provide possible explanations for the diverse evolutionary patterns, we first consider the growth conditions in the previous and our studies.

With the serial transfer method used in the previous study9, a bacterial population experiences exponential and stationary phases during overnight culture. While the growth rate is a main factor for resistance evolution in exponential phase, other factors such as cell–cell interaction and resource competition would come into play in stationary phase. Evolution using the serial transfer method would thus be slower than the case when the exponential phase was maintained (that is, morbidostat) because the bacteria need to adapt to both conditions in exponential and stationary phases. Antibiotic cycling would further slow the adaptation process because of collateral sensitivity, additional fitness costs by resistance to a second antibiotic and other factors9.

In contrast, the rate of resistance evolution in the morbidostat is primarily determined by the growth rate, in particular the growth rate of a new mutant relative to that of the population8. This would explain the evolutionary patterns under single antibiotic stress as the growth rate (that is, fitness) of an antibiotic-resistant mutant may vary depending on the antibiotics used. For example, if only one single mutation is required to become strongly resistant to an antibiotic, such resistant mutant will quickly sweep through the population in the presence of the antibiotic and resistance at the population level rises sharply. On the other hand, if a fitness gain by a mutation is relatively small, the evolution will slowly increase. These mechanisms may explain the rapidly or gradually increasing evolutionary patterns (for example, KAN or CHL, respectively, Supplementary Fig. 5).

The evolution under single POL stress showed a unique pattern with a relatively long silent phase that was not observed in the other antibiotics tested. A simple explanation for the pattern could be to assume that multiple mutations are required for POL resistance. However, it seems implausible because the reproducible evolutionary trajectories of POL resistance (for example, Fig. 1e) cannot be explained considering the inherent stochasticity of genetic mutation. In addition, as described below, the whole-genome sequence data did not support the explanation (Supplementary Table 1). Instead, we speculate a two-step adaptation mechanism of POL resistance because it has been known that there are non-genetic resistance mechanisms to POL, such as Lipid A modification37 and heteroresistance38. Recent studies indicated that heteroresistance to POL was achieved by upregulations of putrescine synthesis and YceI protein without any genetic mutations. In addition, release of the molecules from highly resistant subpopulations of heteroresistant bacteria protects less resistant bacteria from POL and other antibiotic stresses39,40. This evidence suggests that the E.coli populations in this study might have used the non-genetic mechanisms to cope with POL stress at the initial stage of evolution when the drug concentration was low. This could then have switched to genetic mechanisms, that is, beneficial mutations conferring POL resistance as the drug concentration increased, which may explain the step-like evolutionary pattern of POL resistance.

Based on the single drug evolutionary patterns, we then consider the evolutionary patterns of the antibiotic cycling, in particular KAN-POL (Fig. 1c) and POL-CHL (Fig. 1e) cycling. In the former case, a bacterial population may first acquire resistance to KAN as it was demonstrated to occur rapidly during the initial stage of evolution experiments (Fig. 1b). Considering that having KAN resistance generally incur fitness costs, this may then have significantly delayed evolution to acquire additional resistance to POL by extending the silent phase, especially because POL resistance comes with large fitness costs37. Indeed, in a similar case with NAL-POL cycling (Supplementary Fig. 3a), the emergence of POL resistance was largely delayed, suggesting that POL resistance appeared after fitness cost of NAL resistance was compensated by secondary adaptations41,42. This would also suggest that, if the KAN-POL cycling was continued for a longer period, POL resistance may appear at a later stage of the experiment. The latter case of CHL-POL cycling, in which CHL resistance was completely suppressed, can be explained by the reversible nature of CHL resistance (which was also observed with POL resistance, as described below). In Supplementary Fig. 2, it was observed that developed CHL resistance was reversed back to the susceptible level when the antibiotic stress was switched to POL. This indicates that CHL resistant mutants appeared during a period of CHL stress were perished when the antibiotic was switched to POL. If the cycling rate was more frequent, it would keep the CHL resistance level low because CHL resistant mutants would be removed before they become a majority in the population.

Mathematical model of bacterial multi-drug evolution

From above results, we suspected that three factors were involved in the evolutionary patterns under cycled antibiotic stress, that is, (1) the duration of silent phase, (2) the rate of adaptation and (3) cross resistance and collateral sensitivity profiles. To study the obtained results from a broader perspective, we derived a simple model incorporating the three factors to describe the evolutionary patterns (Supplementary Methods). It should be noted that we employed a phenomenological model rather than commonly adopted population genetics models8,42,43. This is because previous studies revealed that, despite diverse underlying genetic alterations, bacterial evolution for stress resistance can exhibit remarkably similar phenotypic trajectories7,8,44,45,46. Additionally, from a prediction point of view, it would be helpful if a model is based on experimental parameters that can be measured relatively easily. Thus, we here developed a theoretical model relying only on relative MIC levels.

The model was based on two assumptions: (1) The development of bacterial resistance to an antibiotic A, termed as R A , is constantly promoted during exposure to antibiotic A. However, the resistance level eventually saturates due to fitness costs47, as we observed above. (2) Increased R A negatively affects resistance to another antibiotic B, or R B largely if the antibiotics are a collateral sensitive pair12,13,29. Otherwise it has small effect (neutral or cross resistant pair). Conceptually, the model can be described as a system with positive autoregulation and double negative feedback loops (Fig. 3a). The model has three key parameters, the duration of silent phase θ, the rate of adaptation α and the collateral sensitivity/cross resistance coefficient β. These model parameters were determined based on the fitting parameters of experimental results under single antibiotic stress.

Figure 3: Theoretical models of multi-drug resistance evolution. (a) Schematic diagram of the model consisting of positive autoregulation and double negative feedback loops. (b) Phase space of multi-drug evolutionary dynamics with POL and KAN. Two solid circles indicate stable steady states, corresponding to single and multi-drug-resistant states, respectively. A broken circle is an unstable steady state. MICs were normalized to one by the maximum MIC of single antibiotic evolution. (c) Phase space of multi-drug evolutionary dynamics with NIT and CHL. Only one stable steady state, multi-drug-resistant state, exists in this case. (d–f) Simulated evolutionary trajectories of bacterial resistance with POL and KAN, POL and NIT, and CHL and NIT, respectively. Antibiotics were cycled with 1 day interval. (g–i) Experimental data collected from the morbidostat system. Full size image

The geometric structure of the model shown in Fig. 3b,c explains the origin of the unique evolution patterns of POL resistance. The nullclines for POL resistance (yellow curve in Fig. 3b) and KAN resistance (green) intersect at three points, two of which are stable steady states (indicated as solid circles). Thus, in theory, the evolution of bacterial resistance to POL and KAN can result in either a multi-drug-resistant or single drug-resistant state. However, in the case of cycling with KAN, the bacterial population developed resistance to KAN first due to the rapid adaptation rate (determined by α) and POL resistance did not develop at all despite antibiotic cycling (Fig. 3d,g). In contrast, when POL was cycled with an antibiotic with slow adaptation rate, such as NIT, the evolution ends up in the multi-resistant state (Fig. 3e,h). This is because POL resistance rapidly increased before the system reached a single antibiotic state due to the slow adaptation rate of NIT. The geometric structure is unique to the case with POL which has a large silent phase (determined by θ). In contrast, the nullclines have only one intersection when the antibiotics other than POL were cycled (Fig. 3c). This means that the evolution under antibiotic cycling always results in a multi-drug-resistant state in this system (Fig. 3f,i). Indeed, the cycling experiments using antibiotics other than POL always resulted in a multi-drug-resistant state so far as we tested (Supplementary Fig. 3b).

Evolutionary trajectories simulated by the model were in good agreement with experimental results in terms of the final drug-resistant states (Fig. 4) as well as the evolutionary patterns (Supplementary Fig. 6). In Fig. 4a,b, the normalized drug resistance level for one of the cycled antibiotics was plotted against that for another antibiotic. For example, KAN-POL cycling where KAN and POL are represented as drug #1 and #2, respectively, was mapped on the x axis because only KAN resistance was developed while POL resistance was completely suppressed. In general, the results can be categorized into two groups, that is, the single and multi-drug-resistant states, as observed above. The former cases (for example, KAN-POL) were plotted on the axes, while the latter ones (for example, NAL-KAN) were in the upper right part of the plot. To illustrate the predictability of the theoretical model, the normalized predicted MIC values were plotted against the normalized observed MIC values (Fig. 4c,d). Overall, the model predictions showed good correlations with the experimental results (r2=0.44 and 0.64, respectively). We also constructed a null model of the mathematical model as a comparison. In the null model, all the collateral sensitivity/cross resistance coefficient term β were set to zero. Comparison of predictions by the null model and experimental results showed poor correlations (r2=0.03 and 0.20. Supplementary Fig. 9). As the null model has no interaction between antibiotic resistance due to β=0, both drug resistance levels approached to the maximum at the rate determined by α. This result illustrates that the interaction between α and β is crucial for diverse evolutionary patterns and hence demonstrates the predictability of the mathematical model.

Figure 4: Comparison of experimental results and model predictions. (a) Scatter plot of observed MICs of the first drug against that of the second drug from the final day samples of cycling experiments with 1 day interval. Note that a label for each point indicates the order of antibiotics cycled. For example, KAN-POL means that KAN was used on the first day (denoted as ‘Drug #1’ in the axis label) and POL was on the second day (‘Drug #2’), and then cycled with the order. Due to this notation, KAN-POL and POL-KAN are plotted on the x and y axes, respectively, even though only KAN resistance was developed in both cases (Supplementary Fig. 3a). (b) Scatter plot of predicted MICs from the final time point of simulated results. Note that the effect of fitness cost alleviation for NAL and POL cycling was incorporated in the results. (c,d) Scatter plots of predicted MICs against observed MICs for the first drug and the second drug, respectively. Note that the same data used for a and b are used here. Full size image

It should be also noted that there were some cases in which the model and experiments did not match. For example, CHL-POL cycling resulted in a single drug resistance (POL resistance only) while the model predicted that the bacterial population would result in multi-drug-resistant states (Supplementary Figs 3a and 7a). As discussed in the previous section, an additional effect of bacterial resistance, such as the reversible CHL resistance, might be involved in the experimental pattern of antibiotic cycling and need to be incorporated in the model to reproduce the behaviour accurately. In the case of cycling using NAL and POL, experimental results showed that POL resistance appeared at a later stage of the experiment. However, the theoretical model predicted a single drug-resistant state (NAL resistance only). This suggests that there was some additional adaptation process that alleviated the fitness costs for NAL resistance occurred during the cycling experiment (for example, compensatory mutations41,42). Indeed, the model reproduced the experimental results when such effect was taken into account by reducing the collateral sensitivity/cross resistance coefficient β to zero at a certain point during the experiment (Supplementary Fig. 6b and Supplementary Methods).

We derived an approximate analytical solution of the model, which describes that the normalized maximum resistance to antibiotic A is determined as in the simplest form (Supplementary Methods). This means that the maximum resistance level is determined by the ratio of cross resistance/collateral sensitivity profile and the rate of adaptation. Such effect of α and β was indeed seen in the experimental results (Supplementary Fig. 3): Antibiotic cycling with collateral sensitive pair (that is, large β) lowered the maximum resistance level compared to that with cross resistant pair (small β), for instance, CHL and KAN, and CHL and NAL pairs, respectively. The equation also indicates that if antibiotic A with large α is cycled with antibiotic B with small α, the evolution results in on the final day, which can be confirmed in the experimental results.

Reversing drug resistance by modulating antibiotic stresses

One of the important implications in the model is that the evolution of bacterial drug resistance can be viewed as a dynamical system because the model proposed here is a simple deterministic system. Such dynamical systems view of biological systems was previously proposed in the context of stem cell differentiation and drug resistance evolution of caner cells48,49. In the current context, this view suggests that drug-resistant states of bacterial population can be directed from one state to another if antibiotic stresses are modulated externally. In fact, it was theoretically indicated that antibiotic cycling with a pair of synergistic drugs can select susceptible bacteria while eliminating drug-resistant ones from a population20. To verify this experimentally, we examined the reversibility of evolved antibiotic resistance as it can be considered as a state transition between resistant states. In particular, we here mainly focused on POL resistance because the cycling with this antibiotic has two resistant states (Fig. 3b) and also because it is one of the ‘last-resort’ antibiotics37. The model predicted that an evolved resistance can be removed from a bacterial population if antibiotic cycling is switched to a single antibiotic stress (Supplementary Fig. 7a). To confirm this experimentally, we challenged bacterial populations from the cycling experiments with consecutive exposure to counteracting antibiotics. Results showed that the resistance can be indeed reversed even from the last day of antibiotic cycling at the expense of increased resistance to the counteracting antibiotic (Fig. 5a). It should be also noted that, not only single, but also multi-drug resistance in a bacterial population can be reversed by exposing to a counteracting antibiotic (POL-NIT case in Supplementary Fig. 7b). The reversibility did not depend on the evolved resistance or counteracting antibiotics used (Supplementary Figs 2b and 7b). However, there were also the cases when an antibiotic resistance was not reversed, especially the cases when a population was evolved under a single antibiotic stress (Supplementary Fig. 7c).

We examined the mechanism behind the evolved resistance and reversibility, and identified conditions that contributed to this phenomenon: Genetic and phenotypic heterogeneity in a bacterial population. First, we sequenced the whole-genome of POL resistant mutants evolved under single or cycling POL stress. While a number of mutations were found in the evolved strains (Supplementary Table 1), we identified mutations in basS and secD genes as key mutations that confer POL resistance (Fig. 5b and Supplementary Fig. 8). BasSR (also called PmrAB) two component system regulates the lipopolysaccharide modification pathways and is known to confer POL resistance37. SecD, a part of the Sec protein translocase complex, belongs to the resistance-nodulation-cell division family of multi-drug exporters50 and is known to play an auxiliary role in antimicrobial peptide resistance51. The sequencing results indicated that the reversal of an antibiotic resistance was possible when there were genetic variations in basS or secD gene within a population (Fig. 5b and Supplementary Fig. 8). In contrast, the reversal was not observed in the genetically homogeneous population in terms of basS or secD allele (Supplementary Fig. 8). This result suggests that the cycling of antibiotics prevented selective sweep and thus maintained the genetic heterogeneity in the population even after 24 day of antibiotic exposure. This result was supported by the fact that the populations evolved under a single antibiotic stress did not show the reversal of drug resistance as the population was genetically homogeneous. This genetic diversity was also reflected at the phenotype level. We observed MICs of individual isolates were diverse and differ from MIC of a whole population (black circles and a cross in Fig. 5c)34,35. These genetic and phenotypic diversities would partly account for the reversal of evolved resistance because the exposure to a counteracting antibiotic selects less POL resistant and more CHL resistant cells due to collateral sensitivity.