Using a Sequential Regimen to Eliminate Bacteria at Sublethal Antibiotic Dosages

We need to find ways of enhancing the potency of existing antibiotics, and, with this in mind, we begin with an unusual question: how low can antibiotic dosages be and yet bacterial clearance still be observed? Seeking to optimise the simultaneous use of two antibiotics, we use the minimal dose at which clearance is observed in an in vitro experimental model of antibiotic treatment as a criterion to distinguish the best and worst treatments of a bacterium, Escherichia coli. Our aim is to compare a combination treatment consisting of two synergistic antibiotics to so-called sequential treatments in which the choice of antibiotic to administer can change with each round of treatment. Using mathematical predictions validated by the E. coli treatment model, we show that clearance of the bacterium can be achieved using sequential treatments at antibiotic dosages so low that the equivalent two-drug combination treatments are ineffective. Seeking to treat the bacterium in testing circumstances, we purposefully study an E. coli strain that has a multidrug pump encoded in its chromosome that effluxes both antibiotics. Genomic amplifications that increase the number of pumps expressed per cell can cause the failure of high-dose combination treatments, yet, as we show, sequentially treated populations can still collapse. However, dual resistance due to the pump means that the antibiotics must be carefully deployed and not all sublethal sequential treatments succeed. A screen of 136 96-h-long sequential treatments determined five of these that could clear the bacterium at sublethal dosages in all replicate populations, even though none had done so by 24 h. These successes can be attributed to a collateral sensitivity whereby cross-resistance due to the duplicated pump proves insufficient to stop a reduction in E. coli growth rate following drug exchanges, a reduction that proves large enough for appropriately chosen drug switches to clear the bacterium.


Batch-transfer protocol
For all the protocols consisting of multiple serial batch transfers referred to in the main text, with each transfer conducted once per 'season', we used the same microtitre plate reader (BioTek) to measure optical densities every 20 minutes at 600nm as a proxy for bacterial population densities in different environments (written as OD 600nm or just OD). We used 96-well plates containing 150µL of liquid in each well incubated at 30 • C to culture bacteria, shaken in a linear manner before each OD measurement was taken.
For an experiment exposing bacteria to antibiotics lasting several seasons (where each season lasts 12h or 24h) an initial inoculation was performed using an isogenic population obtained from a single colony and cultured overnight in M9 minimal media (0.2% glucose, 0.1% casamino acids) at 30 • in a shaker-incubator. At the end of each season, the same 96-pin plate replicator was used to sample the liquid volume (containing bacteria and spent medium) which was then transferred to a new plate containing fresh growth medium and antibiotics, ensuring the same environment for each replicate population was maintained. Every subsequent transfer to a fresh plate containing medium was performed using the same 96-pin replicator, we estimated the volume transferred to be 1.5µL. The so-obtained, timedependent optical densities were then imported into Matlab in order to subtract the background (determined from blank wells containing only medium) and compute the mean optical densities and other statistics. N.B.: media-only wells testing for the presence of potential contamination were used on every 96-well microtitre plate. If any showed turbidity above blank levels, the assay in question was repeated.
1.3 Live cell counts: optical density is a reasonable proxy for cell density  Fig. S1 -Data showing that OD is positively correlated with a live cell count (CFU per ml) in the plate reading devices that we use. (a) E. coli K12(AG100) was cultured in M9 and serially diluted to different OD values, these were then counted to reveal CFU values (all error bars are 95% confidence intervals of the mean, three replicates). A linear regression (blue) is shown next to data. A constrained linear regression that must pass through the datapoint (0, 0) is also shown (red). (b) The latter indicates that OD has the potential to over-estimate live cell numbers at the lowest densities.
Many of the bacterial growth and inhibition experiments described in this article require the continual measurement of bacterial population densities and this cannot be done with by-hand lab techniques, such as colony counting. We therefore use devices for which proxies of population density can be rapidly produced using automated protocols that read densities as light absorbance.
No claim is made in the article on the basis of OD data alone that a zero population density has resulted under any treatment conditions, the measurement error inherent to such devices is too great for that to be reliable. When any claim is made that no cells are present, this is always deduced from first observing an OD 600 value of below 10 −2 units after at least 12h growth, thereafter cells are counted manually following a spot test in order to determine whether cells are present and, if so, how many. Fig. S1 shows that OD is a reasonable indicator of bacterial population density: this measure correlates positively, and indeed linearly, with viable cell counts (measured in units of colony-forming units per ml).

Dose-response curves
We determined dose-responses by measuring bacterial OD 600nm dynamics at different dosages of both antibiotics, the resulting growth data was used to estimate the drug concentration necessary to achieve x% inhibition with respect to the null-antibiotic control (at drug concentrations hereafter denoted IC x ). As illustrated in Fig. S2, the approximate IC 50 values obtained for a 12h experiment were 0.04µg/ml for DOX and 6.13µg/ml for ERY. The degree of variation in these measures is illustrated using a 95% confidence interval for each dose response curve in Fig. S2

Drug interaction curves
We say that synergy is exhibited by a drug pair, at a given time, if the following holds. If B(t, D, E) denotes the time-series of bacterial densities as a function of time and dose, given two drugs denoted D and E, then when (1) B(t 0 , D 0 , 0) = B(t 0 , 0, E 0 ) (an equi-potency assumption on both drugs) (2) i(θ) := B(t 0 , θD 0 , (1 − θ)E 0 ) < B(t 0 , 0, E 0 ) for all θ between 0 and 1, (note the later equals B(t 0 , D 0 , 0) by assumption) the drugs D and E are said to synergise at basal dosages D 0 , E 0 at time t 0 . Under these assumptions, by fitting a quadratic to the interaction function i(θ), so that αθ 2 + βθ + γ ≈ i(θ), we can determine the numerical degree of synergy from the value of α. When the inequality α > 0 holds for all E 0 , D 0 and t 0 , a criterion that can be established in practise using standard numerical fitting routines and a t-test, we then say D and E synergise. This procedure, called the α-test, is detailed fully in [18]. The quadratic fits to bacterial density data in Fig. S3 establish that this statistical test is passed for ERY, DOX in the experimental described conditions described above using the strain E.coli K12(AG100).
1.6 The (n + 1)-protocol: measuring collateral sensitivity and cross resistance In order to quantify how acclimation (UK English equivalent 'acclimatisation') to one antibiotic impacts population growth after a switch to another antibiotic occurs, and to measure this as a function of the duration of the acclimation period, we evolved five replicates of an isogenic population of E. coli for n + 1 seasons at IC 70 (see Fig. S4), where each season lasts 24h, as follows. Bacteria were exposed to a constant environment (DOX or ERY) for n ≥ 3 seasons before being transferred into an environment containing the other drug: those growing in DOX were transferred into an environment containing only ERY, those growing in ERY were then treated with DOX. We compared the optical densities of the population that underwent a change in drug at the end of the n-th season to the corresponding optical density of the population growing under constant conditions. This allows us to determine whether acclimation to one drug is associated with increased sensitivity to the other. Fig. S3 -Establishing statistically significant synergy at both 12h and 24h using the α-test for synergy defined in [18] (p-values indicated in the figure) with both erythromycin and doxycycline at their respective IC50 dosages. Superimposed fits of a quadratic function is shown, the shaded region represents an estimated 95% confidence interval. The property that the fit is a convex function (and not concave) establishes the ERY-DOX synergy because a 50-50 combination reduces the optical density (indicated on the y-axis) relative to both ERY-only and DOX-only monotherapies.
(a)  The treatments implemented using the (n + 1)-protocol are illustrated in Fig. S5.  Differences in population densities produced by different sequential treatments were measured as follows. An 8season batch transfer protocol was implemented, with each season having 12h duration. At the end of each season, a 1% sample of the final population was transferred using a pin replicator into a fresh 96-well plate containing a replenished environment (containing both fresh M9 media and antibiotic).
Drug treatments consisted of the sequential deployment of DOX and ERY whereby only one drug was used within each season. The order in which the drugs were deployed defined each treatment and, in order to perform a like-withlike comparison, we introduced the constraint that all treatments were balanced: by the end of an 8-season treatment, each drug would be used for 4 seasons (ERY for 4, DOX for 4) and this is true for all treatments implemented. Thus, these sequential treatments are illustrated using the schematic in Fig. S6. Each such treatment was initially replicated three times. Different results are obtained in this protocol depending on the choices of dose for ERY and DOX. All treatments, including monotherapies (ie. treatments that use one antibiotic for the entire duration), were implemented at IC 50 and IC 70 dosages. As a control, we also implemented the 50-50 combination treatments where, to clarify what this means, if monotherapies used dosages of 3µg/ml for ERY and 0.1µg/ml of DOX, the 50-50 treatment would have 1 2 × 3µg/ml of ERY and 1 2 × 0.1µg/ml of DOX. By design, the 50-50 combination treatment must yield a lower OD than the monotherapies in the first season to be consistent with previous studies reporting ERY and DOX (and our prior data) as a synergistic drug pair [9,18]. §2 Typical growth data and the rate of adaptation Typical raw data typical are shown in Fig. S7. If no drug is deployed (yellow data) the OD increases each day with the bacterium entering stationary phase progressively earlier each transfer (a.k.a. season). Antibiotic use reduces population growth and delays entry into stationary phase. Calculating the phenotypic rates of adaptation using the measure defined in [9] for the data presented in Fig. S7 produces Fig. S8. This shows that growth rate adaption to treatment is different for each condition, being fastest for the 50-50 antibiotic treatment and slowest for the DOX monotherapy. Interestingly, the rate of adaptation for populations cultured in no-drug conditions are between these two extremes. This data shows that antibiotic inhibitory effect and the rate of adaptation are not necessarily correlated and that adaptation to a media-only environment (i.e. no drug) can be just as rapid as adaptation to antibiotic-inhibited growth.   Fig. S7. There is a positive correlation between rate of adaptation and final growth rate in the treated populations. The no-drug conditions produce higher growth rates than the drug-treated conditions and adaptation is just as rapid as in drug-treated conditions. (b) Adaptation is fastest to the synergistic 50-50 combination treatment. Importantly, the rate of adaptation to the no drug control environment is not the lowest of all the conditions tested. This data indicates no positive correlation between rate of growth adaptation and the inhibitory effect of the antibiotics because adaptation can be just as fast in the absence, as in the presence, of antibiotics. (c) Growth rate is, of course, greatest in the last season in the populations with no antibiotic. Moreover, growth rate is not lowest on the last season in the 50-50 combination treatment, despite the ERY-DOX synergy (vertical bars are s.e., n = 6). §3 Additional Data 3.1 The combination treatment eventually loses out to all sequential treatments at IC 50 The data presented in the figures of this section (Fig. S9, Fig. S10) demonstrate that a treatment producing the greatest inhibitory effect when measured over a single season need not continue to inhibit growth maximally as treatment proceeds. At these dosages, the 50-50 combination treatment is optimal at reducing population growth in season 1, but it is close in performance to the worst treatment of all by season 8. Vertical thick black lines denote the mean OD of the 50-50 combination treatment, normalised and written as zero, vertical dotted lines represent the mean OD taken over all the sequential treatments and the coloured dots show the deviation of the OD of each sequential treatment from the 50-50 combination (blue dots: DOX, green dots: ERY). The histograms represent the distributions of performances of all sequential treatments relative to the combination treatment in each season. This series of so-called forest plots show how at the beginning of the experiment all dots are to the right-hand side of the 50-50 line: their performance is worse than the multidrug combination. However, as the number of seasons increases, more treatments cross the vertical line and thus achieve higher inhibition than the 50-50 combination. By the end of the experiment all sequential treatments outperform the multidrug combination.  We first determined the ODs produced by all the balanced sequential treatments at IC 70 dosages and, as shown in Fig. S11, we sought a negative correlation between the number of switches in an 8-season treatment and the number of bacteria it produces. This was done in two ways: we measured the total number of bacteria produced by each treatment (using OD as the proxy) and we measured the final OD produced, namely at 96h. We then compared the number of drug switches with these population density measures produced by each treatment. This figure shows no evidence of the correlation we sought.  in OD from one season to the next (namely (y − x)/(x + y) where x is OD at 12h on day n and y is OD at 12h on day n + 1, expressed as a percentage) provides no evidence to reject H0 using the same dataset as used in Fig. S10. Here (y − x)/(x + y) (that, note, necessarily lies between −1 and +1) is computed for two treatment classes: 'twist' treatments occur when ERY follows DOX or DOX follows ERY from day n to day n + 1, 'stick' treatments occur when DOX follows DOX or ERY follows ERY. (b) Stratifying the same treatments further is not sufficient to reject H0.    Fig. S13. This illustrates the short-term synergy, the loss of synergy, it also highlights that these treatments lead to increasing population densities at 96h (error bars are s.e. with n = 3). Note how the population density of one replicate can eventually collapse for a given treatment whilst a different replicate can recover, leading to large betweenreplicate variations. Note also that this replicated dataset has 21 successful treatments from the full set of 48 replicates. This is a different outcome to the first set of replicates shown in Fig. S13, however the so-called 'red square' treatments behave consistently in all replicates by producing a zero CFU count at some time. §4 A whole genome sequencing analysis In order to investigate antibiotic resistance adaptation based on prior and de novo mechanisms, we performed the following analysis. We cultured bacteria with 12h seasons using a batch transfer protocol implemented in a shaken flask with 5ml of M9 liquid medium, implementing one environmental condition without drugs (called the 'no drug control') and two different antibiotic conditions. All were replicated three times. The two antibiotic treatments were the 50-50 combination at IC 50 dosages and the sequential treatment of the from 'EDEDEDED'. Measuring a 150µL sample from the flask in a plate reader at 600nm, the cumulative optical densities produced were similar for both these drug treatments as shown in Fig. S16. All replicate populations from 24h and 96h were then sequenced using the paired end technology on a Illumina 7500 machine at the Exeter sequencing service. Both single-drug monotherapies and the 'DEDEDEDE' treatments were implemented as controls, but were not sequenced (see Fig. S16).

Library preparation method
DNA was fragmented by sonication using a Biorupter for 30s on, 90s off, using low power for 10 minutes on ice. Libraries were prepared using SPRIworks cartridges for Illumina (Beckman Coulter) and Nextflex indexed adapters, with 300-600 bp size selection, amplified with 8 cycles PCR using Kapa HiFi DNA polymerase and purified using GeneRead kit (Qiagen). Concentrations were determined using a Bioanalyser 7500 DNA chip. Libraries were pooled in equimolar amounts, denatured, diluted to 6.5 pMol and clustered on a flowcell using a cBot (Illumina). 100 paired end sequencing with a custom barcode read was completed on a HiSeq 2500 using Truseq SBS v3 reagents (Illumina).

Construction of a Local Reference
To facilitate the genomic analysis we first constructed a local reference genome of E. coli K12 (AG100), constructed by modifying the publicly available annotated genome of strain MG1655 1 [19]. This choice was motivated by its close relation with AG100 and the quality of the existing annotation.
Reads were processed with fastq-mcf [2] to remove adapters from the sequencing data and to trim and filter lowquality reads. In particular, cycles with at least 1% of 'N's were removed (command-line parameter: -x 0.01). The remaining reads were mapped to MG1655 using the Burrows-Wheeler aligner BWA [13] with standard parameters. The resulting alignments were processed with Samtools 1.0 [14], with pair/trio calling enabled (command-line parameter: -T). Subsequently, alignments were sorted, artifacts and duplicates were removed, and finally the alignments was indexed. Unaligned reads were stored separately.
Structural Variations (SVs) were detected using Pindel [25]. Its pattern-growth algorithm detects breakpoints of large deletions and medium-sized insertions by identifying paired reads for which only one of the reads can be mapped to the reference. It then attempts to break the unmapped read into two and maps both shorter fragments to the reference. If successful, the breakpoints of deletions or insertions can thus be determined. We restricted the detection to SVs of maximum 2,071,552bp (command-line parameter: -x 8). The detected SVs were visually inspected with IGV [21,23]. We validated 5 deleted regions, 9 putative breakpoints and 9 indels, nucleotide polymorphisms were intentionally ignored at this step.
Using a Python script, we removed deleted regions and putative breakpoints detected by Pindel, resulting in a set of 15 non-overlapping sequences free of major SVs. Detected indels were subsequently applied to these sequences. With Samtools we extracted the reads previously mapped to removed regions or within 2/3 of the library insert size from their borders. These and the reads without a feasible alignment were assembled with Velvet [26] resulting in a set of contigs. The union of this set and the set of the 15 sequences was used to generate a set of scaffolds using SSPACE [3]. This software uses the distance information of paired-ends to assess the order, distance and orientation of the submitted contigs and combines them into scaffolds. The resulting scaffolds were extended with GAPFILLER [4] that attempts to close gaps between scaffolds using the distance information of paired-read data; 21,942 scaffolds were thus obtained.
These scaffolds were aligned to MG1655 using Mauve [20]. This software can align a set of sequences to one genome from a related one that differs from the former for the presence of SVs and only the 13 longest sequences (the shortest was 109,280 bp) could be aligned to the reference. The remaining sequences, not more than 2,375 bp each, were discarded. The 13 scaffolds were concatenated into one sequence, here called the 'intermediate' reference according to an order and orientation given by Mauve.
We used Pindel to verify if SVs were still present in the intermediate reference and to validate the multiple nucleotide polymorphisms detected in the previous execution of Pindel. No SV was observed visually with IGV, while 6 polymorphisms of 2 consecutive nucleotides each, detected even before by Pindel, were still present. We applied these polymorphisms to the intermediate reference using the script vcf-consensus of the vcf-tools utility [6]. Finally, SNPs were detected using the software VarScan.
VarScan uses both heuristic and statistical methods to call SNPs based on read depth, base quality, significance and variant frequency. We used VarScan with the following parameters: p-value threshold of 0.05 for calling variants, minimum read depth of 20 to make a call at a position, at least 8 supporting reads at a position to call variants, base quality not less than 20 at a position to count a read, minimum variant allele frequency of 0.03; frequency to call homozygote of at least 0.9 (command-line parameters: --p-value 0.05 --min-coverage 20 --min-reads2 8 --min-avg-qual 20 --min-var-freq 0.03 --min-freq-for-hom 0.9) The detected SNPs were applied to the intermediate reference with vcf-consensus and the resulting sequence was our local reference. The local reference was mostly annotated using RATT [17], transferring the annotation from MG1655 to this sequence for regions with high synteny. Using the RAST annotation server, we annotated the two longest regions with low synteny (37,193 bp and 7,166 bp, respectively).

Mapping onto local reference
All 18 replicates of three drug treatments were mapped to the local AG100 reference genome using the same method adopted for mapping the control replicate to the MG1655. As above, reads were processed with fastq-mcf and mapped to the local reference with BWA. Again, SNPs and SVs were called with VarScan and Pindel, respectively, using the same parameters previously adopted. Additionally, we used CNVnator [1] to discover copy number variations (CNVs). This software detects CNVs through an analysis of read mapping density (coverage) within different bins along the genome. A bin size of 60 was chosen for all CNVnator analyses.

SNP-detection heuristic
Let '0' denote a known wild-type allele at a given locus and suppose that Illumina sequencing produces n aligned reads covering this locus that is denoted by a sequence of alleles, {X j } n j=1 , where each X j is a Bernoulli random value that only takes the values '0' or '1', the latter symbol (taking a numerical of unity) denoting any synonymous or non-synonymous mutation. Now define Y n = 1 n n j=1 X j . Suppose that the systemic per-nucleotide error rate at each locus is , this is the probability that one Illumina read reports an incorrect allele following alignment.
A putative SNP acceptance rate of α per read (expressed as a value between 0 and 1) is here said to produce a false positive rate p that is defined to be the probability that n aligned reads reports at least α × n occurrences of the allele '1' due to Illumina read error. This is the value p := P(nY n > nα) = P(Y n > α). As each X j can be modelled as a Bernoulli trial with parameter , for sufficiently large n the probability distribution of Y n can be approximated by a normal distribution with the following mean and variance parameters (following standard notation) µ := E(Y n ) = and σ 2 := var(Y n ) = (1 − )/n. Thus following the change of variable, q = (α − µ)/(σ √ 2). Using we deduce the following approximation for the false positive SNP acceptance rate for given α: where α > . Using a value = 10 −2 [16,15] and a mean observed value for coverage of n ≈ 200 (this is representative of values given in Fig. S17), with this calculation, shows a putative variant frequency of α = 2% yields the value p ≈ 0.0776 whereas using α = 3% yields p ≈ 2.24 × 10 −3 . The requirement that p < 0.05 thus permits us only to report SNPs with an estimated frequency of 3% or greater.

WGS Results
Our analysis indicated a 412,380bp duplication (from 273,601 to 685,980) in all 6 drug treatment replicates at 96h. After a list of putative SNPs was first generated in the manner detailed above, Table S1,Table S3 and Table S4 then summarise only those SNPs for which significant longitudinal or between-treatment differences are observed using a one or two-way anova with a threshold of p < 0.1. Putative SNPs that do not pass this acceptance criterion have not been reported. Genes annotated using red font are contained within the duplicated region. For completeness, Table  S5 contains basic statistics used in the whole-genome analysis.  Table S1 are known tetracycline or erythromycin resistance mutations. 23S erythromycin resistance mutations in the database: http://goo.gl/f7WrCH The database has a tetracycline resistance mutation but we do not observe it here: http://goo.gl/0NsnJf (2) atoB is a short-chain fatty acid degradation enzyme (thiolase II) [11] in the ato operon.

Tables of SNPs
These summarise significant SNPs observed at 24h and 96h under the different treatment conditions. The mean frequency in the populations where the SNP was detected is indicated and a superscript denotes the number of replicates in which the SNP was observed. A blank entry denotes it was not detected.  Table S2 contains a list of genes within the duplicated genomic region that are known to interact with antibiotics in some way, whether through drug binding or efflux, or else are implicated in the regulation of stress pathways that may be associated with increased resistance.   gene  position  24h  96h  24h  96h  24h  96h annotation ampE 119,553 b.p. 8 (2) 6 (1) 6 (1) 13 (1) regulator; ampicillin resistance inner membrane protein duplicated region (red text denotes duplicated genes) paoC 298,353  The mathematical model featured in the main text and implemented using the parameter values in Table S6 shows that an asymmetric collateral sensitivity between ERY and DOX can arise for the following reason. The efflux pump is predicted to have different affinities for ERY and DOX, this causes different rates of selection for cells that express the pump and, also, for cells that duplicate the pump. This means different population structures are supported by the use of either drug in monotherapy.
To show this can be the basis of an ACS using the model, we first choose antibiotic supply concentrations for both drugs at their IC 50 (measured at 24h) in the model (main text, third column in Fig. 5A). This same figure from the main text then shows that the population of bacteria adapts during five days of treatment with ERY, with increases in the number of cells both expressing the efflux pump and having duplications of the pump genes.
In the mutation-selection equilibrium indicated by this figure, the population converges to a stable configuration with population densities below ones that are achieved under no drug conditions. However, a switch to DOX in the model sees a change in the relative frequency of pump duplications but, for instance, because the pump in the model has different efficiencies at removing the two drugs from within the cell, this exchange in drug sees population densities rapidly increase. In the figure, after the switch the population soon achieves densities close to those observed without the drug. The third column, bottom plot of Fig. 5A and then Fig. 5B show analogous, but opposite, effects when the drugs are exchanged in the opposite sense.
Thus there are a number of factors that can contribute to an asymmetric collateral sensitivity: different efficiencies of the pumps at effluxing the drugs, the different affinities of the drug for their targets, the different reductions in absolute fitness attributable to each drug, the mutation rates associated with the duplications and the selection for each of those mutations resulting from these factors that control the rate of sweep of the duplications. Given such physical and evolutionary asymmetries, the model exhibits both a cross resistance and a collateral sensitivity, as the terms are defined in the text. Fig. S18 shows that the population densities produced by each treatment depends, even in this theoretical model, on the way in which drugs are sequenced. Although both drugs are calibrated in the figure to an IC 50 dose with respect to a population consisting almost of entirely of wild-type cells that do not carry duplications of the pump, as treatments progress the population structures do change. This can mean, as in Fig. S18, that some sequential treatments produce more bacteria than the 50-50 combination, whereas some produce fewer. Indeed, of the treatments shown in Fig. S18, it is a sequential treatment that produces fewest cells of all.   Table S6 give rise to optical densities in the mathematical model (at 0.2% glucose) that are substantially larger than those reported in the experimental data throughout this article. This is not an error, rather the stated parameters were determined by calibrating the mathematical model in the main text against optical density data read at 600nm in 384 well plates [18]. With similar liquid volumes as the 96 well plates we use here, the path length taken for the absorbance measurement is necessarily larger and, therefore, similar population densities will result in very different optical densities between that article and this. Therefore, an optical density conversion factor of approximately 1 4 OD units is required to compare data from this article and the mathematical model defined in [18] that we use here.