Collateral sensitivity is contingent on the repeatability of evolution

Antibiotic resistance represents a growing health crisis that necessitates the immediate discovery of novel treatment strategies. One such strategy is the identification of sequences of drugs exhibiting collateral sensitivity, wherein the evolution of resistance to a first drug renders a population more susceptible to a second. Here, we demonstrate that sequential multi– drug therapies derived from in vitro evolution experiments can, in some cases, have overstated therapeutic benefit – potentially suggesting a collaterally sensitive response where cross resistance ultimately occurs. The evolution of drug resistance need not be genetically or phenotypically convergent, and where resistance arises through divergent mechanisms, the efficacy of a second drug can vary substantially. We first quantify the likelihood of this occurring by use of a mathematical model parametrised by a set of small combinatorially complete fitness landscapes for Escherichia coli. We then verify, through in vitro experimental evolution, that a second–line drug can indeed stochastically exhibit either increased susceptibility or increased resistance when following a first. Genetic divergence is confirmed as the driver of this differential response through targeted and whole genome sequencing. These results indicate that the present methodology of designing drug regimens through experimental collateral sensitivity analysis may be flawed under certain ecological conditions. Further, these results suggest the need for a more rigorous probabilistic understanding of the contingencies that can arise during the evolution of drug resistance.

, the x − y plane represents the genotypes and the height of the landscape above this plane represents fitness. Two evolutionary trajectories, both starting from a wild-type genotype (yellow circle), are shown. These trajectories diverge at an evolutionary saddle point (blue triangle) and terminate at distinct local optima of fitness (purple pentagon, green star). As the saddle point exists, evolutionary trajectories need not be repeatable. B) Schematic landscapes for a potential follow-up drug are shown, the collateral response can be (i) always cross-resistant, (ii) always collaterally sensitive or (iii) dependent on the evolutionary trajectory that occurs stochastically under the first drug. C) A potential evolutionary branching point in the TEM gene of E. coli identified in the fitness landscape for cefotaxime derived by Mira et al. [21].
known genotypes to form a fitness landscape. However, to derive fitness landscapes through this 91 method, the number of strains that must be engineered grows exponentially with the number of  A schematic of the model used to derive collateral response. Sequential mutations are simulated to fix in the population until a local optimum genotype arises. The fitness of this resultant genotype is compared to the fitness of the wild-type genotype for each of the panel of antibiotics.
B) The landscape for ampicillin derived by Mira et al. [21] represented as a graph of genotypes. Arrows indicate fitness conferring mutations between genotypes represented as nodes. Cyan nodes indicate genotypes from which evolution can stochastically diverge, grey nodes indicate genotypes from which there is only a single fitness conferring mutation. Squares indicate local optima of fitness with colour indicating the ordering of fitness amongst these optima (darker red indicates higher fitness). Two divergent evolutionary trajectories, in the sense of the model shown schemaically in A, are highlighted by coloured arrows. C) The best, worst, most likely and mean tables of collateral response derived through stochastic simulation of the experimental protocol. Columns indicate the drug landscape under which the simulation was performed and rows indicate the follow-up drug under which the fold-change from wild-type susceptibility is calculated. Bar charts indicate, for each labelled first drug, the number of follow-up drugs exhibiting collateral sensitivity (blue) or cross resistance (red) in each case. • Pairs with reported collateral sensitivity induce collateral resistance with probability p=0.52 extent to which the fitness 3 to ensure our model is a Markov Chain. In this case the step t to t + 1 can n to take some fixed arbitrary time. istribution of a population at time t is related to its initial distribution, µ (0) , by Markov chain is absorbing we know that there exists some k such that k [73]. Consequently, we know that the matrix P * = lim t→∞ P t (6) d in fact this limit is reached after only finitely many matrix multiplications. ively see that this limit is reached in finitely many steps note that all paths the Markov chain are strictly increasing in fitness and there are only finitely tes (corresponding to the genotypes). Thus a given initial population ion µ (0) will converge to a stationary distribution µ * after a finite number of our model. Furthermore, if P * is known then we compute the stationary * Peaks s case the step t to t + 1 can ts initial distribution, µ (0) , by xists some k such that (6) any matrix multiplications. y steps note that all paths and there are only finitely en initial population µ * after a finite number of compute the stationary ng to ensure that the disease re the effects of applying * for the associated fitness arkov chain we can pective. In particular, as the Chain. In this case the step t to t + 1 can is related to its initial distribution, µ (0) , by w that there exists some k such that t the matrix only finitely many matrix multiplications. n finitely many steps note that all paths sing in fitness and there are only finitely s). Thus a given initial population y distribution µ * after a finite number of own then we compute the stationary sufficiently long to ensure that the disease , we can explore the effects of applying he matrices P * for the associated fitness ynamics in a Markov chain we can algebraic perspective. In particular, as the l is a Markov Chain. In this case the step t to t + 1 can rbitrary time. tion at time t is related to its initial distribution, µ (0) , by rbing we know that there exists some k such that we know that the matrix reached after only finitely many matrix multiplications. it is reached in finitely many steps note that all paths strictly increasing in fitness and there are only finitely the genotypes). Thus a given initial population to a stationary distribution µ * after a finite number of re, if P * is known then we compute the stationary is applied for sufficiently long to ensure that the disease y equilibrium, we can explore the effects of applying considering the matrices P * for the associated fitness volutionary dynamics in a Markov chain we can ocess from an algebraic perspective. In particular, as the sure our model is a Markov Chain. In this case the step t to t + 1 can e some fixed arbitrary time. ion of a population at time t is related to its initial distribution, µ (0) , by v chain is absorbing we know that there exists some k such that Consequently, we know that the matrix t this limit is reached after only finitely many matrix multiplications. e that this limit is reached in finitely many steps note that all paths kov chain are strictly increasing in fitness and there are only finitely responding to the genotypes). Thus a given initial population will converge to a stationary distribution µ * after a finite number of el. Furthermore, if P * is known then we compute the stationary s ovided a drug is applied for sufficiently long to ensure that the disease es evolutionary equilibrium, we can explore the effects of applying equentially by considering the matrices P * for the associated fitness ncoding the evolutionary dynamics in a Markov chain we can volutionary process from an algebraic perspective. In particular, as the case the step t to t + 1 can s initial distribution, µ (0) , by xists some k such that This fitness function represents a genotype-phenotype map in the simplest sense -assigning to each 68 genotype a single real-valued fitness. Gillespie [1983,1984] showed that if the mutation rate u and 69 population size M of a population satisfy Mu log M << 1, and if we assume that each mutation 70 is either beneficial or deleterious, then each beneficial mutation in the population will either reach 71 fixation or become extinct before a new mutation occurs. Further, selection will be sufficiently Strong selection weak mutation This fitness function represents a genotype-phenotype map in the simplest sense -assigning to each 68 genotype a single real-valued fitness. Gillespie [1983,1984] showed that if the mutation rate u and 69 population size M of a population satisfy Mu log M << 1, and if we assume that each mutation 70 is either beneficial or deleterious, then each beneficial mutation in the population will either reach 71 fixation or become extinct before a new mutation occurs. Further, selection will be sufficiently Engineered E. coli ) and Ham(i, j) = 1 ,

Strong selection weak mutation
al neighbors , This fitness function represents a genotype-phenotype map in the simplest sense -assigning to each 68 genotype a single real-valued fitness. Gillespie [1983,1984] showed that if the mutation rate u and 69 population size M of a population satisfy Mu log M << 1, and if we assume that each mutation 70 is either beneficial or deleterious, then each beneficial mutation in the population will either reach 71 fixation or become extinct before a new mutation occurs. Further, selection will be sufficiently 72 strong that any deleterious mutation will become extinct with high probability and hence we may 73 assume that this always occurs. In the case that  we do not allow deleterious mutations to fix in the population.

105
Using this Markov Chain we can explore the possible evolutionary trajectories of a p

Predicting Effective Drug Sequences Through in vitro Experiments
d in the case of the fusidic acid isolates more rger. Exceptions to this trend were observed llin, ciprofloxacin-ampicillin, and fusidic in isolates where IC90 values were only value. Increased resistance differed among to the same drug(s) and in some cases this nsiderable ( fig. 3). We attributed the differithin a given drug(s) group to be the result hanges acquired by the isolates through cid-amikacin isolates (antagonistic interactary data S1, Supplementary Material online) increase in resistance improvement followed s adapted to doxycycline-ciprofloxacin (synon, supplementary data S1, Supplementary . Isolates evolved to ciprofloxacin-ampicillin tion, supplementary data S1, Supplementary had the least resistance improvement, an the WT MIC value. These results contrast ports based on sub-MIC adaptations, which Exposure 3 re of WT S. aureus was used to inoculate microtiter Three replicate populations were recreated for each and then used to inoculate the next concentration tion. phenotype map, implicit in their results is the observation that the environment has a minor role in Big Bang tumors (Fig. 1). New clones arising after the establishment of the early tumor likely have little effect on the fitness of the presence of abnormally motile cells early in the tumor's development. Tumors with these mutation patterns, which are described as "born to be bad, " have increased risk of both invasive growth and metastatic spread.

The genotype-phenotype map
An ecological view of cancer has emerged in recent years, one that explicitly considers cancer as much more than a collection of mutated cells and embraces a more dynamic dialog between tumor and host [6][7][8] . Critical to this view are interactions between tumor cells, between tumor and stroma, and between the tumor and its environment. Evolution and ecology are intimately entwined through mutation and selection, and both have a central role in tumor progression. However, we must understand not only the identity of the mutations that drive evolution but also the mechanisms through which these mutations manifest themselves in phenotypic change, that is, the genotype-phenotype map ( Fig. 1). This mapping is not one to one but many to many and is fundamentally the junction at which both genes and the environment meet to produce phenotypes. Recent work has shown that the complexity inherent within the genotype-phenotype map is responsible for the difficulty in predicting evolution 9 ; many genotypes produce identical phenotypes, and many phenotypes can emerge from a single genotype 10 .    Figure 3.4. The remainder of the scatter plots are presented in Appendix 1.
not well founded although drugs within a class, for example the cephalosporins, show more correlation that those between groups. An ideal pair of drugs for use in an alternating fashion would have a high negative correlation such that the evolution of resistance under one drug would induce sensitivity to a second. Unfortunately, no such drug pair exists and we must employ more sophisticated methods to identify viable sequential drug strategies. We next performed an in silico derivation of tables of collateral response, or collateral sensitivity matrices (CSMs), by simulation of evolution with the model described by Equation 3.4. These simulations mirror the experimental techniques used to derive empirical collateral response, for example those used by Imammovic and Sommer [131] to determine drug cycling protocols. Evolutionary trajectories in each drug fitness landscape, f x , were stochastically simulated from the wild-type starting genotype (g 0 = 0000) by sampling the associated Markov chain defined by P x (Equation 3.4). The simulation was terminated when the evolutionary trajectory encountered a local optimum genotype, g * x . The fitness of this final genotype in the second drug landscape, f y , was then recorded and collateral response was calculated as These collateral response measures we collated to form a table.
We can count the total number of CSMs that can be generated through this simulation. There exist 3 landscapes with only one peak accessible from the genotype g 0 = 0000, 6 in which two peaks are accessible, 4 in which three peaks are accessible and 2 in which four peaks are accessible. Assuming that for each landscape, f x , evolution is simulated from g 0 a single time to determine g * x and then the collateral sensitivity of g * x in each of the other landscapes, f y is recorded. Then there exist 1 3 × 2 6 × 3 4 × 4 2 = 82944 Drug: l pair of drugs for use in an alternating fashion hat the evolution of resistance under one drug nately, no such drug pair exists and we must y viable sequential drug strategies.
of tables of collateral response, or collateral evolution with the model described by Equarimental techniques used to derive empirical Imammovic and Sommer [131] to determine ries in each drug fitness landscape, f x , were tarting genotype (g 0 = 0000) by sampling the ion 3.4). The simulation was terminated when optimum genotype, g * x . The fitness of this final as then recorded and collateral response was that can be generated through this simulation. accessible from the genotype g 0 = 0000, 6 in hree peaks are accessible and 2 in which four landscape, f x , evolution is simulated from g 0 ollateral sensitivity of g *    Figure 3.4. The remainder of the scatter plots are presented in Appendix 1.
not well founded although drugs within a class, for example the cephalosporins, show more correlation that those between groups. An ideal pair of drugs for use in an alternating fashion would have a high negative correlation such that the evolution of resistance under one drug would induce sensitivity to a second. Unfortunately, no such drug pair exists and we must employ more sophisticated methods to identify viable sequential drug strategies. We next performed an in silico derivation of tables of collateral response, or collateral sensitivity matrices (CSMs), by simulation of evolution with the model described by Equation 3.4. These simulations mirror the experimental techniques used to derive empirical collateral response, for example those used by Imammovic and Sommer [131] to determine drug cycling protocols. Evolutionary trajectories in each drug fitness landscape, f x , were stochastically simulated from the wild-type starting genotype (g 0 = 0000) by sampling the associated Markov chain defined by P x (Equation 3.4). The simulation was terminated when the evolutionary trajectory encountered a local optimum genotype, g * x . The fitness of this final genotype in the second drug landscape, f y , was then recorded and collateral response was calculated as These collateral response measures we collated to form a table.
We can count the total number of CSMs that can be generated through this simulation. There exist 3 landscapes with only one peak accessible from the genotype g 0 = 0000, 6 in which two peaks are accessible, 4 in which three peaks are accessible and 2 in which four peaks are accessible. Assuming that for each landscape, f x , evolution is simulated from g 0 a single time to determine g * x and then the collateral sensitivity of g * x in each of the other landscapes, f y is recorded. Then there exist 1 3 × 2 6 × 3 4 × 4 2 = 82944

Drug holidays stochastically induce collateral s with few conserved motifs
In the clinical setting, drug holidays have been suggested as resistance, as resistance may not be preserved throughout t substantial time period in which no drug is given, after the ad prior to the administration of the next. This drug holiday m sequencing protocol, but is often neglected in experimental a alike. It is therefore of critical importance to assay the stabilit those between groups. An ideal pair of drugs for use in an alternating fashion igh negative correlation such that the evolution of resistance under one drug ensitivity to a second. Unfortunately, no such drug pair exists and we must phisticated methods to identify viable sequential drug strategies. rformed an in silico derivation of tables of collateral response, or collateral ices (CSMs), by simulation of evolution with the model described by Equae simulations mirror the experimental techniques used to derive empirical nse, for example those used by Imammovic and Sommer [131] to determine otocols. Evolutionary trajectories in each drug fitness landscape, f x , were imulated from the wild-type starting genotype (g 0 = 0000) by sampling the ov chain defined by P x (Equation 3.4). The simulation was terminated when trajectory encountered a local optimum genotype, g * x . The fitness of this final second drug landscape, f y , was then recorded and collateral response was response measures we collated to form a table. nt the total number of CSMs that can be generated through this simulation. ndscapes with only one peak accessible from the genotype g 0 = 0000, 6 in s are accessible, 4 in which three peaks are accessible and 2 in which four sible. Assuming that for each landscape, f x , evolution is simulated from g 0 determine g * x and then the collateral sensitivity of g * x in each of the other s recorded. Then there exist 1 3 × 2 6 × 3 4 × 4 2 = 82944 d. Unfortunately, no such drug pair exists and we must to identify viable sequential drug strategies. derivation of tables of collateral response, or collateral ulation of evolution with the model described by Equathe experimental techniques used to derive empirical se used by Imammovic and Sommer [131] to determine ry trajectories in each drug fitness landscape, f x , were ild-type starting genotype (g 0 = 0000) by sampling the P x (Equation 3.4). The simulation was terminated when ed a local optimum genotype, g * x . The fitness of this final ape, f y , was then recorded and collateral response was ponse of Y to X = log 2 f y (g * x ) f y (g 0 ) e collated to form a table. of CSMs that can be generated through this simulation. one peak accessible from the genotype g 0 = 0000, 6 in n which three peaks are accessible and 2 in which four t for each landscape, f x , evolution is simulated from g 0 hen the collateral sensitivity of g * x in each of the other re exist 2 6 × 3 4 × 4 2 = 82944 Drug: d an in silico derivation of tables of collateral response, or collatera SMs), by simulation of evolution with the model described by Equa lations mirror the experimental techniques used to derive empirica r example those used by Imammovic and Sommer [131] to determin ls. Evolutionary trajectories in each drug fitness landscape, f x , wer ed from the wild-type starting genotype (g 0 = 0000) by sampling th ain defined by P x (Equation 3.4). The simulation was terminated when tory encountered a local optimum genotype, g * x . The fitness of this fina d drug landscape, f y , was then recorded and collateral response wa nse measures we collated to form a table. total number of CSMs that can be generated through this simulation pes with only one peak accessible from the genotype g 0 = 0000, 6 in accessible, 4 in which three peaks are accessible and 2 in which fou Assuming that for each landscape, f x , evolution is simulated from g mine g * x and then the collateral sensitivity of g * x in each of the othe ded. Then there exist 1 3 × 2 6 × 3 4 × 4 2 = 82944 y Imammovic and Sommer [131] to determine ories in each drug fitness landscape, f x , were starting genotype (g 0 = 0000) by sampling the tion 3.4). The simulation was terminated when optimum genotype, g * x . The fitness of this final as then recorded and collateral response was that can be generated through this simulation. accessible from the genotype g 0 = 0000, 6 in hree peaks are accessible and 2 in which four landscape, f x , evolution is simulated from g 0 ollateral sensitivity of g * x in each of the other We exhaustively explored all evolutionary trajectories in the 15 small antibiotic landscapes to identify potentially divergent collateral response. We found a total of 82944 unique CSMs. The most common CSM occurs with probability 0.0023.

186
We note that X7 exhibits an increase in resistance to cefotaxime without any associated 187 genomic alterations. Similarly X1, X5, X9 and X12 exhibit mutations, but none that are known 188 to be associated with antibiotic resistance. Thus, we can infer that physiological adaptation or 189 epigenetic adaptation is also driving resistance to cefotaxime.     replaced by a fitter (as determined by the landscape) neighbouring genotype (defined as any 283 genotype whose Hamming distance from the population genotype is equal to one). This process 284 is stochastic and the likelihood of a genotype, j, replacing the present population genotype, i, is 285 given by Where no such fitter neighbour exists, the process is terminated. The value of r determines the 287 extent to which the fitness benefit of a mutation biases the likelihood that it becomes the next 288 population genotype. We take r = 0, corresponding to fixation of the first arising resistance termination at a local optimum of fitness, say g * . Simulated collateral response was calculated 294 as the fold difference between g 0 and g * in a second fitness landscape.

295
The code used to implement the model, produce the figures and analyse the experimental 296 data is available upon request and will be made publicly available upon publication. replicates and drugs. Note the assumption that we never erroneously take a measurement that 319 differs from the true MIC by greater than a factor of two. This is justified by noting that in no 320 instance do the maximum and minimum MICs measured in our analysis differ by greater than 321 4× (Supplementary Table X). ( where δ denotes the Kronecker delta function. By observation, the mle for each m d i is given by the median of x d i,1 , x d i,2 and x d i,3 , except in the case that two of these values are precisely 4× the other, in which case the mle is the mid-point between the maximum and minimum. Letting r denote the number of replicate/drug combinations in which all three measurements equal the mle, s denote the number in which 2/3 measurements equal the mle, t the number in which 1/3 equal the mle and u the number in which 0/3 equal the mle. Then the mle for p is given by p = s + 2t + 3u 3(r + s + t + u) .
This identity can be verified by first principles (by taking the derivative of the likelihood function) 323 but is also quite intuitive -it is simply the proportion of measurements that differ from the 324 inferred mle for the MIC. In our experiment, r = 338, s = 196, t = 11 and u = 4, which yields an 325 mle for the measurement error rate of p = 0.14.

326
The full data set, along with the inferred MIC values, are presented in Supplementary  Affairs.