Driven Disordered Systems Approach to Biological Evolution in Changing Environments

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 trade-offs 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.

Popular Summary

Parallels between biological evolution and models of statistical physics are well known and have successfully been exploited to yield quantitative and predictive understanding of evolving populations. However, this approach has mostly been restricted to evolution in a fixed environment, the physical analogs of which often resemble physical systems in equilibrium. Generalizing this analogy to evolution in changing environments is a challenging task. Here, we investigate a class of models describing the microbial evolution to changing concentrations of an antibiotic.

We find that the adaptive behavior of the population displays a close analogy to the physics of disordered systems driven by external fields, such as sheared amorphous materials or ferromagnetic materials subject to magnetic fields. Building on ideas from the physics of driven disordered systems, we show that evolution under a varying drug concentration is generally irreversible and dependent on history, exhibiting effects of memory formation such as hysteresis and retention of partial information about past concentration changes. Moreover, this analogy allows us to derive several precise rules regarding adaptive evolution in such conditions.

While we illustrate the power of our approach by focusing on the evolution of drug resistance, the framework we develop here is more general and can be adapted to a variety of problems in evolutionary biology. Building analogies between the adaptive dynamics of biological populations and the dynamics of driven disordered systems opens the possibility of achieving a unified perspective that cuts across disciplines.

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 Ref. [21] for the antibiotic piperacillin at three different concentrations (128, 256, and $512\mu \mathrm{g}/\mathrm{ml}$). Each node is a local fitness maximum at one of these concentrations. Black and gray 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 \mathrm{g}/\mathrm{ml}$, but when the concentration is switched to $512\mu \mathrm{g}/\mathrm{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.

Figure 2

Trade-off-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 gray 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 (b). The nodes of the transition graph are the stable states which, in this simple case, comprise all genotypes. The gray 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 (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.

Figure 3

Transition graph of a realization of the TIL model (a) and its Preisach analog (b) with $L=5$ sites. The symbolic ordering sequence of this realization is given in Eq. (5). 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 gray 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.

Figure 4

TIL transition graph with $L=4$ and no secondary mutations. The ordering sequence is of the form given in Eq. (6). The transition graph is unique for this ordering sequence, and it is identical to the graph for the Preisach analog.

Figure 5

TIL transition graph with $L=5$ and nested hysteresis loops. This graph was generated from a system with trade-off 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 wild type 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.

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. (c1) and (c2) in Appendix pp3. 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 $(\ln x)/b$, where $b=⟨\ln m⟩$. The inset shows rescaled values of the mutation number in comparison to the Preisach analog (dashed lines). The rescaling is done by the mean resistance level $⟨m(x)⟩$. (b) Mean resistance level $⟨m(x)⟩$ 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 =-(a/b)$ and $a=-⟨\ln r⟩$. Angular brackets denote averages over realizations.

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 $⟨m{⟩}^{\alpha }$. The quasistatic case corresponds to the limit $\alpha \to 0$. The inset shows an enlarged 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.

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 main text and Appendix pp3. (a) 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 $⟨d_U(x)⟩$ 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$.