# A Driven Disordered Systems Approach to Biological Evolution in Changing
Environments

Suman G Das  Institute for Biological Physics, University of Cologne,
Zülpicher Straße 77, D-50937 Köln, Germany  Joachim Krug  Institute for
Biological Physics, University of Cologne, Zülpicher Straße 77, D-50937 Köln,
Germany  Muhittin Mungan  Institute for Biological Physics, University of
Cologne, Zülpicher Straße 77, D-50937 Köln, Germany  Institut für Angewandte
Mathematik, Universität Bonn, Endenicher Allee 60, D-53115 Bonn, Germany

###### Abstract

Biological evolution of a population is governed by the fitness landscape,
which is a map from genotype to fitness. However, a fitness landscape depends
on the organism’s environment, and evolution in changing environments is still
poorly understood. We study a particular model of antibiotic resistance
evolution in bacteria where the antibiotic concentration is an environmental
parameter and the fitness landscapes incorporate tradeoffs between adaptation
to low and high antibiotic concentration. With evolutionary dynamics that
follow fitness gradients, the evolution of the system under slowly changing
antibiotic concentration resembles the athermal dynamics of disordered
physical systems under external drives. Exploiting this resemblance, we show
that our model can be described as a system with interacting hysteretic
elements. As in the case of the driven disordered systems, adaptive evolution
under antibiotic concentration cycling is found to exhibit hysteresis loops
and memory formation. We derive a number of analytical results for quasistatic
concentration changes. We also perform numerical simulations to study how
these effects are modified under driving protocols in which the concentration
is changed in discrete steps. Our approach provides a general framework for
studying motifs of evolutionary dynamics in biological systems in a changing
environment.

Fitness landscapes $|$ Evolution in changing environment $|$ Driven disordered
systems $|$ Interdisciplinary evolutionary theory $|$ Statistical physics and
biological evolution

##  I Introduction

The concept of the fitness landscape of a biological species, introduced by
Sewall Wright [1], is a useful tool for understanding evolutionary processes.
According to this picture, evolving populations are driven uphill along
fitness gradients by natural selection. Mathematically, a fitness landscape is
a map which assigns fitness values to genetic sequences. In recent decades, it
has become feasible to empirically determine fitness landscapes comprising
several mutations, and a wealth of new work has illuminated various aspects of
evolutionary dynamics on different classes of fitness landscapes [2, 3, 4, 5,
6, 7, 8, 9, 10]. At the same time, the availability of empirical data has
renewed the interest in studying fitness landscapes theoretically [11, 12, 13,
14, 15, 16].

A less well-studied topic in this field is evolution in changing environments.
Fitness landscapes are a function of environment, and can change in systematic
ways as environmental parameters change. Whereas the fitness landscape
provides information about $G\times G$ (gene-gene) interactions, the
introduction of the environmental parameter furnishes information about
$G\times G\times E$ (where $E$ stands for environment) interactions, i.e.,
about how the environment modifies the gene-gene interactions [17, 18, 19,
20]. A few studies on microbial growth have measured or interpolated fitness
values as a function of environmental parameters [21, 22, 23], but systematic
theoretical work in this field is still limited.

Understanding and predicting the effect of the environment on fitness
landscapes has important practical applications. A pertinent example is the
case of antibiotic resistance in bacteria, where it has been shown that the
fitness landscape depends strongly on the antibiotic concentration [21, 22].
Uncontrolled variation in antibiotic concentration, both in clinical settings
and elsewhere [24, 25], is a cause for the rise in antibiotic resistance,
which is a major clinical challenge today. Figure 1 shows an empirical example
of the kind of processes we are interested in. The fitness values of the
genotypes (i.e. genetic sequences) in the figure were measured in [21], and
based on them, one can predict transitions between genotypes under
concentration increase (black/gray arrows) or decrease (red/orange arrows).
Notice that this small system already exhibits some interesting properties,
such as a hysteresis loop under antibiotic concentration cycling and transient
genotypes that are not part of the loop.

Figure 1: State transition graph of antibiotic resistance evolution. The nodes
depict genotypes composed of four mutations in the antibiotic resistance
enzyme TEM-50 $\beta$-lactamase. Genotypes are represented as binary strings
where a $1$ denotes the presence and $0$ the absence of a specific mutation.
The growth rates of bacteria expressing these mutant enzymes were reported in
[21] for the antibiotic piperacillin at three different concentrations ($128$
$\mu$g/ml, $256$ $\mu$g/ml and $512$ $\mu$g/ml). Each node is a local fitness
maximum at one of these concentrations. Black and grey arrows connect nodes
that would be reached under adaptive evolution when the concentration is
increased, and red and orange arrows represent the dynamics under
concentration decrease. For example, $0001$ is a local maximum at $256$ $\mu$g
ml, but when the concentration is switched to $512$ $\mu$g/ml, it is no longer
a fitness maximum. Evolution through a greedy adaptive walk (where every step
is maximally fitness increasing) leads to the new maximum $1101$. The graph
displays a hysteresis loop $0101\to 0001\to 1101\to 1100\to 0101$. The green
nodes are transient and cannot be reached under cyclic concentration changes.
Gray and orange arrows mark transitions out of transient states.

Our focus is primarily on antibiotic resistance evolution, where the
environmental dependence of the fitness landscape is governed by a tradeoff
between two phenotypes, bacterial growth rate and resistance [23]. While we
will mostly use the language of antibiotic resistance evolution in the
following, the theory developed is more generally applicable, as will become
clear from the mathematical model. Our work uses tools from statistical
physics, specifically the physics of disordered systems. Concepts and methods
from statistical physics have been used in the theory of evolution for a long
time [26, 27]. Precise quantitative analogies with evolutionary phenomena have
been found with equilibrium statistical physics [28], the theory of random
walks [29], spin glasses [30, 31, 32], and many more. Most of these, however,
focus on static fitness landscapes.

Here we investigate evolution on rugged landscapes, i.e. landscapes with a
large number of local fitness maxima, that vary with changes in an external
parameter. This setting is naturally reminiscent of the physics of driven
disordered systems, particularly in the athermal quasistatic (AQS) regime
[33], where thermal activation processes are absent or negligible. The primary
effect of the external forcing is then to alter the set of stable equilibria
or their locations. As a result, under a time-varying external forcing such
systems remain in a given equilibrium until it becomes unstable, and a fast
relaxation process leads to a new equilibrium. Despite the absence of thermal
activation processes, the resulting dynamics can nevertheless be rather
complex, exhibiting memory effects [34, 35] as well as dynamic phase
transitions, such as the jamming transition in granular materials [36], or the
yielding transition in amorphous solids [37].

In particular, we find that evolutionary genotypic change has close parallels
with systems such as cyclically sheared amorphous solids [38, 39], where a
changing environmental parameter is analogous to an external shear, and
transitions to new genotypes are similar to localized plastic events inside
the solid which can exhibit hysteresis. As was shown recently [40, 41, 35],
the AQS conditions permit a rigorous description of the dynamics of such
systems in terms of a directed state transition graph. Since the transition
graph represents the response of the system to any possible deformation
protocol, it provides a bird’s-eye view of the possible dynamics [35, 42, 43].

The main goal of this paper is to show that the driven disordered systems
approach leads to new insights into evolution in changing environments, such
as the prevalence of hysteresis, precise rules for the substitution of
mutations along fitness-increasing paths in the rare-mutation regime, and the
encoding of the evolutionary past in the genome. We establish new quantitative
results, such as the number of fitness maxima across the entire permissible
range of environmental parameters, the mean number of mutations that fix after
an instability occurs, and the extent of reversibility of the adaptive
evolution. Finally, we extend our results by going beyond the quasistatic
limit, considering the adaptive response to discrete jumps in concentration.
We find that several of our results carry over to this more realistic setting,
but there are interesting differences as well when large jumps in
concentration are involved.

Our analysis is carried out on a model of antibiotic resistance evolution.
This model is based on empirical observations on the generic properties of
dose-response curves obtained in the literature on drug resistance and
provides a principled way of describing the environment-dependence of fitness
landscapes. This introduces a new dimension into the traditional study of
fitness landscapes, which has mostly been concerned with a fixed environment.
Exploiting the analogy with disordered systems, we find that the transition
graphs describing the evolution of antibiotic resistance have a structure that
bears strong resemblance to the Preisach model [44] of hysteresis in magnets,
but in a generalized setting where the elementary units of hysteresis interact
with each other [45, 43, 46, 47, 48]. Lastly, we believe that the disordered
systems view point, in particular the state transition graph approach, is a
useful addition to the mathematical repertoire of evolutionary theory.

The paper is structured as follows. In Section II, we introduce the
mathematical model and the central concepts essential to its analysis. Section
III contains the results and is divided into seven subsections. Subsections
A-C develop the formalism of state transition graphs and construct the set of
local fitness maxima. Here we also derive several general results related to
hysteresis and memory formation under quasistatic environmental change using
this formalism. Subsection D provides statistical results regarding the number
of local fitness maxima across environments, while subsection E reports
numerical results on phenotypic reversibility, with statistical averages
performed over evolutionary trajectories with quasistatically changing
concentration. Subsection F goes beyond the quasistatic approximation by
studying dynamics under discrete changes in concentration. Subsection G
focuses on reversibility at the genotypic level and compares and contrasts it
with phenotypic reversibility. Section IV summarizes the main results and
discusses future directions for research.

##  II Model

We define a genotype $\bm{\sigma}$ as a binary string of length $L$, i.e.
$\sigma_{i}\in\\{0,1\\}$, where $i=1,2,\dots,L$ denotes the sites where
mutations can occur, and $\sigma_{i}=1$ indicates the presence of a mutation.
An equivalent and useful way of thinking about $\bm{\sigma}$ is as a set of
mutations drawn from a total of $L$ possible mutations. The genotype without
mutations, commonly referred to as the wild type, is then the empty set,
whereas the all-mutant is the set with all the $L$ mutations. We will use the
notation $\bm{\sigma}$ both as a string and as a set, and clarifications on
the notation will be provided wherever necessary.

Our focus is on the tradeoff induced landscapes (TIL) model introduced in
[23], which is defined through three key properties that are motivated by
empirical observations. 1) The fitness of each genotype $\bm{\sigma}$ is a
function of an environmental parameter $x\geq 0$, the antibiotic drug
concentration, and is described by a fitness curve of the form

$f_{\bm{\sigma}}=r_{\bm{\sigma}}w(x/m_{\bm{\sigma}}).$ (1)

The fitness curve thus has the same shape for different genotypes except for a
rescaling of the axes by the genotype-specific parameters $r_{\bm{\sigma}}$
and $m_{\bm{\sigma}}$. This is a common observation for various bacterial
strains and antibiotics [49, 23, 50]. We call $r_{\bm{\sigma}}$ the null-
fitness and $m_{\bm{\sigma}}$ the resistance of a genotype $\bm{\sigma}$,
following terminology used for bacterial dose-response curves that represent
the population growth rate as a function of drug concentration [23, 51, 49].
We choose units such that for the wild type $\bm{\sigma}=\bm{0}$,
$r_{\bm{0}}=1$ and $m_{\bm{0}}=1$, so that $f_{\bm{0}}(x)=w(x)$. Further,
$w(x)$ is a monotonic decreasing function, reflecting the decreasing fitness
of a bacterial cell with increasing drug concentration. For numerical
purposes, we will choose the widely-used Hill function form:
$w(x)=1/(1+x^{n})$, where $n$ is called the Hill exponent. Empirically
obtained dose-response curves are frequently fitted through Hill functions
[51], and the scaling property expressed in (1) shows up as a common value of
$n$ shared by various mutants of the same strain exposed to the same drug [49,
50]. 2) Every mutation comes with two parameters $r_{i}$ and $m_{i}$, and for
any genotype, $r_{\bm{\sigma}}=\exp[\sum_{i}\sigma_{i}\ln r_{i}]$ and
$m_{\bm{\sigma}}=\exp[\sum_{i}\sigma_{i}\ln m_{i}]$. Thus, the effects of
individual mutations combine in a simple multiplicative manner. This is based
on empirical observation that phenotypes of null-fitness and resistance in
fact exhibit limited or no epistasis for several microbial species and drugs
[52, 23, 53]. 3) The mutations exhibit tradeoff between adaptation to low and
high drug concentrations [25, 52, 54], i.e. $r_{i}<1$ and $m_{i}>1$. This
means that every mutation enhances the resistance, but this comes at the cost
of reduced null-fitness. The fitness curves of a specific realization of the
TIL model with $L=2$ mutations are shown in Fig 2(a).

The problem of analyzing this model has two components. First, one needs to
understand the topography of the fitness landscape, i.e. the set of local
fitness maxima and the paths that lead to the maxima, for a fixed $x$. The
second part involves questions about evolutionary dynamics between maxima
under changing drug concentrations. The first part has been addressed in
detail in [23], and we describe some of the salient features of landscape
topography here. The landscape of the TIL model is highly rugged (except at
very low and very high $x$), i.e. the number of fitness maxima is
asymptotically exponential in $L$ [23]. To describe evolutionary dynamics at
fixed $x$, it is useful to introduce the notion of a fitness graph. The nodes
of the fitness graph are the genotypes, and edges connect mutational
neighbors, i.e. genotypes that differ by a single mutation. The fitness graph
is an acyclic oriented graph, where the edges point towards increasing fitness
[16, 55]. The fitness graph depends on $x$: a fitness maximum for a certain
value of $x$ may not be a fitness maximum for another (see Fig. 2(b) for an
example). Note that the fitness graphs change only when the fitness curves of
two mutational neighbors intersect. Evolution is assumed to proceed through
adaptive walks, i.e. the entire population moves along the edges of the
fitness graph respecting their orientation [30, 56, 57, 58, 59] (see Appendix
A for further details). While this is an idealization, adaptive walks have
been found useful in the analysis of microbial evolution experiments [60, 61,
62].

Along the path taken by an adaptive walk, the fitness increases monotonically,
and such paths are therefore (evolutionarily) accessible [2, 11, 63]. The
adaptive walk terminates once a local fitness maximum is reached. In general,
there are multiple accessible paths starting from a genotype. A greedy
adaptive walk is an adaptive walk in which every step is maximally fitness
increasing [56]. A more realistic dynamics is obtained by assuming that the
probability of a transition $\bm{\sigma}\to\bm{\sigma}^{\prime}$ is
proportional to $1-e^{-2s}$ when $s>0$ and $0$ otherwise, where $s\equiv
f_{\bm{\sigma}^{\prime}}/f_{\bm{\sigma}}-1$ denotes the selection coefficient
[57, 58]. We have used this version of the adaptive walk dynamics for all
simulations, and the greedy walk has been used in the discussion of some
topological properties of specific transition graphs.

Figure 2: Tradeoff-induced fitness landscapes (TIL) model. (a) Fitness curves
for four genotypes in a TIL model with two sites ($L=2$), with parameters
$r_{1}=0.8,~{}m_{1}=1.3$ and $r_{2}=0.5,~{}m_{2}=2$. The shape of the curve is
the Hill function $w(x)=1/(1+x^{4})$. [49, 51]. The figure is divided into
five regions A-E, corresponding to different fitness graphs. A new fitness
graph occurs when the fitness curves of two mutational neighbors intersect.
The $x$-values of the intersection points are marked by the letters a-d. The
elements $x_{1}$, $x_{2}$, $\overline{x_{1}}$, $\overline{x_{2}}$ of the
ordering sequence (see main text) are indicated by solid vertical lines. (b)
Fitness graphs in the regions A-E. Concentration $x$ increases in the downward
direction. In each fitness graph, the local fitness maxima (LFMs) are marked
in red. Evolution in a fitness graph follows the oriented edges until a
fitness maximum is reached. The curved grey arrows follow the evolution of the
system under quasistatic increase of $x$ starting from the stable state $00$
at $x=0$ until the all-mutant $11$ is reached; the curved red arrows continue
the trajectory as the concentration is quasistatically decreased until $00$ is
reached again. (c) Transition graph for the two-site system is shown. This
should be distinguished from the fitness graphs in panel (b). The nodes of the
transition graph are the stable states which, in this simple case, comprise
all genotypes. The grey arrows are the U transitions, i.e. the transitions
under concentration increase, and the red arrows are the D transitions, i.e.
transitions under concentration decrease. The transition graph can be read off
from the sequence of fitness graphs in panel (b). (d) The transition between
states is shown schematically. Each horizontal level is a genotype, and the
vertical lines denote transitions. The black lines correspond to genotypes
reached under U transitions (starting from $00$ at $x=0$), while the line
traced out by the red dots indicates the genotypes reached under D transitions
(starting from $11$ at large $x$). The hysteresis loop pqrs is marked in the
figure.

We also make use of the notion of mutationally directed (or simply directed)
paths, which are paths in the fitness graph along which the number of
mutations relative to the wild type increases or decreases monotonically. The
following interesting property about the TIL landscapes at fixed $x$ was
established in [23]. It is worth discussing, since we will use it to prove
certain results in the following sections.
Directed Path Accessibility: Given a fixed concentration $x$, every
(mutationally) directed path ending at a local fitness maximum is accessible.
In other words, every local maximum $\bm{\sigma}$ is evolutionarily accessible
from the subsets (respectively supersets) of $\bm{\sigma}$ by a sequential
gain (resp. loss) of mutations; the mutations may be gained (or lost) in any
order. This property has remarkable consequences. For example, the wild type
is a subset of every genotype, and therefore can access every fitness maximum
through all directed paths. Whenever the wild type is a fitness maximum, it
must be the only fitness maximum in the landscape, since it can be accessed
from all genotypes. The same two properties also hold for the all-mutant. With
this background, we move on to the main focus of this article, which is
evolutionary dynamics under slow changes in $x$. It is here that the relation
with driven disordered systems, in particular the Preisach model, will become
apparent.

##  III Results

###  III.1 Stable states

We consider an evolutionary dynamic where the system is driven by changing the
parameter $x$, and at each value of $x$ we wait long enough so that the system
reaches a local fitness maximum through an adaptive walk. An estimate of the
necessary waiting time in terms of population-genetic parameters is provided
in Appendix A. We call a genotype a stable state if it is a local fitness
maximum (LFM) at some concentration $x$. As $x$ changes, the fitness graph is
altered by flipping the direction of one edge every time the fitness curves of
two mutational neighbors intersect (see Fig 2(a) and (b)). A stable genotype
$\bm{\sigma}$ ceases to be a LFM once its fitness curve intersects that of a
mutational neighbor, and the system transitions to a new stable state by
moving along the oriented edges of the new fitness graph.

Given a state $\bm{\sigma}$, we define the two disjoint sets
$I^{+}[\bm{\sigma}]=\\{i:\sigma_{i}=1\\},~{}I^{-}[\bm{\sigma}]=\\{j:\sigma_{j}=0\\}$.
For $i\in I^{+}[\bm{\sigma}]$, we denote by $\bm{\sigma}^{-i}$ the
configuration obtained form $\bm{\sigma}$ by setting $\sigma_{i}=0$. Likewise
for $j\in I^{-}[\bm{\sigma}]$, let $\bm{\sigma}^{+j}$ denote the configuration
obtained from $\bm{\sigma}$ by setting $\sigma_{j}=1$. Let $x_{i}$ be the
intersection point of the dose-response curves of the wild-type
$\bm{\sigma}=\bm{0}$ and the genotype with a single mutation at site $i$, i.e.
$\bm{0}^{+i}$. Hence, $x_{i}$ is the solution of the equation

$w(x)=r_{i}{w(x/m_{i})}.$ (2)

By a suitable choice of the function $w(x)$ this solution can be guaranteed to
be unique (see [23] and Appendix B). It then follows that for $\bm{\sigma}$
and $i\in I^{-}[\bm{\sigma}]$, the fitness curves of $\bm{\sigma}$ and
$\bm{\sigma}^{+i}$ intersect at $m_{\bm{\sigma}}\,x_{i}$. Likewise, for $j\in
I^{+}[\bm{\sigma}]$ the fitness curves of $\bm{\sigma}$ and $\bm{\sigma}^{-j}$
intersect at $m_{\bm{\sigma}}\,\overline{x}_{j}$, where we have defined
$\overline{x}_{j}=\frac{x_{j}}{m_{j}}$. We now see that a necessary and
sufficient condition for a genotype $\bm{\sigma}$ to be a stable state is that

$\max_{j\in I^{+}[\bm{\sigma}]}\,\overline{x}_{j}<\min_{i\in
I^{-}[\bm{\sigma}]}\,x_{i}.$ (3)

If this holds, let the index $\ell$ ($u$) correspond to the site where the
maximum (minimum) on the left (right) hand side of the inequality is attained.
The stability range of $\bm{\sigma}$ is then
$(m_{\bm{\sigma}}\overline{x}_{\ell},m_{\bm{\sigma}}{x}_{u})$. When starting
from a genotype $\bm{\sigma}$ that is a LFM at $x$, the sites $\ell$ and $u$
are the first sites that undergo a mutation when decreasing, respectively,
increasing the concentration. We will refer to sites $\ell$ and $u$ as the
least-stable sites.

At this point, we introduce a fruitful analogy with the standard Preisach
model, which is comprised of a set of non-interacting two-level systems
referred to as hysterons [44, 64]. The mutation variable
$\sigma_{i}\in\\{0,1\\}$ is analogous to the $i$-th hysteron, and its states
$0$ and $1$ correspond to the up and down states of the hysteron. The
parameter $x$ plays the role of an external magnetic field $H$ that drives the
system. In the Preisach model, each hysteron has an upper and lower threshold
$h^{+}_{i}$ and $h^{-}_{i}$, respectively. The $i$th hysteron remains in state
$0$ as long as $H<h^{+}_{i}$, and transits to $1$ otherwise. Likewise for
$H>h^{-}_{i}$ it remains in state $1$, transitioning to $0$ when this
condition does not hold. Imposing for each hysteron that
$h^{-}_{i}<h^{+}_{i}$, implies that in the range $(h^{-}_{i},h^{+}_{i})$
hysteron $i$ can be in either of its two states. The particular state chosen
is history-dependent, giving thereby rise to hysteresis. Thus a necessary and
sufficient condition for a hysteron configuration $\bm{\sigma}$ to be
attainable at some magnetic field $H$ is that

$\max_{j\in I^{+}[\bm{\sigma}]}\,h^{-}_{j}<\min_{i\in
I^{-}[\bm{\sigma}]}\,h^{+}_{i},$ (4)

which is identical to the TIL stability condition, Eq. (3). We thus define the
Preisach analogue of the TIL model as composed of $L$ hysterons, where the
upper and lower thresholds of the $i$-th hysteron are $h^{+}_{i}=x_{i}$ and
$h^{-}_{i}=\overline{x}_{i}$, respectively. Note that since $m_{i}>1$, we have
$\overline{x}_{i}<x_{i}$, so that the conditions Eq. (3) and (4) are in fact
equivalent, and we arrive at our first key result:
The stable states of the TIL model and its Preisach analogue are identical.

One immediate result following from the Preisach-TIL equivalence is the
following proximity property of stable states which holds for both models: if
$\bm{\sigma}$ is a stable state, and $\ell$ and $u$ are its least-stable
sites, then the genotypes $\bm{\sigma}^{+u}$ and $\bm{\sigma}^{-\ell}$, must
be stable states as well. The proof follows by noting that if $\bm{\sigma}$ is
stable and hence the inequality Eq. (3) holds, then by virtue of
$\overline{x}_{i}<x_{i}$, this inequality must also hold for
$\bm{\sigma}^{+u}$ and $\bm{\sigma}^{-\ell}$. The proof of the proximity
property is given in Appendix B.

For the Preisach model the proximity property implies the no-avalanche
condition [64]: when $H=h^{+}_{u}$, we have a transition from $\bm{\sigma}$ to
$\bm{\sigma}^{+u}$, such that (i) $\bm{\sigma}^{+u}$ is a stable state by peak
proximity, and moreover, (ii) $\bm{\sigma}^{+u}$ is a LFM at the field
$H=h^{+}_{u}$ that triggered the transition and hence no further state changes
occur. An analogous result holds when $H=h^{-}_{\ell}$. However, while the
proximity property (i) holds for the TIL model as well, the additional
dependence on $m_{\bm{\sigma}}$ of the stability range
$(m_{\bm{\sigma}}\overline{x}_{\ell},m_{\bm{\sigma}}{x}_{u})$ of $\bm{\sigma}$
implies that (ii) will not hold in general. As we will show next, this leads
to significant differences in the dynamical properties of the two models. In
particular, the TIL dynamics generically includes avalanches.

###  III.2 Dynamics and the transition graph

In this subsection, we discuss dynamics under quasistatically changing $x$.
While this assumption is not realistic in clinical applications, it is
amenable to analytical treatment and provides a stepping stone towards more
realistic protocols, which we will consider in subsection III.6. Crucial
aspects of the dynamics under quasistatically changing $x$ can be described by
transitions among stable states, such that the concentration is changed just
enough so that the state ceases to be a LFM. We call such transitions under
concentration increase a U-transition, and that under concentration decrease a
D-transition. Then the dynamics can be described by a transition graph (see
Fig 2(c) for an example) where the nodes are the stable states, and each node
has outgoing U and D edges corresponding to increasing and decreasing
concentrations respectively. This must be distinguished from a fitness graph,
which is defined at a fixed concentration, and where the directed edges
connecting mutational neighbors indicate the direction along which the fitness
increases.

Figure 3: Transition graph of a realization of the TIL model (left) and its
Preisach analogue (right) with $L=5$ sites. The symbolic ordering sequence of
this realization is given in Eq. 6. Each genotype is assigned an integer
label, placed within the nodes, by interpreting the genotype string as a
binary code where the leftmost digit is the least significant. The grey arrows
are U-transitions and red arrows are D transitions. The yellow nodes are the
genotypes that cannot be reached starting from the wild type. When multiple
outgoing arrows are present from a state, the solid ones correspond to greedy
walks, whereas the dashed lines represent fitness-increasing walks but with
one or more steps that are not maximally fitness-increasing.

While the TIL model and its Preisach analogue share the same set of stable
states, the dynamical properties are in general different. To illustrate this,
in Fig 3 we show a particular realization of the TIL model with $L=5$ along
with its Preisach analogue. In the Preisach model, each transition comprises
of a single switching of the least stable hysteron, which leads to a new
stable state. In the TIL model, a change at a single site need not lead to a
stable state. For example, the U-transition $15\to 27$ in the TIL model in Fig
3 involves changes at the third and fifth sites.

To understand why, first notice that under concentration changes in the U
direction, the first change in a stable state $\bm{\sigma}$ (which must
satisfy Eq. 3) is a flip $0\to 1$ at the site $u$ which has the smallest
$x_{i}$ among sites with $\sigma_{i}=0$. This flip occurs at
$x=m_{\bm{\sigma}}x_{u}$. By the proximity-property of the TIL model, the new
state $\bm{\sigma}^{+u}$ is also stable, satisfying Eq. 3. However, in order
for $\bm{\sigma}^{+u}$ to be a LFM at $x=m_{\bm{\sigma}}x_{u}$ we must require
that its lower stability threshold is less or equal to $x$. This threshold is
$m_{\bm{\sigma}^{+u}}\overline{x}_{j}=m_{\bm{\sigma}}m_{u}\overline{x}_{j}$,
where $j\in I^{+}[\bm{\sigma}^{+u}]$ and $\overline{x}_{j}>\overline{x}_{k}$
for all $k\neq j$ with $k\in I^{+}[\bm{\sigma}^{+u}]$. Therefore the new state
is a fitness maximum if and only if
$m_{u}\frac{\overline{x}_{j}}{x_{u}}={\frac{\overline{x}_{j}}{\overline{x}_{u}}\leq
1}$, in which case the dynamics terminates at $\bm{\sigma}^{+u}$. When the
condition is violated, additional secondary mutations occur until a fitness
maximum is reached. In the terminology of driven disordered systems, this
corresponds to the occurrence of an avalanche. While many detailed properties
of the secondary mutations depend on system parameters, certain features are
general and in particular do not depend on the choice of $w(x)$. In the
following we will mention some of these.

In the Preisach model, exactly $L$ U-transitions are required to get from the
wild type to the all mutant, and exactly $L$ D-transitions to go from the all-
mutant to the wild type, as is seen in Fig 3. In the TIL model, due to the
existence of secondary mutations, these numbers are generally different. The
TIL model in Fig 3 requires $6$ U-transitions to go from the wild type to the
all-mutant, and $6$ D-transitions in the reverse direction, even though $L=5$.
One important consequence of the secondary mutations is that the number of
mutations does not always increase monotonically under U or decrease
monotonically under D. The first secondary mutation is always of a
complementary kind to the initial mutation, where the changes $0\to 1$ and
$1\to 0$ are defined to be of complementary kind to each other (see Appendix B
for a proof). Further mutations may also continue to be complementary to the
initial mutation, leading to a (temporary) decrease in the number of mutations
under U or an increase under D, as shown in a typical trajectory for $L=20$
mutations in Fig 6(a). This seems counterintuitive, but arises from the state-
dependent pre-factor $m_{\bm{\sigma}}$ in the stability thresholds of stable
states.

Moreover, when secondary mutations are present, the state
$\bm{\sigma}^{\prime}$ to which a transition occurs from a state $\bm{\sigma}$
need not be unique, due to the possible presence of multiple adaptive paths.
In the TIL graph of Fig 3, the state $15$ can transition either to the state
$27$ or the state $23$ under concentration increase. It can also be shown,
using the property of direct path accessibility, that secondary mutations
cannot cause a transition to a subset or superset of $\bm{\sigma}$ (see
Appendix B), i.e. both the initial and the final state must contain at least
one mutation not contained in the other. Another related consequence of the
secondary mutations is that $\bm{\sigma}$ may transition to the same state
$\bm{\sigma}^{\prime}$ under U and D-transitions. For example, the state $5$
in the TIL graph in Fig 3 transitions to the state $3$ under both U and
D-transitions. This also appears counterintuitive from a biological
standpoint, but it can occur when the stability range of $\bm{\sigma}$ is
contained in that of $\bm{\sigma}^{\prime}$.

Figure 4: TIL transition graph with $L=4$ and no secondary mutations. The
ordering sequence is of the form given in Eq. 7. The transition graph is
unique for this ordering sequence, and is identical to the graph for the
Preisach analogue.

To understand the transition graph of the TIL model in a more systematic way,
we adopt a strategy that has been fruitful for the Preisach model [64]. We
construct a symbolic sequence $p$ that specifies the total order among all the
elements of $\\{x_{i}\\}_{1}^{L}$, $\\{\overline{x}_{j}\\}_{1}^{L}$. First,
without loss of generality, we order our indices $i$ such that
$x_{1}<x_{2}<\dots<x_{L}$. Next, it is useful to define the permutation $\rho$
of $(1,2,\ldots,L)$ that orders the $\bar{x}_{i}$ among themselves from
largest to smallest, so that
$\bar{x}_{\rho_{1}}>\bar{x}_{\rho_{2}}>\ldots>\bar{x}_{\rho_{L}}$. Since
$\overline{x}_{i}<x_{i}$ for each $i$, what remains is the specification of
the ordering relation between $\overline{x}_{i}$ and $x_{j}$ for $j\neq i$.
Given the sets $\\{x_{i}\\}_{1}^{L}$, $\\{\overline{x}_{j}\\}_{1}^{L}$ and the
ordering prescribed by $\rho$, we can describe the total ordering in terms of
a symbolic sequence $p$ of elements $\bar{i}$ and $i$ by making the
correspondence $\bar{i}\leftrightarrow\bar{x}_{i}$ and $i\leftrightarrow
x_{i}$, so that the sequence specifies the increasing order of $x_{i}$ and
$\bar{x}_{i}$. Since $\bar{x}_{i}<x_{i}$ and the permutation $\rho$ have to be
respected, the sequence has to be such that the following hold: for each $i$,
$\bar{i}$ is to the left of $i$; the subsequence of sites without overbars is
$1,2,\ldots,L$; the subsequence of sites with overbars is
$\rho_{L},\rho_{L-1},\ldots,\rho_{1}$. As an example, consider the TIL Model
in Fig 3, which has $L=5$, $\rho=(43521)$ and the ordering

$\overline{x_{1}}<\overline{x_{2}}<x_{1}<\overline{x_{5}}<\overline{x_{3}}<x_{2}<\overline{x_{4}}<x_{3}<x_{4}<x_{5}.$
(5)

The corresponding symbolic ordering sequence $p$ is then

$p=\quad\overline{1}\quad\overline{2}\quad
1\quad\overline{5}\quad\overline{3}\quad 2\quad\overline{4}\quad 3\quad 4\quad
5.$ (6)

Figure 5: TIL transition graph with $L=5$ and nested hysteresis loops. This
graph was generated from a system with tradeoff among all pairs of mutations,
i.e requiring that $r_{i}<r_{j}\iff m_{i}>m_{j}$ for all $i,j$, and
considering only transitions under greedy adaptive walks. The nesting of graph
loops implies that genotypes can partially encode past changes of
concentration. For example, starting from the wildtype $0$ at concentration
$x=0$, the genotype $19$ can be realized only by increasing first the
concentration enough in order to reach at least $23$, followed by a decrease
leading to at least $17$ (but not further than $16$) and a final increase of
concentration.

Because of the Preisach-TIL correspondence, whether a genotype/Preisach state
is stable or not, and what the least stable sites of a stable state are, can
be read off from $p$, since the condition in Eq. 3 is easy to check by
inspecting $p$ (details are given in Appendix B). In the case of the Preisach
model this implies that $p$ completely determines the transition graph [64].
While this is not the case for the TIL model, the sequence $p$ nevertheless
contains considerable information about the TIL transition graph. In
particular, this representation provides a precise condition for the existence
of secondary mutations (see Appendix B):
Secondary mutations are absent from all transitions in the TIL model if and
only if the ordering sequence is of the form

$p=\overline{1}\quad 1\quad\overline{2}\quad 2\quad\dots\quad\overline{L}\quad
L.$ (7)

In the absence of secondary mutations, the transition graph of the TIL model
becomes identical to that of its Preisach analogue. The transition graph in
this case has a simple chain structure (see Fig 4), and the number of stable
states is $L+1$, which is the lowest possible in a TIL (or Preisach) model.
Note that despite identical transition graphs in this case, some dynamical
differences are still present. Each Preisach element is hysteretic, and
therefore forward and reverse transitions between two states do not occur at
the same concentration; in the TIL model satisfying Eq. 7, however, they occur
at the same concentration, namely the one at which the dose-response curves of
the two genotypes intersect.

###  III.3 Hysteresis, reversibility and memory

The reversibility of evolution under a reversal of environmental conditions is
an important question in evolutionary biology [65, 66, 67]. In the specific
case of antibiotic resistance evolution, to what extent resistance is reversed
in a drug-free environment is a question of considerable clinical importance
[68, 69, 70]. One should note that different notions of reversion are used
here. One common definition refers to a sudden (rather than slow) change in
environment to a new state, followed by a switch back to the original state
[68, 71, 72]. In the context of our model, we consider first reversion under
quasistatic environmental changes, which would appear to be most conducive for
approximately reversible behavior. The phenomenon of reversion is naturally
linked to the notion of hysteresis under a slow and continuous change of an
external field, and indeed the Preisach model was originally proposed as a
simplified, tractable model of hysteresis in magnetic materials [44, 64]. We
then turn in section III.6 to the possibility of reversion under jump-like
concentration changes.

The TIL model also exhibits hysteresis and irreversibility. The highest degree
of reversibility is exhibited by systems with chain-like transition graphs,
such as in Fig 2(c) or Fig 4, where each transition
$\bm{\sigma}\to\bm{\sigma}^{\prime}$ is accompanied by the transition
$\bm{\sigma}^{\prime}\to\bm{\sigma}$ in the reverse direction, and there are
no states with multiple outgoing edges in either direction. This means that
under a reversal of the direction of concentration change, the same genotypes
occur in reversed sequence. However, the transitions
$\bm{\sigma}\to\bm{\sigma}^{\prime}$ and $\bm{\sigma}^{\prime}\to\bm{\sigma}$
need not occur at the same concentration. For example, in Fig 2(d), the
transition $10\to 01$ occurs at the point $x=c$ during concentration increase,
but $01\to 10$ occurs at $x=b$ during concentration decrease. On the other
hand, for systems of the type shown in Fig 4, the forward and reverse
transitions occur at the same concentration. However, such perfect
reversibility is not typical of TIL models. The degree of reversibility
depends on parameter choices. Figure 5 shows a realization of the TIL model
with $L=5$ loci and a high-degree of irreversibility, i.e. forward transitions
with no corresponding reverse transitions, such as the D-transition $16\to 1$
or the U-transition $15\to 23$.

Based on the observation of the systems described in Fig 2 and 4, we need to
distinguish between two kinds of hysteresis loops. We say that two states
$\bm{\sigma}$ and $\bm{\sigma^{\prime}}$ form a concentration loop if one can
go from $\bm{\sigma}$ to $\bm{\sigma^{\prime}}$ under quasistatic
concentration increase and from $\bm{\sigma^{\prime}}$ to $\bm{\sigma}$ under
concentration decrease, and there is some range of $x$ over which the forward
and reverse trajectories do not share the same genotype. The system in Fig 2
exhibits a concentration loop, as shown by the rectangle with corners marked
by the points $p,q,r$, and $s$ in Fig 2(d). We say that $\bm{\sigma}$ and
$\bm{\sigma}^{\prime}$ form a graph loop
$(\bm{\sigma}$,$\bm{\sigma}^{\prime})$ if one can go from $\bm{\sigma}$ to
$\bm{\sigma^{\prime}}$ under U-transitions and from $\bm{\sigma}$ to
$\bm{\sigma^{\prime}}$ under D-transitions, and if there is at least one
genotype contained in either the forward or reversed sequence of states that
is not contained in the other. A graph loop between two states implies a
concentration loop, but a concentration loop does not imply a graph loop. For
example, the case shown in Fig 2 has a concentration loop, but does not have a
graph loop, as can be seen in Fig 2(c). An example of a graph loop in the TIL
model in Fig 3(a) is the one formed by the sequence of U-transitions from $1$
to $23$ and the D-transitions leading from $23$ back to $1$. A necessary
condition for the existence of graph loops in the TIL model is the presence of
secondary mutations; for otherwise, every transition must be among mutational
neighbors, and such that the upper stability threshold of one coincides with
the lower stability threshold of the other, causing the transitions to be
reversible.

Hysteresis is also linked to the notion of memory [34, 35]. A genotype
encountered along a trajectory not only contains information about the
concentration, but also about the history of environmental change. At the
simplest level, it may contain information about whether one is on the U or D
boundary of a loop. For example, in the region between $x=b$ and $x=c$ in Fig
2(d), the state $10$ indicates that the concentration has been increasing,
while $01$ indicates that it was decreasing. But there is more information
available than this, in general. The subloops seen in the TIL graphs in Figs.
3 and 5 contain (partial) information about extreme values of $x$ reached in
previous rounds of concentration cycling. For example in the graph shown in
Fig. 5, if the dynamics started with the wild type at $x=0$ and reached the
genotype $19$ at some point, one infers that the last transition happened by
increasing the concentration to above the stability threshold of $17$; but we
also see from the transition graph that on some previous upward path the
concentration must have exceeded the upper stability threshold of $15$,
followed by some sequence of transitions that brought it to the lower
stability threshold of $16$ for the first time since this happened.

In this context, it is important to mention the phenomenon of return point
memory (RPM) [73, 40] possessed by certain systems. In our setting of adaptive
evolution, RPM implies that genotypes at which the direction of the
concentration change has been reversed can be returned to with a subsequent
reversal and hence remembered. The RPM property is universally present in the
Preisach model [34, 40, 64]. In the context of state transition graphs, one
talks about the loop-RPM property, which ensures that the system cannot escape
any loop between two states $\bm{\sigma}$ and $\bm{\sigma^{\prime}}$ without
passing through one of these states (see [40, 64] for a detailed exposition).
The RPM property implies the loop-RPM property, and is therefore possessed by
the Preisach model and can be checked for the Preisach graph in Fig 3.

Return-point memory is a mechanism by which a memory of local extremes of the
driving parameter can be retained. For example, in the TIL model in Fig 3(a)
and under greedy dynamics, starting with the wild type at $x=0$, and
increasing the concentration until state $15$ is reached, any decrease of
concentration followed by a subsequent increase will eventually lead again to
state $15$. However, if the concentration continues to increase, so that state
$23$ is reached, then a concentration decrease to say $19$ followed by an
increase will not lead to $15$ anymore. Thus the memory of $15$ as the
genotype at a local extreme event has been erased and replaced by $23$. While
we have found many realizations of the TIL model that possess the loop-RPM
property under greedy transitions, such as the transition graphs shown in
Figs. 3 and 5, it is not universally present. Additionally, the existence of
alternative fitness increasing transitions, such as the ones shown in Fig 3 by
the dashed lines, can cause a loss of this kind of memory. To see this,
consider the loop formed by the greedy U transitions leading from state $7$ to
$23$ and the greedy D transitions going from $23$ to state $7$. On the
downward trajectory from $23$, it is possible to escape the loop without going
first through $7$ by making the transition $19\to 11$.

###  III.4 Expected number of stable states

Figure 6: Average properties of evolutionary trajectories along the main loop
for $L=20$. The dose-response curve is $w(x)=1/(1+x^{2})$. The mutation
parameters $(r_{i},m_{i})$ are drawn randomly from a joint probability density
described by Eqs. (28) and (29) in Appendix C. Adaptive walks (based on
selection coefficients) were used for evolutionary dynamics. Greedy walks
produce very similar results and are not shown. A total of $10^{4}$
realizations were used to calculate averages. (a) Mean number of mutations in
the genotypes along the U and D boundaries, which describe the behavior under
increasing (U) and decreasing (D) concentrations. The dashed purple and orange
lines show a typical sample trajectory. The dashed brown line is the curve
$\frac{\ln x}{b}$, where $b=\langle\ln m\rangle$. The inset shows rescaled
values of the mutation number in comparison to the Preisach analogue (dashed
lines). The rescaling is done by the mean resistance level $\langle
m(x)\rangle$. (b) Mean resistance level $\langle m(x)\rangle$ of the genotypes
along the boundaries. The brown dashed curve is $x$. The inset shows the
fitness of the genotypes along the boundaries. The dashed green line is the
fitness of the wild type $w(x)$, and the dashed blue line is the power law
$x^{-\beta}$, where $\beta=-\frac{a}{b}$ and $a=-\langle\ln r\rangle$. Angular
brackets denote averages with respect to the distribution
$Q(\\{(r_{i},m_{i})\\})$ (see main text for further details).

While many properties of the TIL model depend only on the ordering sequences
$\rho$ and $p$, for a more detailed study of the concentration-dependent
evolutionary dynamics the fitness values of the model need to be explicitly
assigned. In the standard Preisach model, one usually considers the thresholds
to be independent random variables. Similarly for the TIL model, we assume
that $r_{i}$ and $m_{i}$ follow a joint probability distribution with density
$P(r,m)$ and the ordered pairs $(r_{i},m_{i})$ for $i=1,2,\ldots L$ are
independently and identically distributed. Their joint probability density is
then given by $Q(\\{(r_{i},m_{i})\\})=\prod_{j=1}^{L}P(r_{j},m_{j})$.

Since $x_{i}$ and $\bar{x}_{i}$ are functions of $r_{i}$ and $m_{i}$ only, the
pairs $(x_{i},\bar{x}_{i})$ for $i=1,2,\ldots L$ are independently and
identically distributed as well. Let the (marginal) cumulative distribution
function (CDF) of $x_{i}$ be $F_{x}(x_{i})$ and that of $\bar{x}_{i}$ be
$F_{\bar{x}}(\overline{x}_{i})$, and let
$P_{\bar{x}}(z)=F_{\bar{x}}^{\prime}(z)$ denote the probability density
function of $\overline{x}_{i}$. The probability that a genotype is a fitness
maximum is the probability that Eq. 3 holds. The calculation of this
probability is facilitated by the fact that $I^{+}[\bm{\sigma}]$ and
$I^{-}[\bm{\sigma}]$ are disjoint sets. One can show that in the limit of
large $L$ the average number of stable states $\langle N_{ss}\rangle$ is given
by

$\langle{N_{ss}}\rangle\simeq\sqrt{\frac{{2\pi}}{{(L-1)|G^{\prime\prime}(z_{0})|}}}{\,{\rm
e}^{G(z_{0})(L-1)}},$ (8)

where the average is taken with respect to $Q(\\{r_{i},m_{i}\\})$,
$G(z)=\ln(1+F_{\bar{x}}(z)-F_{{x}}(z))$, and the global maximum of this
function is at $z_{0}$. Thus the mean number of stable states is
asymptotically exponential in $L$, showing the highly rugged nature of these
landscapes (see Appendix C for the derivation of Eq. 8).

Figure 7: Average properties of evolutionary trajectories along the main loop
under finite driving rates $\alpha$. (a) Simulation results for the mean
number of mutations along the main hysteresis loops have been plotted for
different values of the jump sizes in concentration. The concentration has
been changed in discrete steps by factor of $\langle m\rangle^{\alpha}$. The
quasistatic case corresponds to the limit $\alpha\to 0$. The inset shows a
zoomed-in version for clarity. (b) The length of adaptive walks has been
plotted as a function of concentration for various values of $\alpha$.
Adaptive walk lengths exceeding one step imply that secondary mutations have
occurred.

###  III.5 Reversibility in the main hysteresis loop

An important dynamical question is understanding the evolutionary sequence of
genotypes as the concentration is cycled between very low and very high
values. We numerically generate trajectories starting from the wild type at
$x=0$ and increasing $x$ until the all-mutant is reached, and then decreasing
$x$ until the wild type is reached again. At each $x$, the system is evolved
through the adaptive walk (based on selection coefficients, as discussed in
Section II) until an LFM is reached. We call the trajectory comprised of the
sequences of LFMs the main hysteresis loop, and the upward and downward parts
of it the U- and D-boundary respectively. The mean number of mutations in the
genotype under quasistatic change in $x$ on both the U and D boundaries are
shown in Fig 6(a), which clearly shows hysteresis. The inset shows a
comparison with the Preisach model. Recall that in the TIL model the range of
concentrations over which a genotype $\bm{\sigma}$ is a local maximum has an
overall scale factor $m_{\bm{\sigma}}$. Therefore, in order to facilitate
comparison of the data for the TIL model and its Preisach analogue, we have
rescaled the concentration axis for the former by $\langle m(x)\rangle$, the
average scale factor $m_{\bm{\sigma}}$ of the states $\bm{\sigma}$ that are
stable at concentration $x$. From the inset of 6(a) we see that for the
Preisach model the curve for the number of mutations on the U-boundary lies
below the corresponding curve for the D-boundary. This effect can be
understood qualitatively as follows. The intersection points along the
U-boundary are governed by the distribution of $x_{i}$ and those along the
D-boundary by the distribution of $\bar{x}_{i}=x_{i}/m_{i}<x_{i}$. As a
result, the $i$-th mutation is acquired at a larger $x$ along the U-boundary
compared to where it is lost on the D boundary. More generally, for any
randomly chosen pair of mutations $i$ and $j$, $x_{j}$ tends to be higher than
$\bar{x}_{i}$ since all the $m_{i}$’s are larger than $1$. The consequence is
that the intersections along the U-boundary tend to occur at larger values of
$x$ compared to the D-boundary, making the curve for the number of mutations
along the U-boundary lower. Essentially the same effect is visible for the TIL
model when the $x$ values are rescaled by the value of $\langle m(x)\rangle$
on the boundaries. When the rescaling is not done, the U-boundary becomes
higher, as seen in the main Fig 6(a).

Figure 8: Results on secondary mutations and genotypic reversibility.
Simulations were conducted along the main hysteresis loop, and the mutation
parameters were generated randomly as described in the text and Appendix C.
(a) This figure shows the mean number of secondary mutations in an
evolutionary transition under quasistatic driving as a function of the number
of mutations in the originating background genotype. We used $10^{5}$
realizations for averaging. Black symbols correspond to the U boundary and red
to the D boundary. (b) Mean path irreversibility $\langle d_{U}(x)\rangle$ is
shown for various values of the driving rate $\alpha$. Averages were performed
over $10^{4}$ realizations. The dashed black line corresponds to the
quasistatic limit $\alpha\to 0$.

The clue to understanding this comes from Fig 6(b), which shows that the
average resistance level $\langle m(x)\rangle$ at given $x$ is lower for the
U-boundary. Since the intersection points have the pre-factor
$m_{\bm{\sigma}}$ in the TIL model, this effect tends to make the intersection
points along the U-boundary occur at lower values of $x$. For our system, this
effect is apparently strong enough to shift the curve of the number of
mutations along the U boundary above that of the D-boundary.

Generally, the changes of resistance level and mutation number have a complex
mutual dependence, and these can vary between systems depending on the dose
response curve and the parameter distribution. However, certain asymptotic
features that hold generically for stable states can be computed to leading
order. For example, as shown in Appendix C, the mean number of mutations in a
state that is stable at $x$ scales asymptotically as $\simeq\frac{\ln x}{b}$,
where the parameter $b=\langle{\ln m}\rangle$, and $m$ is the resistance of an
individual mutation. This is indicated by the brown dashed line in Fig 6(a).
At the same level of approximation, the mean level of resistance satisfies the
relation $\ln\langle m(x)\rangle\simeq\ln x$, which is shown as a dashed brown
line in Fig 6(b). The inset of Fig 6(b) displays the fitness as a function of
$x$, which is seen to decline at a much lower rate than that of the wild type,
as a consequence of the increasing level of resistance. Detailed derivations
of these results are given in Appendix C.

###  III.6 Beyond the quasistatic approximation

For many applications, a more realistic driving protocol is one where $x$ is
changed discontinuously. Experiments on antibiotic resistance evolution often
consider protocols where the concentration is changed stepwise by a constant
factor [74, 75]. To simulate this scenario, we numerically implemented a
driving protocol where concentration is sequentially increased in discrete
steps from very low to very high values and then reversed. Both the increase
and decrease occur by a factor of $\langle m\rangle^{\alpha}$ in each step,
and several values of $\alpha$ were used. The limit $\alpha\to 0$ corresponds
to quasistatic driving. The results for the mean number of mutations along the
U- and D-boundaries are plotted in Fig 7(a). The hysteresis loops are thinner,
implying a higher phenotypic reversibility. When the concentration changes by
a finite amount, the subsequent adaptive walk can move stochastically along
multiple paths involving larger numbers of mutations (see Fig 7(b)). One
therefore expects the system to reach a final distribution of phenotypes that
depends less strongly on the starting point, consistent with higher
reversibility on the phenotypic level (but not necessarily on the genotypic
level, as we will see in the next subsection).

###  III.7 Genotypic reversibility

Figure 8(a) shows simulation results with quasistatically changing $x$ for the
mean number of secondary mutations as a function of the number of background
mutations (i.e., the number of mutations in the genotype from which the
transition originates) along the main hysteresis loop. For large $L$, the
number of secondary mutations depends weakly on the number of background
mutations (unless the latter is close to $0$ or $L$). Moreover, as $L$
increases linearly, the gaps between the curves decrease, indicating a
possible asymptotic convergence towards a limiting shape (although this could
not be verified conclusively). The number of secondary mutations is seen to be
small, implying that adaptive walks at fixed $x$ are short for quasistatic
changes in $x$.

As was shown above, secondary mutations are also the source of genotypic
irreversibility under quasistatic drive in the TIL model. We now describe a
measure of genotypic irreversibility, adapted from a distance measure for
evolutionary paths introduced in [76]. Let ${\bm{\sigma}}$ be the genotype at
$x$ on the U boundary of the main loop. Then $d_{U}(x)$ is defined as the
minimum of the Hamming distance between ${\bm{\sigma}}$ and the genotypes on
the D boundary (for any concentration). The quantity $d_{D}(x)$ can be defined
in an analogous way. The quantity $\langle d_{U}(x)\rangle$ is plotted in Fig
8(b). The distance measures vanish at very low and high concentrations, which
is expected since the wild type and all-mutant are on both the U and D
boundaries of the loop. The maximum value is reached close to the
concentration at which the estimated mean number of mutations is $L/2$, as
shown by the vertical dotted line. The level of reversibility is quite high
for quasistatic driving, consistent with the low number of secondary mutations
seen in Fig 8(a). The quantity $\langle d_{D}(x)\rangle$ shows very similar
behavior (not shown).

We have also quantified the genotypic irreversibility for finite $\alpha$
using the measure $\langle d_{U}(x)\rangle$, which is shown in Fig 8(b). For
large $\alpha$, the genotypic reversibility is lower, in contrast to the
phenotypic reversibility in Fig 7(a). The walk lengths increase with
increasing $\alpha$ as shown in Fig 7(b). The longer walks lead to higher
genotypic divergence, even though the phenotypic properties of the LFMs that
are accessed are very similar, as demonstrated in Fig 7(a). Thus, even when
there is substantial hysteresis at the genotypic level, one may still observe
very similar levels of drug resistance evolving along the U and D boundaries.
For small $\alpha$, we see from Fig 8(b) that the maximum of $\langle
d_{U}\rangle$ is slightly lower than in the quasistatic limit. While the
origin of this effect is not entirely clear, it is likely related to the fact
that a small jump in concentration helps the system move beyond the
quasistatic boundary and explore a slightly larger genotypic space. The
trajectories would then be attracted preferentially towards a small number of
nearby local maxima of high fitness, leading to these states evolving
frequently along both the U and D boundaries. As $\alpha$ increases, a larger
number of maxima become accessible through the adaptive walk, leading to a
larger dispersion and reduced recurrence of the final state.

##  IV Summary and Discussion

We have investigated a class of models of bacterial evolution under changing
drug concentrations, and shown that their dynamics are closely related to that
of driven disordered systems, exhibiting dynamical phenomena such as
hysteresis and memory formation. As in the case of driven disordered systems,
a state transition graph captures the dynamics of adaptive evolution in a
changing environment. As a result, questions about evolution can be cast as
questions about the graph topology. We have shown that in such models partial
information about the changing environment is stored in evolving genotypes.
Although our analysis has focused on the relatively simple TIL model, we
should emphasize that the memory effects described here are generic features
of disordered systems.

Specifically, we have found that the transitions between genotypes in a
homogeneous population exhibit a number of generic properties which we have
described in detail. In particular, adaptive walks are found to be short,
i.e., they involve a small number of mutations, when a slow change in the drug
concentration renders a LFM unstable. In this regime, the dynamics are not
qualitatively altered by the precise nature of the adaptive walk. Moreover,
the systems generically exhibit hysteresis loops, i.e. a lack of reversibility
when the direction of concentration change is reversed. Conceptually, our work
highlights the distinction between genotypic and phenotypic hysteresis.
Genotypic reversibility is found to be low for concentration changes with
large jumps. However, the degree of phenotypic reversibility is found to be
rather high under driving protocols involving both small and large jumps in
the environmental parameter. Thus, the phenotypes are similar along both the
upward and downward boundaries of drug concentration cycling. We have obtained
asymptotic scaling approximations to phenotypes such as the number of
mutations and the fitness and resistance levels as a function of drug
concentration. These are relatively robust to modifications in the driving
protocol, making them easily accessible to empirical testing.

Under slow concentration changes, we find that the genotypic hysteresis loops
in several instances of our model exhibit perfect or near-perfect return-point
memory (RPM), i.e. the capability to return to a previously visited genotype
at which the direction of concentration change was reversed. Although not
universally present, existence of RPM is relevant to applications such as the
emergence of antibiotic resistance, since it implies that drug concentration
cycling can lead to reversal of resistance. The degree of reversibility,
including the presence of RPM, is affected by the distribution of the effects
of mutations on the null-fitness and resistance phenotypes. However, empirical
knowledge regarding the distribution of these mutational effects is still
limited. Moreover, our model assumes as an approximation that the mutations
combine non-epistatically, which rules out phenomena such as the occurrence of
compensatory mutations that reduce the cost of resistance. Further work is
needed to understand how various aspects of reversibility and memory are
affected by these factors.

From a broader perspective, our work introduces a systematic approach to
understanding how the information about a changing environment is encoded in
the genotypes of a population. A possible direction for future research is the
development of algorithms that infer features of the environmental history
from the knowledge of evolved genomes using the state transition graph.
Conversely, our approach can be used to design treatment protocols that are
optimized to avoid or slow down the evolution of drug resistance by
controlling the drug concentration or by cycling different antibiotics [77,
78].

Acknowledgments. We wish to thank two anonymous reviewers for their useful
comments on an earlier version of the manuscript. SGD and JK acknowledge
support by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) within SFB 1310 Predictability in evolution. MM was supported by
the DFG under project No. 398962893, project No. 211504053 - SFB 1060, and
under Germany’s Excellence Strategy - GZ 2047/1, project No. 390685813.

## Appendix

Here we include details of the mathematical models and computations used in
our work.

##  Appendix A Adaptive walks and changing environment

Adaptive walks describe the evolutionary dynamics in a limiting regime of
strong selection and weak mutation (SSWM) [57, 79]. Here strong selection
implies that only beneficial mutations can fix, and weak mutation implies that
no more than two genotypes are present in the population at any time. The
time-scale over which new beneficial mutations arise is $t_{\mbox{\small
mut}}=1/(N\mu)$, where $N$ is the population size and $\mu$ is the mutation
rate per genome per generation. The time-scale for fixation of a beneficial
mutation is $t_{\mbox{\small fix}}=s^{-1}\ln(N^{2}s)$ [80], where $s>0$ is the
selection coefficient of the beneficial mutation. The weak mutation condition
requires that $t_{\mbox{\small fix}}\ll t_{\mbox{\small mut}}$, which implies

$N\mu\ln(N^{2}s)\ll 1.$ (9)

In addition, we use the strong selection condition, where the fixation
probability $1-e^{-2s}$ for beneficial mutations is derived from the well-
known Kimura formula [81] in the limit

$Ns\gg 1.$ (10)

In our work, we have studied protocols in which we hold the concentration
fixed for a waiting time $t_{\mbox{\small wait}}$ until a LFM is reached.
Since this involves several evolutionary steps, a necessary condition on the
waiting time is $t_{\mbox{\small wait}}\gg t_{\mbox{fix}}$. This condition is
approximately satisfied in experimental evolution protocols where the drug
concentration is increased by a constant factor whenever the bacterial cell
density has recovered a predetermined level, signaling that a new fitness
maximum has been reached; see [75] for a recent example. The adaptive walk
approximation requires in addition the condition (9), which was not the case
in [75], but could be achieved by working at smaller population sizes.

##  Appendix B The TIL-Preisach equivalence

We derive first some results for the TIL model which will be used in the
following subsections. These results are general in the sense that they do not
depend on the shape of the dose-response function $w(x)$ as long as it
satisfies the following properties:

  * (W1)

$w(x)$ is a continuous, strictly decreasing function for $x\geq 0$.

  * (W2)

For all pairs of permissible values $(r,m)$ such that $r<1$ and $m>1$, the
curves $w(x)$ and $rw(x/m)$ intersect at precisely one point.

As noted in [23] these conditions are satisfied for a large class of dose-
response functions considered in the literature, including the exponential and
half-Gaussian functions. When considering statistical results, we will
specialize to the $n=2$ Hill-type dose-response function $w(x)=1/(1+x^{2})$.
In this case, in order for (W2) to be satisfied, the permissible pairs $(r,m)$
have to satisfy also $m^{2}r>1$.

As given in the main text, the stability condition of a state of the TIL model
is

$\max_{j\in I^{+}[\bm{\sigma}]}\,\overline{x}_{j}<\min_{i\in
I^{-}[\bm{\sigma}]}\,x_{i},$ (11)

In the following, we shall assume that we are given a stable genotype
$\bm{\sigma}$ such that

$x^{-}[\bm{\sigma}]=m_{\bm{\sigma}}\,\overline{x}_{\ell}\quad\mbox{and}\quad
x^{+}[\bm{\sigma}]=m_{\bm{\sigma}}\,x_{u},$ (12)

hold with $x^{-}[\bm{\sigma}]<x^{+}[\bm{\sigma}]$, implying that the stability
condition Eq. 11 is satisfied, and the sites $\ell$ and $u$ are those at which
the maximum, respectively minimum of the terms on the left and right hand side
of the inequality are attained. We say that sites $\ell$ and $u$ are the least
stable sites in the sense that the first mutation will occur there under
concentration decreases or increases, respectively. The stability range of
$\bm{\sigma}$ is then given as $(x^{-}[\bm{\sigma}],x^{+}[\bm{\sigma}])$.

By setting $h^{-}_{i}=\overline{x}_{i}$ and $h^{+}_{i}=x_{i}$, the condition
(11) becomes the stability condition for hysteron configuration $\bm{\sigma}$
of the Preisach model, as derived in the main text. The stability range
$(h^{-}[\bm{\sigma}],h^{+}[\bm{\sigma}])$ is then given by

$h^{-}[\bm{\sigma}]=\overline{x}_{\ell}\quad\mbox{and}\quad
h^{+}[\bm{\sigma}]=x_{u},$ (13)

which formally can be obtained from the stability range of its corresponding
TIL state, (12), by setting $m_{\bm{\sigma}}=1$. Thus while a TIL model and
its Preisach analogue have identical sets of stable states, the stability
ranges of the two are different in general.

###  B.1 The proximity-property of stable states

Given a stable TIL-Preisach state $\bm{\sigma}$, we let again $\ell$ and $u$
denote its least stable sites. The following result holds for the stability
ranges $h^{\pm}[\bm{\sigma}]$ of the Preisach model:

$\displaystyle h^{-}[\bm{\sigma}]\leq h^{-}[\bm{\sigma}^{+u}]$
$\displaystyle<h^{+}[\bm{\sigma}]<h^{+}[\bm{\sigma}^{+u}],$ (14)
$\displaystyle h^{-}[\bm{\sigma}^{-\ell}]<h^{-}[\bm{\sigma}]$
$\displaystyle<h^{+}[\bm{\sigma}^{-\ell}]\leq h^{+}[\bm{\sigma}].$ (15)

Note that these inequalities establish in particular that
$h^{-}[\bm{\sigma}^{+u}]<h^{+}[\bm{\sigma}^{+u}]$ and
$h^{-}[\bm{\sigma}^{-\ell}]<h^{+}[\bm{\sigma}^{-\ell}]$, so that
$\bm{\sigma}^{+u}$ and $\bm{\sigma}^{-\ell}$ are stable Preisach states as
well. But by the TIL-Preisach equivalence they must be stable TIL states, too.
This can also be seen by noting that the corresponding stability ranges are
obtained from those of its Preisach analogue by multiplication by
$m_{\bm{\sigma}^{+u}}$, respectively $m_{\bm{\sigma}^{-\ell}}$. This is the
proximity-property of the stable states for the TIL model and its Preisach
equivalent.

We only prove the inequalities in (14), the proof of (15) is similar. First
note that $I^{+}[\bm{\sigma}]$ is a proper subset of
$I^{+}[\bm{\sigma}^{+u}]$, and likewise $I^{-}[\bm{\sigma}^{+u}]$ is a proper
subset of $I^{-}[\bm{\sigma}]$. Since each of the two sets
$\\{h^{\pm}_{i}\\}_{i=1}^{N}$ is assumed to have distinct elements, it follows
that

$\min_{j\in I^{-}[\bm{\sigma}^{+u}]}\,h^{+}_{j}>\min_{j\in
I^{-}[\bm{\sigma}]}\,h^{+}_{j}.$ (16)

Thus $h^{+}[\bm{\sigma}^{+u}]>h^{+}[\bm{\sigma}]$ and the rightmost inequality
of (14) has been proven. Next, consider $h^{-}[\bm{\sigma}^{+u}]$ and
particularly its least stable element $k$ under field decreases. Recall that
we denoted the corresponding element for $\bm{\sigma}$ as $\ell$. Therefore,
either (i) $k=\ell$ or (ii) $k=u$. In the former case, it must have been that
$h^{-}_{u}<h^{-}_{\ell}$ and therefore
$h^{-}[\bm{\sigma}^{+u}]=h^{-}[\bm{\sigma}]=h^{-}_{\ell}$. In the letter case,
the opposite must be true, i.e. $h^{-}_{u}>h^{-}_{\ell}$, and therefore
$h^{-}[\bm{\sigma}^{+u}]=h^{-}_{u}>h^{-}[\bm{\sigma}]$. Combining these two
cases, it follows that $h^{-}[\bm{\sigma}^{+u}]\geq h^{-}[\bm{\sigma}]$,
thereby establishing the left most inequality of (14). However, since by
definition $h^{-}_{u}<h^{+}_{u}$, in both cases it must be that
$h^{-}[\bm{\sigma}^{+u}]<h^{+}[\bm{\sigma}]=h^{+}_{u}$, thereby establishing
the middle inequality of (14).

###  B.2 Construction of stable states from the symbolic order sequence $p$

Here we provide an explicit construction of the set of stable states of the
TIL model (and its Preisach equivalent) from the symbolic order sequence $p$,
which represents the ordering of the $2L$ concentrations $\overline{x}_{i}$
and $x_{j}$. Denoting the set of stable states associated with $p$ as
$\mathcal{S}_{p}$, this set can be partitioned into the subset of states
$\mathcal{S}_{p,u}$ whose least stable site under concentration increases is
$u$, with $u=1,2,\ldots,L$ (this subdivision leaves out the all-mutant
$\bm{\sigma}=\bm{1}$, which is always stable, and we assign it to the
singleton set $\mathcal{S}_{p,L+1}$).

To illustrate the construction of $\mathcal{S}_{p,u}$, consider the order
sequence (6) of the example given in the main text:

$\displaystyle\overline{x_{1}}<\overline{x_{2}}<x_{1}<\overline{x_{5}}<\overline{x_{3}}$
$\displaystyle<x_{2}<\overline{x_{4}}<x_{3}<x_{4}<x_{5},$
$\displaystyle\Leftrightarrow$ (17) $\displaystyle
p=\quad\overline{1}\quad\overline{2}\quad 1\quad$
$\displaystyle\overline{5}\quad\overline{3}\quad 2\quad\overline{4}\quad
3\quad 4\quad 5.$

Let us construct $\mathcal{S}_{p,2}$, the set of stable states with least
stable site $u=2$. All stable states must satisfy the inequality (11), and in
particular the right hand side of it must be equal to $x_{2}$. This can only
be the case, if $\sigma_{2}=0$ and $\sigma_{1}=1$, i.e. the site $1$ must
belong to $I^{+}[\bm{\sigma}]$. In terms of the order sequence $p$ this
condition is equivalent to requiring that any element $i$ without an overbar
which is located to the left of element $u$ must be assigned as
$\sigma_{i}=1$.

Likewise, in order to ensure that the left hand side of the inequality (11) is
strictly less than its right hand side, we require that any site $j$ with an
overbar to the right of element $u$ must have $\sigma_{j}=0$. In our example
this requires that $\sigma_{4}=0$. Any site $k$ left undetermined by these two
conditions can be assigned as $\sigma_{k}=0$ or $1$. In the above example
these conditions leave the sites $k=3$ and $5$ undetermined so that
$\mathcal{S}_{p,2}$ has four elements given by
$\mathcal{S}_{p,2}=\\{(10000),(10100),(10001),(10101)\\}$. Repeating the
construction for all values of $u$, the reader may verify that
$\mathcal{S}_{p}$ has $14$ states in total. These are the states shown in Fig.
3.

In the following we will make repeated use of two results whose validity is a
direct consequence of the construction of stable states given above: (i) Given
a symbolic ordering sequence $p$ and using Eq. 11, the set of all stable
genotypes $\bm{\sigma}$ can be inferred from it, and hence the possible pairs
of least stable sites $(\ell,u)$ associated with these. (ii) Given any pair of
sites $\ell$ and $u$ such that $\overline{x}_{\ell}<x_{u}$, there exists a
symbolic ordering sequence $p$ such that Eq. 12 holds and thus $\ell$ and $u$
are the least stable sites for some stable state $\bm{\sigma}$. Whenever we
assume that Eq. 12 holds, this will either imply that we are given a specific
order sequence $p$ and that with respect to $p$ the state $\bm{\sigma}$ is
stable with $(\ell,u)$ being the pair of least-stable sites, or alternatively,
we are given $(\ell,u)$ and consider the set of order sequences $p$ and stable
states $\bm{\sigma}$ compatible with this choice of least stable sites. The
particular point of view will be clear from the context.

###  B.3 Properties of secondary mutations

####  B.3.1 Secondary mutation must be complementary

We first derive some simple inequalities for the TIL model. The two properties
given below follow immediately from assumptions (W1) and (W2) made above for
$w(x)$:

$\displaystyle x$ $\displaystyle<x_{i}\quad\Leftrightarrow\quad
w(x)>r_{i}w\left(\frac{x}{m_{i}}\right),$ (18) $\displaystyle x$
$\displaystyle>x_{i}\quad\Leftrightarrow\quad
w(x)<r_{i}w\left(\frac{x}{m_{i}}\right).$ (19)

Now, let ${\bm{\sigma}}$ be a stable state, such that Eq. 12 holds. Then, we
can show from the previous results that for all $i\in
I^{-}[\bm{\sigma}]\setminus\\{u\\}$,

$\frac{f_{\bm{\sigma}^{+u,+i}}(x^{+}[\bm{\sigma}])}{f_{\bm{\sigma}^{+u}}(x^{+}[\bm{\sigma}])}<1,$
(20)

and for all $j\in I^{+}[\bm{\sigma}]\setminus\\{\ell\\}$,

$\frac{f_{\bm{\sigma}^{-\ell,-j}}(x^{-}[\bm{\sigma}])}{f_{\bm{\sigma}^{-\ell}}(x^{-}[\bm{\sigma}])}<1.$
(21)

The last two inequalities together assert that first secondary mutations which
are in the same direction as the primary mutation, are fitness-decreasing.
Thus either the first mutation leads to a LFM, and hence there will be no
further mutations, or the first secondary mutation must be complementary to
the original mutation.

####  B.3.2 Locations of complementary secondary mutations

Let ${\bm{\sigma}}$ be a stable state such that Eq. 12 holds. Then for $i\in
I^{+}[\bm{\sigma}]$,

$\displaystyle\overline{x}_{i}$
$\displaystyle>\overline{x}_{u}\quad\Leftrightarrow\quad\frac{f_{\bm{\sigma}^{+u,-i}}(x^{+}[\bm{\sigma}])}{f_{\bm{\sigma}^{+u}}(x^{+}[\bm{\sigma}])}>1.$
(22)

Therefore, subsequent to an initial mutation under concentration increase at
site $u$, fitness increasing complementary mutation sites are those sites
$i\in I^{+}[\bm{\sigma}]$ for which the symbol $\overline{i}$ is located to
the right of $\overline{u}$ in the order sequence $p$. Note in particular,
that the initial mutation site $u$ itself cannot be also the site for a
subsequent secondary mutation, as this would have implied that $\bm{\sigma}$
has a higher fitness than $\bm{\sigma}^{+u}$ at the triggering concentration.

Likewise, for $j\in I^{-}[\bm{\sigma}]$,

$\displaystyle x_{j}$
$\displaystyle<x_{\ell}\quad\Leftrightarrow\quad\frac{f_{\bm{\sigma}^{-\ell,+j}}(x^{-}[\bm{\sigma}])}{f_{\bm{\sigma}^{-\ell}}(x^{-}[\bm{\sigma}])}>1.$
(23)

Any secondary mutation following an initial mutation under concentration
decrease at site $\ell$, must be a site $j\in I^{-}[\bm{\sigma}]$ located to
the left of $\ell$ in the symbolic order sequence $p$. The statements Eq. 22
and Eq. 23 are proven by repeated application of Eq. 18, Eq. 19, and the
properties of ordering sequences $p$ that are compatible with the assumption
Eq. 12.

####  B.3.3 Secondary mutations cannot cause transitions to a subset or
superset

Assume the contrary. Then, according to the property of Directed Path
Accessibility, a path must exist where the first secondary mutation is in the
same direction as the original mutation. But this is not possible according to
the previous result.

####  B.3.4 Conditions for the absence of secondary mutations

Consider a realization of the TIL model with $L$ sites and let $\bm{\sigma}$
be a stable state satisfying Eq. 12. Assume that we are given a symbolic
ordering sequence $p$ compatible with Eq. 12. For any site $u=1,2,\ldots L$,
we will be interested in the interval of elements of $p$ that is bounded to
the left by $\overline{u}$ and to the right by $u$. Denote by
$\mathcal{I}_{u}$ the set of sites $j$ that appear in this interval without
overbars. Likewise, let $\overline{\mathcal{I}}_{u}$ be the set of sites that
appear in this interval with overbars. Our definition is such that neither of
the two sets of sites $\mathcal{I}_{u}$ and $\overline{\mathcal{I}}_{u}$
contain $u$.

Consider now transitions out of $\bm{\sigma}$ under concentration increases.
By assumption, under a concentration increase to (a value slightly above)
$x^{+}[\bm{\sigma}]$, the site $u$ will mutate first, $\sigma_{u}=0\to 1$,
leading to $\bm{\sigma}^{+u}$, and as a result, the upper limit of the
stability range of $\bm{\sigma}^{+u}$ increases to
$x^{+}[\bm{\sigma}^{+u}]>x^{+}[\bm{\sigma}]$. In order to assert the stability
of $\bm{\sigma}^{+u}$ at the concentration $x^{+}[\bm{\sigma}]$ which
triggered the mutation at $u$, we must require that

$x^{-}[\bm{\sigma}^{+u}]\leq x^{+}[\bm{\sigma}].$ (24)

If this condition is not satisfied, then
$x^{\pm}[\bm{\sigma}^{+u}]>x^{+}[\bm{\sigma}]$ and at least one secondary
mutation occurs. We thus need to find conditions under which Eq. 24 holds.

Now in terms of the ordering sequence $p$ , the site $\overline{\ell}$ must be
located to the left of $u$, as must be the site $\overline{u}$. Moreover,
since $\ell$ and $u$ have to be distinct, $\overline{\ell}$ is either to the
left or right of $\overline{u}$. In the former case we have
$\overline{x}_{\ell}<\overline{x}_{u}$, and hence

$x^{-}[\bm{\sigma}^{+u}]=m_{\bm{\sigma}^{+u}}\overline{x}_{u}=m_{\bm{\sigma}}x_{u}=x^{+}[\bm{\sigma}].$
(25)

Since Eq. 24 is satisfied, genotype $\bm{\sigma}^{+u}$ is a local fitness
maximum at this concentration and there will therefore be no secondary
mutations. Suppose next that $\overline{x}_{\ell}>\overline{x}_{u}$. In this
case

$x^{-}[\bm{\sigma}^{+u}]=m_{\bm{\sigma}^{+u}}\,\overline{x}_{\ell}=m_{\bm{\sigma}}\,\frac{\overline{x}_{\ell}}{\overline{x}_{u}}\,x_{u}>x^{+}[\bm{\sigma}].$
(26)

Therefore there will be at least one complementary secondary mutation at some
site $i\in I^{+}[\bm{\sigma}]$. Condition Eq. 22 asserts that in order for
such a mutation to be fitness increasing, $i$ must be such that
$\overline{x}_{i}>\overline{x}_{u}$. Using Eq. 11, it follows that the
stability condition of $\bm{\sigma}$, as given by Eq. 12, implies that
$\overline{x}_{i}<x_{u}$, so that the secondary mutation site must be
contained in the set $\overline{\mathcal{I}}_{u}$. Note in particular that the
site $\ell$ itself satisfies these conditions and hence is a possible
candidate for the first secondary mutation.

Combining all of the above results, under increasing concentration, a
secondary mutation will occur, if and only if the set
$\overline{\mathcal{I}}_{u}$ is non-empty, and the state $\bm{\sigma}$ is such
that, for some $j\in\overline{\mathcal{I}}_{u}$ we have $\sigma_{j}=1$.
Conversely, a secondary mutation will not occur if and only if one of the
following two conditions holds:

  * (U1)

The set $\overline{\mathcal{I}}_{u}$ is empty.

  * (U2)

The set $\overline{\mathcal{I}}_{u}$ is non-empty, and the state $\bm{\sigma}$
is such that, for each $j\in\overline{\mathcal{I}}_{k}$ we have
$\sigma_{j}=0$.

In a similar manner, one can show that under decreasing concentration a
secondary mutation will not occur, if and only if one of the following two
conditions holds:

  * (D1)

The set $\mathcal{I}_{\ell}$ is empty.

  * (D2)

The set $\mathcal{I}_{\ell}$ is non-empty, and the state $\bm{\sigma}$ is such
that, for each $j\in\mathcal{I}_{\ell}$ we have $\sigma_{j}=1$.

Observe now that in order for secondary mutations to be absent from all
transitions in a TIL model, the sets $\mathcal{I}_{k}$ and
$\overline{\mathcal{I}}_{k}$ have to be empty for each $k=1,2,\ldots,L$, since
otherwise there will exist stable states for which conditions (U2) or (D2) can
be made not to hold. The only ordering sequence for which both of these sets
are empty for each $k$ is the sequence

$p=\overline{1}\quad 1\quad\overline{2}\quad 2\quad\dots\quad\overline{L}\quad
L.$ (27)

##  Appendix C Statistical Results

In the following subsections, were derive the statistical results discussed in
the main body of the paper.

###  C.1 Probability density function used in the numerics

We assume that the dose-response curve is of Hill-type with $n=2$. In order to
satisfy the requirement (W2) for the dose-response function, the parameters
$(r_{j},m_{j})$ must be chosen such that $m_{j}^{2}r_{j}>1$ for each
$j=1,2,\ldots,L$. We further assume that the pairs $(r_{j},m_{j})$ are
independently and identically distributed, so that their joint density is
given by $Q(\\{r_{i},m_{i}\\})=\prod_{j=1}^{L}P(r_{j},m_{j})$. We write
$P(r_{j},m_{j})=P_{1}(r_{j})P_{2}(m_{j}|r_{j})$. We chose

$\displaystyle P_{1}(r)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2}{\pi}}\frac{e^{-\frac{(\ln r)^{2}}{2}}}{r}$ (28)
$\displaystyle P_{2}(m|r)$ $\displaystyle=$
$\displaystyle\mathcal{N}~{}\frac{e^{-\frac{(\ln
m)^{2}}{2}}}{m}\Theta\big{(}m-\frac{1}{\sqrt{r}}\big{)},$ (29)

where $\Theta(\cdot)$ is the Heaviside step function, and $\mathcal{N}$ is the
appropriate normalization constant. This choice is for easy of implementation.
A similar model was used in [23].

###  C.2 Asymptotic number of stable states

Consider a genotype $\bm{\sigma}$ with $n$ mutations, i.e.
$\sum_{i}\sigma_{i}=n$. The number of such genotypes is $\binom{L}{n}$ and
such a $\bm{\sigma}$ is a stable state if Eq. 11 holds. Since $x_{i}$ and
$\bar{x}_{i}$ are independent for distinct sites, the probability density that
the left hand side of Eq. 11 is less than $z$ and the right hand side is
greater than $z$ is $\frac{d}{dz}F^{L-n}_{\bar{x}}(z)[1-F_{x}(z)]^{n}$. Then
the mean number of stable states is

$\displaystyle\langle N_{ss}\rangle$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{L}{\left(\begin{array}[]{l}L\\\
n\end{array}\right)}\int
dz~{}\left[\frac{d}{dz}F^{L-n}_{\bar{x}}(z)\right]\left[1-F_{x}(z)\right]^{n}$
(32) $\displaystyle=$ $\displaystyle L\int
dz\,{\left[1-F_{x}(z)+F_{\bar{x}}(z)\right]^{L-1}\,{F^{\prime}_{\bar{x}}(z)}},$
(33)

from which the result in the main text follows using a saddle point
approximation for $L$ large.

###  C.3 Asymptotic approximation for number of mutations

The fitness $f$ of a genotype ${\bm{\sigma}}$ can be expressed as

$\ln f=\sum_{i}\sigma_{i}\ln
r_{i}-\ln\left(1+x^{2}\,e^{-2\,\sum_{i}\sigma_{i}\ln m_{i}}\right)$ (34)

The number of mutations in the genotype is $n=\sum_{i}\sigma_{i}$. A simple
heuristic that produces good approximations for the mean of various quantities
at large $L$ is a s follows: we consider the fitness of a genotype to be a
function of $x$ and $n$ only, and replace the parameters associated with the
mutations by suitable averages. Thus, we write Eq. 34 as

$\ln f(n)\simeq-na-\ln\left(1+x^{2}\,e^{-2nb}\right),$ (35)

where $a=-\langle\ln r\rangle$ and $b=\langle\ln m\rangle$. For any given $x$,
one can now maximize Eq. 35 with respect to $n$, yielding an approximation to
the mean mutation number at $x$ for stable maxima. Taking the derivative of
the above with respect to $n$ and setting it to zero produces the equation:

$\frac{2bx^{2}e^{-2nb}}{1+x^{2}e^{-2nb}}=a.$ (36)

The solution to this is

$n=\frac{\ln x}{b}+\frac{1}{2b}\ln\left(\frac{2b}{a}-1\right).$ (37)

For large $x$ and therefore large $n$, the leading order is

$n\simeq\frac{\ln x}{b}.$ (38)

This estimate works well when $L$ is large and $1\ll\frac{\ln(x)}{\langle\ln
m\rangle}\ll L$.

## References

  * Wright [1932] S. Wright, The roles of mutation, inbreeding, crossbreeding, and selection in evolution, in  _Proc. Sixth. Int. Cong. Genet._ , Vol. 1 (1932) pp. 356–366.
  * Weinreich  _et al._ [2006] D. M. Weinreich, N. F. Delaney, M. A. DePristo, and D. L. Hartl, Darwinian evolution can follow only very few mutational paths to fitter proteins, Science 312, 111 (2006).
  * Lozovsky  _et al._ [2009] E. R. Lozovsky, T. Chookajorn, K. M. Brown, M. Imwong, P. J. Shaw, S. Kamchonwongpaisan, D. E. Neafsey, D. M. Weinreich, and D. L. Hartl, Stepwise acquisition of pyrimethamine resistance in the malaria parasite, Proc. Natl. Acad. Sci. USA 106, 12025 (2009).
  * Schenk  _et al._ [2013] M. F. Schenk, I. G. Szendro, M. L. M. Salverda, J. Krug, and J. A. G. M. de Visser, Patterns of epistasis between beneficial mutations in an antibiotic resistance gene, Molecular biology and evolution 30, 1779 (2013).
  * de Visser and Krug [2014] J. A. G. M. de Visser and J. Krug, Empirical fitness landscapes and the predictability of evolution, Nature Reviews Genetics 15, 480 (2014).
  * Palmer  _et al._ [2015] A. C. Palmer, E. Toprak, M. Baym, S. Kim, A. Veres, S. Bershtein, and R. Kishony, Delayed commitment to evolutionary fate in antibiotic resistance fitness landscapes, Nature communications 6, 1 (2015).
  * Bank  _et al._ [2016] C. Bank, S. Matuszewski, R. T. Hietpas, and J. D. Jensen, On the (un)predictability of a large intragenic fitness landscape, Proc. Natl. Acad. Sci. USA 113, 14085 (2016).
  * Domingo  _et al._ [2018] J. Domingo, G. Diss, and B. Lehner, Pairwise and higher-order genetic interactions during the evolution of a tRNA, Nature 558, 117 (2018).
  * Fragata  _et al._ [2019] I. Fragata, A. Blanckaert, M. A. D. Louro, D. A. Liberles, and C. Bank, Evolution in the light of fitness landscape theory, Trends in ecology & evolution 34, 69 (2019).
  * Pokusaeva  _et al._ [2019] V. O. Pokusaeva, D. R. Usmanova, E. V. Putintseva, L. Espinar, K. S. Sarkisyan, A. S. Mishin, N. S. Bogatyreva, D. N. Ivankov, A. V. Akopyan, S. Y. Avvakumov, I. S. Povolotskaya, G. J. Filion, L. B. Carey, and F. A. Kondrashov, An experimental assay of the interactions of amino acids from orthologous sequences shaping a complex fitness landscape, PLoS Genetics e1008079, 15 (2019).
  * Franke  _et al._ [2011] J. Franke, A. Klözer, J. A. G. M. de Visser, and J. Krug, Evolutionary accessibility of mutational pathways, PLoS Computational Biology 7, e1002134 (2011).
  * Szendro  _et al._ [2013] I. G. Szendro, M. F. Schenk, J. Franke, J. Krug, and J. A. G. M. de Visser, Quantitative analyses of empirical fitness landscapes, Journal of Statistical Mechanics: Theory and Experiment 2013, P01005 (2013).
  * Neidhart  _et al._ [2014] J. Neidhart, I. G. Szendro, and J. Krug, Adaptation in tunably rugged fitness landscapes: the rough mount fuji model, Genetics 198, 699 (2014).
  * Ferretti  _et al._ [2016] L. Ferretti, B. Schmiegelt, D. Weinreich, A. Yamauchi, Y. Kobayashi, F. Tajima, and G. Achaz, Measuring epistasis in fitness landscapes: the correlation of fitness effects of mutations, Journal of theoretical biology 396, 132 (2016).
  * Blanquart and Bataillon [2016] F. Blanquart and T. Bataillon, Epistasis and the structure of fitness landscapes: are experimental fitness landscapes compatible with fisher’s geometric model?, Genetics 203, 847 (2016).
  * Crona  _et al._ [2017] K. Crona, A. Gavryushkin, D. Greene, and N. Beerenwinkel, Inferring genetic interactions from comparative fitness data, eLife 6, e28629 (2017).
  * Murugan  _et al._ [2021] A. Murugan, K. Husain, M. J. Rust, C. Hepler, J. Bass, J. M. Pietsch, P. S. Swain, S. G. Jena, J. E. Toettcher, A. K. Chakraborty, _et al._ , Roadmap on biology in time varying environments, Physical Biology 18, 041502 (2021).
  * de Vos  _et al._ [2018] M. G. J. de Vos, S. E. Schoustra, and J. A. G. M. de Visser, Ecology dictates evolution? About the importance of genetic and ecological constraints in adaptation, Europhysics Letters 122, 58002 (2018).
  * Gorter  _et al._ [2018] F. A. Gorter, M. G. M. Aarts, B. J. Zwaan, and J. A. G. M. de Visser, Local fitness landscapes predict yeast evolutionary dynamics in directionally changing environments, Genetics 208, 307 (2018).
  * Anderson  _et al._ [2021] D. W. Anderson, F. Baier, G. Yang, and N. Tokuriki, The adaptive landscape of a metallo-enzyme is shaped by environment-dependent epistasis, Nature Communications 12, 3867 (2021).
  * Mira  _et al._ [2015a] P. M. Mira, J. C. Meza, A. Nandipati, and M. Barlow, Adaptive landscapes of resistance genes change as antibiotic concentrations change, Molecular biology and evolution 32, 2707 (2015a).
  * Ogbunugafor  _et al._ [2016] C. B. Ogbunugafor, C. S. Wylie, I. Diakite, D. M. Weinreich, and D. L. Hartl, Adaptive landscape by environment interactions dictate evolutionary dynamics in models of drug resistance, PLoS computational biology 12, e1004710 (2016).
  * Das  _et al._ [2020] S. G. Das, S. O. Direito, B. Waclaw, R. J. Allen, and J. Krug, Predictable properties of fitness landscapes induced by adaptational tradeoffs, Elife 9, e55155 (2020).
  * Kolpin  _et al._ [2004] D. W. Kolpin, M. Skopec, M. T. Meyer, E. T. Furlong, and S. D. Zaugg, Urban contribution of pharmaceuticals and other organic wastewater contaminants to streams during differing flow conditions, Science of the Total Environment 328, 119 (2004).
  * Andersson and Hughes [2014] D. I. Andersson and D. Hughes, Microbiological effects of sublethal levels of antibiotics, Nature Reviews Microbiology 12, 465 (2014).
  * Fisher [1958] R. A. Fisher,  _The genetical theory of natural selection_ (Dover, 1958).
  * Koonin [2011] E. V. Koonin,  _The logic of chance: the nature and origin of biological evolution_ (FT press, 2011).
  * Sella and Hirsh [2005] G. Sella and A. E. Hirsh, The application of statistical physics to evolutionary biology, Proceedings of the National Academy of Sciences 102, 9541 (2005).
  * Manhart and Morozov [2014] M. Manhart and A. V. Morozov, Statistical physics of evolutionary trajectories on fitness landscapes, in  _First-passage phenomena and their applications_ (World Scientific, 2014) pp. 416–446.
  * Hwang  _et al._ [2018] S. Hwang, B. Schmiegelt, L. Ferretti, and J. Krug, Universality classes of interaction structures for NK fitness landscapes, Journal of Statistical Physics 172, 226 (2018).
  * Stein [1992] D. L. Stein,  _Spin glasses and biology_ (World Scientific, 1992).
  * Franz  _et al._ [1993] S. Franz, L. Peliti, and M. Sellitto, An evolutionary version of the random energy model, J. Phys. A: Math. Gen. 26, L1195 (1993).
  * Maloney and Lemaître [2006] C. E. Maloney and A. Lemaître, Amorphous systems in athermal, quasistatic shear, Phys. Rev. E 74, 016118 (2006).
  * Keim  _et al._ [2019] N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Memory formation in matter, Rev. Mod. Phys. 91, 035002 (2019).
  * Mungan  _et al._ [2019] M. Mungan, S. Sastry, K. Dahmen, and I. Regev, Networks and hierarchies: How amorphous materials learn to remember, Phys. Rev. Lett. 123, 178002 (2019).
  * Behringer and Chakraborty [2018] R. P. Behringer and B. Chakraborty, The physics of jamming for granular materials: a review, Reports on Progress in Physics 82, 012601 (2018).
  * Bonn  _et al._ [2017] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys. 89, 035005 (2017).
  * Regev  _et al._ [2013] I. Regev, T. Lookman, and C. Reichhardt, Onset of irreversibility and chaos in amorphous solids under periodic shear, Phys. Rev. E 88, 062401 (2013).
  * Fiocco  _et al._ [2013] D. Fiocco, G. Foffi, and S. Sastry, Oscillatory athermal quasistatic deformation of a model glass, Phys. Rev. E 88, 020301 (2013).
  * Mungan and Terzi [2019] M. Mungan and M. M. Terzi, The structure of state transition graphs in hysteresis models with return point memory: I. general theory, Ann. Henri Poincaré 20, 2819 (2019).
  * Mungan and Witten [2019] M. Mungan and T. A. Witten, Cyclic annealing as an iterated random map, Phys. Rev. E 99, 052132 (2019).
  * Regev  _et al._ [2021] I. Regev, I. Attia, K. Dahmen, S. Sastry, and M. Mungan, Topology of the energy landscape of sheared amorphous solids and the irreversibility transition, Phys. Rev. E 103, 062614 (2021).
  * Keim and Paulsen [2021] N. C. Keim and J. D. Paulsen, Multiperiodic orbits from interacting soft spots in cyclically sheared amorphous solids, Science Advances 7, eabg7685 (2021).
  * Preisach [1935] F. Preisach, Über die magnetische Nachwirkung, Zeitschrift für Physik 94, 277 (1935).
  * Hovorka and Friedman [2005] O. Hovorka and G. Friedman, Onset of reptations and critical hysteretic behavior in disordered systems, Journal of magnetism and magnetic materials 290, 449 (2005).
  * Lindeman and Nagel [2021] C. W. Lindeman and S. R. Nagel, Multiple memory formation in glassy landscapes, Science Advances 7, eabg7133 (2021).
  * van Hecke [2021] M. van Hecke, Profusion of transition pathways for interacting hysterons, Phys. Rev. E 104, 054608 (2021).
  * Bense and van Hecke [2021] H. Bense and M. van Hecke, Complex pathways and memory in compressed corrugated sheets, Proceedings of the National Academy of Sciences 118, e2111436118 (2021).
  * Chevereau  _et al._ [2015] G. Chevereau, M. Dravecká, T. Batur, A. Guvenek, D. H. Ayhan, E. Toprak, and T. Bollenbach, Quantifying the determinants of evolutionary dynamics leading to drug resistance, PLOS Biology 13, e1002299 (2015).
  * Lukačišinová  _et al._ [2020] M. Lukačišinová, B. Fernando, and T. Bollenbach, Highly parallel lab evolution reveals that epistasis can curb the evolution of antibiotic resistance, Nature communications 11, 1 (2020).
  * Regoes  _et al._ [2004] R. R. Regoes, C. Wiuff, R. M. Zappala, K. N. Garner, F. Baquero, and B. R. Levin, Pharmacodynamic functions: a multiparameter approach to the design of antibiotic treatment regimens, Antimicrobial agents and chemotherapy 48, 3670 (2004).
  * Marcusson  _et al._ [2009] L. L. Marcusson, N. Frimodt-Møller, and D. Hughes, Interplay in the selection of fluoroquinolone resistance and bacterial fitness, PLoS pathogens 5, e1000541 (2009).
  * Knopp and Andersson [2018] M. Knopp and D. I. Andersson, Predictable phenotypes of antibiotic resistance mutations, mBio 9, e00770 (2018).
  * Melnyk  _et al._ [2015] A. H. Melnyk, A. Wong, and R. Kassen, The fitness costs of antibiotic resistance mutations, Evolutionary applications 8, 273 (2015).
  * Crona  _et al._ [2013] K. Crona, D. Greene, and M. Barlow, The peaks and geometry of fitness landscapes, Journal of Theoretical Biology 318, 1 (2013).
  * Kauffman and Levin [1987] S. Kauffman and S. Levin, Towards a general theory of adaptive walks on rugged landscapes, Journal of theoretical Biology 128, 11 (1987).
  * Orr [2002] H. A. Orr, The population genetics of adaptation: the adaptation of dna sequences, Evolution 56, 1317 (2002).
  * Seetharaman and Jain [2014] S. Seetharaman and K. Jain, Adaptive walks and distribution of beneficial fitness effects, Evolution 68, 965 (2014).
  * Agarwala and Fisher [2019] A. Agarwala and D. S. Fisher, Adaptive walks on high-dimensional fitness landscapes and seascapes with distance-dependent statistics, Theor. Pop. Biol. 130, 13 (2019).
  * Rokyta  _et al._ [2005] D. R. Rokyta, P. Joyce, S. B. Caudle, and H. A. Wichman, An empirical test of the mutational landscape model of adaptation using a single-stranded dna virus, Nature genetics 37, 441 (2005).
  * Rokyta  _et al._ [2009] D. R. Rokyta, Z. Abdo, and H. A. Wichman, The genetics of adaptation for eight microvirid bacteriophages, Journal of Molecular Evolution 69, 229 (2009).
  * Schoustra  _et al._ [2009] S. E. Schoustra, T. Bataillon, D. R. Gifford, and R. Kassen, The properties of adaptive walks in evolving populations of fungus, PLoS biology 7, e1000250 (2009).
  * Weinreich  _et al._ [2005] D. M. Weinreich, R. A. Watson, and L. Chao, Perspective: Sign epistasis and genetic costraint on evolutionary trajectories, Evolution 59, 1165 (2005).
  * Terzi and Mungan [2020] M. M. Terzi and M. Mungan, State transition graph of the Preisach model and the role of return-point memory, Phys. Rev. E 102, 012122 (2020).
  * Teotónio and Rose [2001] H. Teotónio and M. R. Rose, Perspective: Reverse evolution, Evolution 55, 653 (2001).
  * Bridgham  _et al._ [2009] J. T. Bridgham, E. A. Ortlund, and J. W. Thornton, An epistatic ratchet constrains the direction of glucocorticoid receptor evolution, Nature 461, 515 (2009).
  * Kaltenbach  _et al._ [2015] M. Kaltenbach, C. J. Jackson, E. C. Campbell, F. Hollfelder, and N. Tokuriki, Reverse evolution leads to genotypic incompatibility despite functional and active site convergence, eLife 4, e06492 (2015).
  * Andersson and Hughes [2010] D. I. Andersson and D. Hughes, Antibiotic resistance and its cost: is it possible to reverse resistance?, Nature Reviews Microbiology 8, 260 (2010).
  * Allen  _et al._ [2017] R. C. Allen, J. Engelstädter, S. Bonhoeffer, B. A. McDonald, and A. R. Hall, Reversing resistance: different routes and common themes across pathogens, Proceedings of the Royal Society B: Biological Sciences 284, 20171619 (2017).
  * Durão  _et al._ [2018] P. Durão, R. Balbontín, and I. Gordo, Evolutionary mechanisms shaping the maintenance of antibiotic resistance, Trends in microbiology 26, 677 (2018).
  * Dunai  _et al._ [2019] A. Dunai, R. Spohn, Z. Farkas, V. Lázár, Á. Györkei, G. Apjok, G. Boross, B. Szappanos, G. Grézal, A. Faragó, _et al._ , Rapid decline of bacterial drug-resistance in an antibiotic-free environment through phenotypic reversion, Elife 8, e47088 (2019).
  * Zur Wiesch  _et al._ [2010] P. S. Zur Wiesch, J. Engelstädter, and S. Bonhoeffer, Compensation of fitness costs and reversibility of antibiotic resistance mutations, Antimicrobial agents and chemotherapy 54, 2085 (2010).
  * Sethna  _et al._ [1993] J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Hysteresis and hierarchies: Dynamics of disorder-driven first-order phase transformations, Phys. Rev. Lett. 70, 3347 (1993).
  * Baym  _et al._ [2016] M. Baym, T. D. Lieberman, E. D. Kelsic, R. Chait, R. Gross, I. Yelin, and R. Kishony, Spatiotemporal microbial evolution on antibiotic landscapes, Science 353, 1147 (2016).
  * Schenk  _et al._ [2022] M. F. Schenk, M. P. Zwart, S. Hwang, P. Ruelens, E. Severing, J. Krug, and J. De Visser, Population size mediates the contribution of high-rate and large-benefit mutations to parallel evolution, Nature Ecology & Evolution 6, 439 (2022).
  * Lobkovsky and Koonin [2012] A. E. Lobkovsky and E. V. Koonin, Replaying the tape of life: quantification of the predictability of evolution, Frontiers in Genetics 3, 246 (2012).
  * Mira  _et al._ [2015b] P. M. Mira, K. Crona, D. Greene, J. C. Meza, B. Sturmfels, and M. Barlow, Rational design of antibiotic treatment plans: A treatment strategy for managing evolution and reversing resistance, PloS ONE 10, e0122283 (2015b).
  * Nichol  _et al._ [2015] D. Nichol, P. Jeavons, A. G. Fletcher, R. A. Bonomo, P. K. Maini, J. L. Paul, R. A. Gatenby, A. R. Anderson, and J. G. Scott, Steering evolution with sequential therapy to prevent the emergence of bacterial antibiotic resistance, PloS Comp. Biol. 11, e1004493 (2015).
  * Gillespie [1984] J. H. Gillespie, Molecular evolution over the mutational landscape, Evolution 38, 1116 (1984).
  * Gerrish and Lenski [1998] P. J. Gerrish and R. E. Lenski, The fate of competing beneficial mutations in an asexual population, Genetica 102, 127 (1998).
  * Kimura [1962] M. Kimura, On the probability of fixation of mutant genes in a population, Genetics 47, 713 (1962).