Exploiting evolutionary non-commutativity to prevent the emergence of bacterial antibiotic resistance

The increasing rate of antibiotic resistance and slowing discovery of novel antibiotic treatments presents a growing threat to public health. Here, we develop a Markov Chain model of evolution in asexually reproducing populations which is an encoding of the Strong Selection Weak Mutation model of evolution on fitness landscapes. This model associates the global properties of the fitness landscape with the algebraic properties of the Markov Chain transition matrix and allows us to derive general results on the non-commutativity and irreversibility of natural selection as well as antibiotic cycling strategies. Utilizing this formalism, we analyse 15 empirical fitness landscapes of E. coli under selection by di.erent beta-lactam antibiotics and demonstrate that the emergence of resistance to a given antibiotic can be both hindered and promoted by di.erent sequences of drug application. Further, we derive optimal drug application sequences with which we can probabilistically ‘steer’ the population through genotype space to avoid the emergence of resistance. This suggests a new strategy in the war against antibiotic.therapy.resistant organisms: drug sequencing to shepherd evolution through genotype space to states from which resistance cannot emerge and by which to maximize the chance of successful therapy. Background: The increasing rate of antibiotic resistance and slowing discovery of novel antibiotic treatments presents a growing threat to public health. Previous studies of bacterial evolutionary dynamics have shown that populations evolving on fitness landscapes follow predictable paths. In this article, we develop a general mathematical model of evolution and hypothesise that it can be used to understand, and avoid, the emergence of antibiotic resistance. Methods and Findings: We develop a Markov Chain model of evolution in asexually reproducing populations which is an encoding of the Strong Selection Weak Mutation model of evolution on fitness landscapes. This model associates the global properties of the fitness landscape with the algebraic properties of the Markov Chain transition matrix and allows us to derive general results on the non-commutativity and irreversibility of natural selection as well as antibiotic cycling strategies. Utilizing this formalism, we analyse 15 empirical fitness landscapes of E. coli under selection by different β-lactam antibiotics and demonstrate that the emergence of resistance to a given antibiotic can be both hindered and promoted by different sequences of drug application. We show that resistance to a given antibiotic is promoted in 61.4%, 68.6% and 70.3% of possible orderings of single, pair or triple prior drug administrations, respectively. Further, we derive optimal drug application sequences with which we can probabilistically ‘steer’ the population through genotype space to avoid the emergence of resistance. Conclusions: Our model provides generalisable results of interest to theorists studying evolution as well as providing specific, testable predictions for experimentalists to validate our methods. Further, these results suggest a new strategy in the war against antibiotic-therapy-resistant organisms: drug sequencing to shepherd an evolving population through genotype space to states from which resistance cannot emerge and from which we can maximize the likelihood of successful therapy using existing drugs. While our results are derived using a specific class of antibiotics, the method is generalisable to other situations, including the emergence of resistance to targeted therapy in cancer and how species could change secondary to changing climate or geographical movement.


Abstract
The increasing rate of antibiotic resistance and slowing discovery of novel antibiotic treatments presents a growing threat to public health. Here, we develop a Markov Chain model of evolution in asexually reproducing populations which is an encoding of the Strong Selection Weak Mutation model of evolution on fitness landscapes. This model associates the global properties of the fitness landscape with the algebraic properties of the Markov Chain transition matrix and allows us to derive general results on the non-commutativity and irreversibility of natural selection as well as antibiotic cycling strategies. Utilizing this formalism, we analyse 15 empirical fitness landscapes of E. coli under selection by di↵erent beta-lactam antibiotics and demonstrate that the emergence of resistance to a given antibiotic can be both hindered and promoted by di↵erent sequences of drug application. Further, we derive optimal drug application sequences with which we can probabilistically 'steer' the population through genotype space to avoid the emergence of resistance. This suggests a new strategy in the war against antibiotic-therapy-resistant organisms: drug sequencing to shepherd evolution through genotype space to states from which resistance cannot emerge and by which to maximize the chance of successful therapy.
Background: The increasing rate of antibiotic resistance and slowing discovery of novel antibiotic treatments presents a growing threat to public health. Previous studies of bacterial evolutionary dynamics have shown that populations evolving on fitness landscapes follow predictable paths. In this article, we develop a general mathematical model of evolution and hypothesise that it can be used to understand, and avoid, the emergence of antibiotic resistance.

Methods and Findings:
We develop a Markov Chain model of evolution in asexually reproducing populations which is an encoding of the Strong Selection Weak Mutation model of evolution on fitness landscapes. This model associates the global properties of the fitness landscape with the algebraic properties of the Markov Chain transition matrix and allows us to derive general results on the non-commutativity and irreversibility of natural selection as well as antibiotic cycling strategies. Utilizing this formalism, we analyse 15 empirical fitness landscapes of E. coli under selection by di↵erent -lactam antibiotics and demonstrate that the emergence of resistance to a given antibiotic can be both hindered and promoted by di↵erent sequences of drug application. We show that resistance to a given antibiotic is promoted in 61.4%, 68.6% and 70.3% of possible orderings of single, pair or triple prior drug administrations, respectively. Further, we derive optimal drug application sequences with which we can probabilistically 'steer' the population through genotype space to avoid the emergence of resistance. Conclusions: Our model provides generalisable results of interest to theorists studying evolution as well as providing specific, testable predictions for experimentalists to validate our methods.
Further, these results suggest a new strategy in the war against antibiotic-therapy-resistant organisms: drug sequencing to shepherd an evolving population through genotype space to states from which resistance cannot emerge and from which we can maximize the likelihood of successful Introduction 2014]. In particular, if we understand which traits are likely to be selected for by which treatments, 23 then we may be able to avoid selecting for those traits which confer resistance. Recent insights 24 into the evolutionary process have yielded some actionable information. Specifically, Weinreich 25 et al. [2005,2006] showed that if the genome of a pathogen exhibits sign epistasis, where a given ]. These findings lead us to hypothesize that one antibiotic could be used to irreversibly steer 33 the evolution of a population of pathogens to a genotype (or combination of genotypes) which 34 is sensitive to a second antibiotic but also from which it is unlikely to acquire resistance to that antibiotic. This hypothesis was partly verified by the work of Imamovic and Sommer [2013] who 36 demonstrated that evolving E. coli to become resistant to certain antibiotics can increase sensitiv-37 ity to others. However, this work is limited in two ways. Firstly, in that the experiments do not 38 exhaustively consider all evolutionary trajectories but instead highlight only those that arose in a 39 small number of replicates and secondly, they do not consider how evolution will proceed once the 40 second drug is applied and whether resistance can then emerge.

41
In this paper we present a Markov Chain model of evolution which we use to derive general 42 results on the non-commutavity of evolution and requirements for cycling strategies. Then, using 43 previously measured landscapes for 15 -lactam anti-biotics, we illustrate that selective pressures 44 are non-commutative, and that the emergence of resistance can be both hindered and promoted 45 by di↵erent orderings of these pressures. These findings suggest new treatment strategies which 46 use rational orderings of collections of drugs to shepherd evolution through genotype space to a configuration which is sensitive to treatment, as in the work of Imamovic and Sommer [2013], but 48 also from which resistance cannot emerge. such as gene insertions, gene deletions and large structural changes to the genotype. This gives a set 56 of 2 N possible genotypes in which individuals of a given genotype, say g, can give rise to mutated 57 o↵spring which take genotypes given by one of the N mutational neighbors of g -precisely those 58 genotypes g 0 for which the Hamming distance [Hamming, 1950], Ham(g, g 0 ), from g is 1. As such, 59 our genotype space can be represented by an undirected N -dimensional (hyper-)cube graph with 60 vertices in {0, 1} N representing genotypes and edges connecting mutational neighbors (Figure 1a). 61 We define a selective pressure on our graph that drives evolution, for example through an 62 environmental change or drug application, as a fitness function This fitness function represents a genotype-phenotype map in the simplest sense -assigning to each 64 genotype a single real-valued fitness. Gillespie [1983,1984] showed that if the mutation rate u and population size M of a population satisfy Mu >> 1 >> M u 2 , and if we assume that each mutation 66 is either beneficial or deleterious, then each beneficial mutation in the population will either reach 67 fixation or become extinct before a new mutation occurs. Further, selection will be su ciently 68 strong that any deleterious mutation will become extinct with high probability and hence we may 69 assume that this always occurs. In the case that Mu 2 ⇡ 1 stochastic tunneling [Iwasa et al., for each i (see Figure 1d). Here the parameter r 0 determines the extent to which the fitness 92 increase of a mutation a↵ects its likelihood of determining the next population genotype. probability proportional to fitness increase (as in Gillespie [1983Gillespie [ , 1984Gillespie [ , 1991). This model di↵ers from the Markov model used by Sella and Hirsh [2005] to study the neutral theory of evolution as 98 we do not allow deleterious mutations to fix in the population.

99
Using this Markov Chain we can explore the possible evolutionary trajectories of a population 100 on a given fitness landscape f . We next define a collection of population row vectors µ (t) for each there are no beneficial mutations that can occur and this definition of a time step is not well defined.

111
In this case there can be no more changes to the population under the SSWM assumptions and 112 for mathematical convenience we define the probability of a local optimum population genotype 113 remaining unchanged as one in equation 3 to ensure our model is a Markov Chain. In this case the 114 step t to t + 1 can be chosen to take some fixed arbitrary time.

115
The distribution of a population at time t is related to its initial distribution, µ (0) , by 116 Since the Markov Chain is absorbing we know that there exists some k such that P k P = P k [Grin-117 stead and Snell, 1998]. Consequently, we know that the matrix exists and in fact this limit is reached after only finitely many matrix multiplications. Thus a given 119 initial population distribution µ (0) will converge to a stationary distribution µ ⇤ after a finite number 120 of steps in our model. Furthermore, if P ⇤ is known then we compute the stationary distribution In particular, provided we assume a drug is applied for su ciently long to ensure that the disease 123 population reaches evolutionary equilibrium, we can explore the e↵ects of applying multiple drugs 124 sequentially by considering the matrices P ⇤ for the associated fitness landscapes. In the following 125 discussion we make this assumption.

126
By encoding the evolutionary dynamics in a Markov Chain we can investigate the evolutionary 127 process from an algebraic perspective. In particular, as the transition matrix P encodes all of the 128 evolutionary dynamics of the associated fitness landscape f , we can explore global properties of f 129 by considering the algebraic properties of P . In the following section we present two simple, yet 130 powerful, consequences of this observation.

131
Non-Commutativity and Cycling of Natural Selection

132
We use the Markov Chain model to formally prove that for a large class of fitness landscape pairs, 133 there is non-commutativity in the evolutionary process as described by the SSWM assumptions.

134
More precisely, consider two drugs, X and Y , with corresponding fitness landscapes x and y. We

135
wish to determine what, if any, di↵erence there is between applying X followed by Y to a population 136 as opposed to applying Y followed by X to that population. If we construct the transition matrices 137 P x and P y corresponding to x and y, respectively, and take the limits P ⇤ x and P ⇤ y , then our model 138 predicts that the ordering makes no di↵erence to the final population distribution on an initial 139 population taking genotype i if, and only if, where µ i is row vector of length 2 N whose i th component is one and all of whose other components 141 are zero.

142
If we do not know the starting population genotype we can only guarantee that the order of 143 application is irrelevant when the outcome is the same regardless of the starting genotype. We 144 require for all possible length 2 N unit vectors Since these unit vectors 145 form a basis of R N this occurs precisely when Hence drug application will only commute when the corresponding limit matrices commute. In  Table 1). Of these 15 antibiotics we found no commuting pairs. We then tested 153 10 6 pairs of random fitness landscapes with varying ruggedness generated according to Kau↵man likely, then cycling strategies can be found e ciently by our model. In particular, for a given 174 starting genotype i the initial population distribution will be µ i and a sequence of drugs X 1 , . . . X k 175 with fitness landscapes x 1 , . . . x k will constitute a cycling strategy precisely when µ This criterion will be satisfied when µ i is a left 1-eigenvector of P ⇤ x 1 . . . P ⇤ x k . As such, we can find 177 cycling strategies using matrix algebra and avoid the graph-search technique used by Goulart et al.  . When Amp is given first any of the three peaks of the landscape are accessible with most resistant genotype 1111 being most likely. If Sam is given first to steer the population to its sole peak 1111, then resistance to Amp will be guaranteed when it is applied. Alternatively, if Sam is given followed by Cpr then the population evolves to the local optima genotype 0110 of the Cpr landscape. If Amp is applied to this primed population the global optimum is inaccessible.
of the fitness landscapes for three of these antibiotics, Ampicillin (Amp), Ampicillin+Sulbactam 202 (Sam) and Cefprozil (Cpr). We will use these three fitness landscapes to demonstrate the steering 203 hypothesis explicitly. In the case of a single peaked landscape, such as that for Ampicillin + 204 Sulbactam, we cannot reduce the likelihood of resistance as all evolutionary trajectories lead to the 205 global fitness optimum. It is only when a drug has a multi-peaked landscape that we may be able is necessary for the landscape to be multi-peaked.

212
In the following we take r = 0 in equation (2) and note that changing the value of r will not 213 change the accessibility of an evolutionary trajectory, hence by taking a di↵erent value of r 0 214 we will only change the result quantitatively (the probabilities may change) but not qualitatively. 215 We begin by supposing that we do not know the initial population genotype. We can model this 216 situation by taking as our prior population distribution µ = [1/2 N , . . . , 1/2 N ] which specifies that each genotype is equally likely to constitute the starting population.

218
If we apply the drug Amp to this population distribution we find that in the expected distribu-  Table 1 shows for each of the 15 antibiotics the combinations which reduce the probability 234 of evolution to a peak fitness genotype to the lowest possible when applied in order followed by 235 the final drug. We found that for 3 of the 15 drugs there exists another which steers an initial 236 population µ to a configuration which prevents evolution to the global fitness optimum of the 237 landscape entirely. This number rose to 6 when pairs of drugs applied sequentially are used to steer 238 the population and to 7 when triples applied in sequence were considered. 239 We then performed a second in silico experiment to find combinations of steering antibiotics 240 that maximize the probability that evolution proceeds to the least fit of the local optima of a final 241 antibiotic. Table 3 show the results of this experiment. We found that, excluding the single peaked for Piperacillin+Tazobactam and Ceftizoxime we found that a random steering combination will 261 increase the probability of the most highly resistance genotype occurring in more than 80% of cases, 262 suggesting that sequential multidrug treatments which use these antibiotics should proceed with 263 caution. These findings suggest that the present system of determining sequential drug orderings 264 without quantitative optimisation based guidelines could in fact be promoting drug resistance and 265 that to avoid resistance we must carefully consider the order in which drugs are applied. through matrix multiplication. We argue then that the ordering in which a collection of drugs is 273 applied can significantly impact the population that exists after the application is complete. 274 We have shown that we can find sequences of drugs that can be applied to both avoid and   Table 1: For each of the 15 antibiotics we have derived the (ordered) sets of one, two and three steering drugs which reduce the probability of evolution to the maximal fitness genotype to the lowest possible. In the case that an ordered set of steering drugs reduced the probability to 0 we have not considered ordered sets of greater length (marked as -in the table). * -as the landscape for SAM is single peaked there can be no combination of steering drugs which reduce the probability of finding the global optimum. In all experiments the initial population distribution is taken as  Table 2: For each of the 15 antibiotics we determined how many of the possible single, ordered double and ordered triple combinations of steering drugs (allowing drugs to appear multiple times in the combination and also allowing the target drug to appear as a steering drug) improved or worsened the probability of the most resistant genotype being found after the steering drugs are applied in order followed by the final target drug. In each case the initial population was given by µ = [1/2 N , . . . , 1/2 N ] and r = 0. As the SAM landscape is single peaked no combination of steering drugs will improve or worsen the outcome. As such, we have computed the overall numbers both with and without the contribution of the SAM row. due to their ine cacy.
However, a major di culty in using sequential drug treatments to steer disease populations is 285 that in order to predict the outcomes we must know the fitness landscapes of the drugs involved. De  Even in the absence of actual fitness landscapes our findings should be taken as a cautionary 301 warning for multi-drug treatments, particularly those used to treat complex diseases such as mul-302 tiple infections or cancers. In the same way that the drug ordering can be used to steer away from 303 resistance we have shown it can also be used to make resistance more likely. Our results show that 304 we may be inadvertently selecting for highly resistant disease populations through arbitrary drug 305 ordering in the same way that highly resistant disease can emerge through irresponsible drug dosing.

306
This result corroborates the findings of Pena-Miller and Beardmore [2010] that antibiotic cycling 307 strategies can vary greatly in their e cacy and be both worse and better than mixing strategies.

308
If we are to avoid resistance to our most e↵ective drugs we must carefully consider how they are 309 used together, both in combination and in sequence, with other drugs and take appropriate steps 310 to reduce the risk.

311
Two major assumptions within our modeling are that drugs are administered for su ciently 312 long that evolution can converge to a local fitness optimum and that this convergence is guaranteed 313 to occur. This assumption poses two potential problems in converting our model predictions to 314 predictions of real-world bacterial evolution. The first is that if selection is strong and mutations 315 are rare then there is a possibility of the population being driven to extinction before an adaptive 316 mutation occurs. We have chosen to ignore this possibility within our modeling as in the context 317 of treating bacterial infections this would constitute a success. The second is that the time to 318 convergence could be prohibitively long for steering to constitute a realistic treatment strategy. 319 We believe that the assumption of reasonable convergence times could be valid as adaptive walks 320 in rugged landscapes are often short [Orr, 2005]. However, it has been shown that for certain  The SSWM model is a simplification of the evolutionary process and given that non-commutativity 351 is present in this highly simplified model it is unlikely that commutativity will emerge as more com-352 plexity is introduced. It follows that the cautionary message regarding sequential drug application 353 which results from our simplified model merits serious consideration. Whether or not measuring 354 fitness landscapes provides su cient information to correctly identify optimal drug orderings in 355 vivo is a question that cannot be answered through mathematical modeling alone. It is only by     Table 3: For each of the 15 antibiotics we have derived the (ordered) sets of one, two and three steering drugs which increase the probability of evolution to the least fit optimum genotype to the highest possible. In the case that an ordered set of steering drugs increases the probability to 1 we have not considered ordered sets of greater length (marked as -in the table). * -as the landscape for SAM is single peaked there can be no combination of steering drugs which change the probability of finding the global optimum. In all experiments the initial population distribution is taken as µ = [1/2 N , . . . , 1/2 N ] and r = 0