Evolutionary non-commutativity suggests novel treatment strategies

In the Autumn of 2011 I received an email from Jacob Scott, now a good friend and better mentor, who was looking for an undergraduate to code an evolutionary simulation. Jake had just arrived in Oxford to start his DPhil in applied mathematics and by chance had dined at St Anne’s College with Peter Jeavons, then a tutor of mine, the evening before. Jake had outlined his ideas, Peter had supplied a number of email addresses, Jake sent an email and I uncharacteristically replied saying I’d give it a shot. These unlikely events would led me to where I am today — a DPhil candidate in the Oxford University Department of Computer Science. My project with Jake was a success and I was invited to speak at the 2012 meeting of the Society of Mathematical Biology in Knoxville, TN. Here I met one of Jake’s supervisors, Alexander Anderson, who invited me to visit the Department of Integrated Mathematical Oncology at the Moffitt Cancer Center and Research Institute for a workshop in December of that year. Here Dr. Anderson and I discussed one of the key issues with the work I will present in this post, issues that now form the basis of my DPhil with Dr. Anderson as one of two supervisors. Fittingly, the other is Peter Jeavons.

Jake was considering the problem of treating and avoiding drug resistance and in his short email provided his hypothesis as a single question: “Can we administer a sequence of drugs to steer the evolution of a disease population to a configuration from which resistance cannot emerge?”

In Nichol et al. (2015), we provide evidence for an affirmative answer to this question. I would like to use this post to introduce you to our result, and discuss some of the criticisms.

Jake’s hypothesis was formed from the results of two papers. In the first, Weinreich et al. (2005) showed that if the genome of a pathogen exhibits sign epistasis then there can exist inaccessible mutational trajectories. These inaccessible trajectories form the latter half Jake’s hypothesis, the inaccessible trajectory to high resistance. The second paper was by Tan et. al. (2011) which showed that adaptive mutations gained under the selective pressure of one drug are often irreversible when a second is applied. Forcing these irreversible mutations to occur was Jake’s concept of steering.

We began with the model used by Weinreich et al. and Tan et al. to derive the results which served as our motivation — the fitness landscape. Fitness landscapes have featured extensively on TheEGG and in fact long enough ago that I used some of the posts here to solidify my own understanding. For simplicity we assume that an individual has a bialleic genome in N loci, given by g \in \{0,1\}^N, and that mutation occurs as point mutations which flip a single bit of g. Thus the genotype space can be represented by an N–dimensional (hyper-)cube with genotypes as vertices and edges representing mutations. A fitness landscape is then a function which assigns to each genotype a single real, positive valued fitness:

f : \{0,1\}^N \rightarrow \mathbb{R}^{\geq 0}   .

Gillespie (1983) showed that if the population size M and mutation rate u of an evolving population satisfy Mu >> 1 >> Mu^2, and if we assume that each mutation is either beneficial or deleterious, then evolution proceeds according to Strong Selection Weak Mutation (SSWM) dynamics. Under SSWM the population can be modelled as monomorphic and inhabiting a single vertex in the genotype space. This population undergoes a stochastic up–hill walk moving at each step adjacent vertex which corresponds to a genotype of higher fitness. When no more moves can be made the population is at a local optimum of fitness.

In our work we model this stochastic walk on the hypercube graph as a Markov chain, P = [\mathbb{P} (i\rightarrow j)]_{i,j \in \{0,1\}^N}, defined by:

\mathbb{P} (i \rightarrow j) = \begin{cases}  \frac{(f(j) - f(i))^r}{\sum_{g \text{ neighbour } i \; { with } f(g) > f(i)}  (f(g) - f(i))^r} & \text{if } f(j) > f(i) \text{ and } j \text{ neighbours } i \\  1 & \text{if } i = j \text{ and } i \text{ has no fitter neighbours} \\  0 & \text{otherwise}  \end{cases}

This is based on the selection rule Artem considered at the end of his discussion of semi-smooth fitness landscapes; the parameter r \geq 0 determines the extent to which the fitness increase of a mutation affects its likelihood of fixing within the population. If r=0 each fitter neighbouring genotype is equally likely, if r \rightarrow \infty then only the fittest mutant can fix and if r=1 then we have Gillespie’s (1983) selection rule with a fitter neighbouring genotype has probability of fixing with probability proportional to the fitness increase it creates. You can see an example Markov Chain produces this way (with r = 0 in the figure below.


Using the Markov Chain encoding we can explore the possible evolutionary trajectories of a population under a given fitness landscape f. We define a collection of population row vectors \mu^{(t)} for each t \in \mathbb{N}, where \mu^{(t)} has length 2^N and k^\text{th} component which gives the probability that the population has the k\textsuperscript{th} genotype at time t (the genotypes are ordered numerically according to their binary value). The distribution of a population at time t is related to its initial distribution, \mu^{(0)}, by \mu^{(t)} = \mu^{(0)}P^t. Since the Markov Chain is absorbing we know that there exists some k such that P^kP = P^k Consequently, we know that the matrix

P^* = \lim_{t \rightarrow \infty} P^t

exists. In particular, we can compute the probability that the population evolves to each of the peaks of f by a simple matrix multiplication \mu^* = \mu^{(0)}P^*.

This Markov Chain encoding associates global properties of the fitness landscape with the algebraic properties of the transition matrix. As a simple example, suppose there two drugs F and G with fitness landscapes f and g and associated transition matrices P_f and P_g respectively. Then our model predicts that the order in which drugs are applied will make no difference to the final population probability distribution if, and only

\mu^{(0)} P^*_f P^*_g = \mu^{(0)} P^*_g

If the initial genotype distribution is unknown then we would require that this condition holds for all \mu^{(0)} and hence that P^*_f P^*_g = P^*_g P^*_f. As we know, matrix multiplication is not commutative and hence different orderings of drugs are likely to produce different results. Specifically, different orderings of the same set of drugs can have different likelihoods of resistance emerging. This result, whist mathematically simple, has profound clinical implications. Prescriptions of sequences of drugs occur frequently in the clinic, and often without any guidelines as to which orderings are preferable. As such, drug orderings are often decided arbitrarily, or from an individual clinician’s historical experience.

We sought to test what effects arbitrary drug orderings have on clinical outcomes. Using the fitness landscapes for 15 commonly used \beta–lactam antibiotics (derived by Mira et. al. 2014) we performed for each of the drugs an exhaustive in silico search of all single, double and triple combinations of steering drugs (allowing the drug itself to appear as a steerer) applied to an unknown starting population, \mu^{(0)} = [1/2^N, \dots, 1/2^N], followed by a final application of that drug. We found that in 57.3% of cases applying a single steering drug before a final drug increased the likelihood of resistance when compared to applying that final drug alone. Further, we found that 64.1% of steering pairs and 65.6% of steering triples increased the likelihood of resistance emerging when compared to application of the final drug alone. In fact, for 12/15 of the antibiotics in our study a random steering combination of one, two or three drugs was more likely to increase the probability of resistance emerging than decrease it. These results serve as a cautionary warning for clinicians prescribing sequential multidrug treatments — that we may be inadvertently selecting for highly resistant disease populations through irresponsible drug ordering in the same way that highly resistant disease can emerge through irresponsible drug dosing.

In our paper we go to test our steering hypothesis by performing for each of the 15 empirically determined drug landscapes a test of steering using combinations of one, two and three drugs to prime the population. We found that for 3 of the 15 drugs there exists another which when applied first steers an initial population \mu^{(0)} to a configuration from which the trajectory to the global fitness optimum is inaccessible. This number rose to 6 when pairs of drugs applied sequentially are used to steer the population and to 7 when triples applied in sequence were considered. These results answer Jake’s original question in the affirmative.

Criticism of the work presented here has come from two angles and primarily from two people. The first criticism, primarily raised by Artem, is a result previously outlined on this blog. Using the limit matrices P^* to simulate evolution assumes that a drug is applied for sufficiently long that the population converges to a local fitness optima. As shown by Kaznatcheev (2013), there can exist evoutionary trajectories through genotype space which have length exponential in the number of loci. As such, this assumption of convergence is a strong one. We offer the retort that the number of loci associated with antibiotic resistance is often low and further that as we get to choose the steering drugs we may be able to avoid long convergence times. Of course there is a simple way to find out — to test the predictions of our model through in vitro experiments.

The second criticism is one which is applied to many results derived from the fitness landscape model, that we are ignoring too much of the complexity at work in evolution. In particular, in fixing a fitness landscape we entirely ignore the phenotype scale and assume a static disease environment. These issues were raised by Dr. Anderson at rooftop bar in the centre of Knoxville a few days before my talk. He argued that, whilst our results are interesting, our assumptions were too strong for our results to translate to biological reality. Sandy argued that our hypothesis could be correct and that our model provided some evidence in favour of it, but that we would be unlikely to find effective steering combinations without considering more of the complexity inherent to the system. Sandy’s criticism has now become my research. In my DPhil we are aiming to build models of the genotype–phenotype map and to examine them through the algorithmic lens. I hope some of my results will find themselves onto this blog in future.

Over the course of this project I’ve learned many things. From our model I’ve learned that conventional wisdom about how to treat disease might not always be correct. This is a theme that continues in my research today where we aim to build evolutionary models and determine better treatment schedules. On a personal level I’ve learned much more. Most importantly to not be afraid of being wrong or engaging with your critics. The two biggest critics of our work so far have been Artem and Sandy. Through Artem I’m able to appear on this blog and to potentially reach a larger audience with our work. Through my conversation with Sandy in Knoxville I earned an invitation to an IMO workshop and eventually a supervisor for a DPhil in which we’re trying to build models which are a little less wrong.


Weinreich, D. M., Delaney, N. F., DePristo, M. A., & Hartl, D. L. (2006). Darwinian evolution can follow only very few mutational paths to fitter proteins. Science, 312(5770): 111-114.

Tan, L., Serene, S., Chao, H.X., & Gore, J. (2011). Hidden randomness between fitness landscapes limits reverse evolution. Physical Review Letters, 106 (19) PMID: 21668204

Gillespie, J. H. (1983). A simple stochastic gene substitution model. Theoretical population biology, 23(2), 202-215.

Mira, P. M., Crona, K., Greene, D., Meza, J. C., Sturmfels, B., & Barlow, M. (2014). Rational Design of Antibiotic Treatment Plans. arXiv preprint arXiv:1406.1564.

Nichol, D., Jeavons, P., Fletcher, A.G., Bonomo, R.A., Maini, P.K., Paul, J.L., Gatenby, R.A., Anderson, A.R.A., & Scott, J.G. (2015). Exploiting evolutionary non-commutativity to prevent the emergence of bacterial antibiotic resistance. biorXic preprint: 007542.

Kaznatcheev, A. (2013). Complexity of evolutionary equilibria in static fitness landscapes. arXiv preprint: 1308.5094.

About Dan Nichol
Presently a DPhil student in computational biology at the University of Oxford and visiting the department of Integrated Mathematical Oncology at the H. Lee Moffitt Cancer Center. Evolution, algorithms, machine learning, foundations of mathematics, good beers and video games.

11 Responses to Evolutionary non-commutativity suggests novel treatment strategies

  1. milosjohnson says:

    Cool project! I’m really interested in this idea of dynamics when you switch between landscapes, and while I agree that this isn’t super biologically applicable (because the dynamics of antibiotics and the size of the fitness landscape are more complex in reality), I think it is a really interesting direction to go in, and I’m interested to see more work like this in the future. I have one clarifying question – you mention an initial distribution of genotypes, but you are also working under SSWM – is that just meant to test a lot of initial states or are you actually talking about a genetically diverse population?

    Also, I made a website that has an interactive description of the Weinreich 2006 paper, you can check it out here if you’d like: http://milosjohnson.com/Scapes/Weinreich/

    • Dan Nichol says:

      Thank you for your comments, the initial distribution is a prior distribution representing our belief about the initial genotype. The population is monomorphic but the genotype is unknown. For this reason drugs with single peaked landscapes make good first drugs in a steering sequence – they remove uncertainty about what the genotype is.

      Your website is really nice. I’ll be releasing the code related to this work in the next couple of week if you’d like to build an interactive steering page!


      • milosjohnson says:

        That makes sense. So I went ahead and took a look at the actual paper which sorted out a lot of my questions, but now I have a comment/criticism –

        I think the condition PfPg = PgPf might be mathematically interesting, but what we care about more biologically would be PfPg = Pg or PgPfPg = Pg. The reason is that PfPg = PgPf won’t ever be true if the two landscapes have different peaks. I caught onto this because in your NK simulation you only get 0.132% that are commutative, even though you should get two single-peaked landscapes 1/25 times. Biologically, the order of two fitness landscapes if they are single peaked and the population has time to get to the peak (as in your model), isn’t really that interesting, even though it is technically not commutative (PfPg does not equal PgPf) as long as the two peaks are different (they are only the same with probability 1/2^5, which leads to the probability of commuting 1/(25*2^5) = 0.00125 – which fits your simulation if you say that sometimes you are single peaked and have the same peaks for K=1,2,3,4).

        What I’m trying to say is that I think the section about non-commutativity (and the statistics about the percentages of empirical and simulated landscape that fulfill the conditions) doesn’t quite match with the stuff later on about how the ordering can change the likelihood of getting to the most resistant state. From a biological perspective, I think this is the story, and it’s a really cool result. I just don’t think the non-commutativity of the limit matrixes really conveys that.

        I’m not sure if PfPg = Pg or PgPfPg = Pg is right, but I think that they are more relevant than PfPg = PgPf to what you’re talking about. Does that make sense?

        Glad you like the website! I’m too busy traveling right now to code something like this up, but I’d love to in the future – just flip a switch to change environments, that would be cool.

        And again, I really like this project, that’s the only reason I’m offering up my criticism unrequested online…



        • Dan Nichol says:

          Hi Milo,

          The logic behind including the non-commutativity result is to show first that order matters and then move onto showing how to exploit that fact. We’re keen to include this result because it is a simple example of the algebraic/landscape property correspondence (along with the cycling result).

          You say that PfPg = PgPf requires the peaks of the landscapes to coincide, but this isn’t quite true. For example, if f is single peaked and g double peaked but having one of those the same as f then PfPg = PgPf. However, certainly the peaks of one landscape must be a subset of the other. Do there exist more complex versions of this scenario? We don’t know. The question we’d really like to answer is an algebraic one — when do stochastic matrices commute? If we knew this then we could work that backwards to find the necessary landscape properties.

          The additional situations you talk about (PfPg = Pg or PgPfPg = Pg) are certainly of interest but our manuscript is getting a little too long to include a thorough analysis. We opted to focus on PgPf = PfPg because of its simplicity and relevance to the story of the paper. When the code is released testing the situations you suggest could be achieved very quickly.

          Again, thanks for your interest and comments on our paper. We’re in the process of looking for future extensions or applications of the model so if you have any suggestions then please get in touch!


          • milosjohnson says:

            Alright, yeah you’re right that they might just have to have one peak overlapping.

            But I gotta say I still think PfPg = Pg or PgPfPg = Pg are more relevant to the story of the paper – that we can “use rational orderings of collections of drugs to shepherd evolution through genotype space to a configuration which is sensitive to treatment… but also from which resistance cannot emerge.”
            I think my issue with emphasizing PfPg != PgPf is that you are saying “order matters” which is definitely technically true, but in a lot of cases (in all cases of single-peaked landscapes for example) all it really means is that only the last antibiotic applied matters. Again, this technically means that order matters, but not in the way that the story of the paper is concerned with, where you can potentially get to a “configuration… from which resistance cannot emerge.” This is the interesting case – where order matters to which peak you end up on in a particular environment!

            Since I was on a roll procrastinating on my actual work here I just took a look at the Mira et al 2014 paper too, and caught a little error in your manuscript – you said they used MIC as a proxy for fitness, but they actually used growth rate (Weinreich used MIC though).

            • Dan Nichol says:

              A previous version of the manuscript (before we had data) emphasised that steering is impossible, and that order doesn’t matter, when the final drug is single peaked. However, we removed this because when we got the data it turns out there’s only a single peaked landscape. Further, PfPg = Pg and PgPfPg = Pg suffer from the same problem, they will always be true for g a single peaked landscape g.

              Good catch on the MIC/Growth rate mix up!


            • milosjohnson says:

              I can’t reply to your comment for some reason so I’m replying here…

              1 – I would vote for putting that comment back in. Even though there was only one single-peaked landscape in that paper, most other papers of this type that I’ve read (the Weinreich one for example) have single peaked landscapes, so I think it’s important to mention.

              2 – I think you read my comment wrong – PfPg = PgPf will actually usually not be true for 2 single peaked landscapes (if they have different peaks), so it’s different from my two criteria in that. This means that PfPg != PgPf isn’t a very specific criteria for the relevant-to-the-story cases you are talking about (which have to involve at least one multi-peaked landscape). So the fact that PfPg = Pg and PgPfPg = Pg will always be true for a single-peaked landscape g makes them more powerful at eliminating the boring single peaked cases…

              I just wonder if you can come up with a more stringent criteria using the limit matrices that is only true (or false) when previous adaptation to one environment will result in a different outcome (different peak) during adaption in a new environment (I haven’t thought it through enough to say if PfPg = Pg or PgPfPg = Pg are really the right criteria to use, but it seems like something like that might work)

  2. Pingback: Markov Chains

  3. Pingback: An update | Theory, Evolution, and Games Group

  4. Pingback: Cataloging a year of blogging | Theory, Evolution, and Games Group

  5. Pingback: Blogging community of computational and mathematical oncologists | Theory, Evolution, and Games Group

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.