Evolution in the weak-mutation limit: Stasis periods punctuated by fast transitions between saddle points on the fitness landscape

Yuri Bakhtin, Mikhail I. Katsnelson, Yuri I. Wolf, Eugene V. Koonin

Contributed by Eugene V. Koonin, December 16, 2020 (sent for review July 24, 2020; reviewed by Sergey Gavrilets and Alexey S. Kondrashov)

Author contributions: Y.B., M.I.K., Y.I.W., and E.V.K. designed research; Y.B. performed research; Y.B., M.I.K., Y.I.W., and E.V.K. analyzed data; and Y.B. and E.V.K. wrote the paper.

Reviewers: S.G., University of Tennessee; and A.S.K., University of Michigan.

https://doi.org/10.1073/pnas.2015665118, Volume: 118, Issue: 4, Pages: 1-9

Article Type: research-article Article History

Publisher: National Academy of Sciences

- Facebook
- Twitter
- Linkedin
- Whatsapp
Altmetric

Table of Contents

Results
Discussion
Data Availability.

Abstract

The gradual character of evolution is a key feature of the Darwinian worldview. However, macroevolutionary events are often thought to occur in a nongradualist manner, in a regime known as punctuated equilibrium, whereby extended periods of evolutionary stasis are punctuated by rapid transitions between states. Here we analyze a simple mathematical model of population evolution on fitness landscapes and show that, for a large population in the weak-mutation limit, the process of adaptive evolution consists of extended periods of stasis, which the population spends around saddle points on the landscape, interrupted by rapid transitions to new saddle points when a beneficial mutation is fixed. Thus, phenomenologically, the default regime of biological evolution seems to closely resemble punctuated equilibrium.

A mathematical analysis of the evolution of a large population under the weak-mutation limit shows that such a population would spend most of the time in stasis in the vicinity of saddle points on the fitness landscape. The periods of stasis are punctuated by fast transitions, in lnN_e/s time (N_e, effective population size; s, selection coefficient of a mutation), when a new beneficial mutation is fixed in the evolving population, which accordingly moves to a different saddle, or on much rarer occasions from a saddle to a local peak. Phenomenologically, this mode of evolution of a large population resembles punctuated equilibrium (PE) whereby phenotypic changes occur in rapid bursts that are separated by much longer intervals of stasis during which mutations accumulate but the phenotype does not change substantially. Theoretically, PE has been linked to self-organized criticality (SOC), a model in which the size of “avalanches” in an evolving system is power-law-distributed, resulting in increasing rarity of major events. Here we show, however, that a PE-like evolutionary regime is the default for a very simple model of an evolving population that does not rely on SOC or any other special conditions.

Keywords

fitness landscapes, low mutation limit, saddle points, punctuated equilibrium, self-organized criticality

License: CC BY

Bakhtin,Katsnelson,Wolf,and Koonin: Evolution in the weak-mutation limit: Stasis periods punctuated by fast transitions between saddle points on the fitness landscape

Phyletic gradualism, that is, evolution occurring via a succession of mutations with infinitesimally small fitness effects, is a central tenet of Darwin’s theory (1). However, the validity of gradualism has been questioned already by Darwin’s early, fervent adept, T. H. Huxley (2), and subsequently many nongradualist ideas and models have been proposed, to account, primarily, for macroevolution. Thus, Goldschmidt (in)famously championed the hypothesis of “hopeful monsters,” macromutations that would be deleterious in a stable environment but might give their carriers a chance for survival after a major environmental change (3). Arguably, the strongest motivation behind nongradualist evolution concepts was the notorious paucity of intermediate forms in the fossil record. It is typical in paleontology that a species persists without any major change for millions of years but then is abruptly replaced by a new one. The massive body of such observations prompted Simpson, one of the founding fathers of the modern synthesis of evolutionary biology, to develop the concept of quantum evolution (4), according to which species, and especially higher taxa, emerged abruptly, in “quantum leaps,” when an evolving population rapidly moved to a new “adaptive zone,” or, using the language of mathematical population genetics, a new peak on the fitness landscape. Simpson proposed that the quantum evolution mechanism involved fixation of unusual allele combinations in a small population by genetic drift, followed by selection driving the population to the new peak.

The idea of quantum evolution received a more systematic development in the concept of punctuated equilibrium (PE) proposed by Eldredge and Gould (5 6 7–8). The abrupt appearance of species in the fossil record prompted Eldredge and Gould to postulate that evolving populations of any species spend most of the time in the state of stasis, in which no major phenotypic changes occur (9, 10). The long intervals of stasis are punctuated by short periods of rapid evolution during which speciation occurs, and the previous dominant species is replaced by a new one. Gould and Eldredge emphasized that PE was not equivalent to the “hopeful monsters” idea, in that no macromutation or saltation was proposed to occur, but rather a major acceleration of evolution via rapid succession of “regular” mutations that resulted in the appearance of instantaneous speciation, on a geological scale. The occurrence of PE is traditionally explained via the combined effect of genetic drift during population bottlenecks and changes in the fitness landscape that can be triggered by environmental factors (11).

PE has been explicitly linked to the physical theory of self-organized criticality (SOC). SOC, a concept developed by Bak (12), is an intrinsic property of dynamical systems with multiple degrees of freedom and strong nonlinearity. Such systems experience serial “avalanches” separated in time by intervals of stability (the avalanche metaphor comes from Bak’s depiction of SOC on the toy example of a sand pile, on which additional sand is poured, but generally denotes major changes in a system). A distinctive feature of the critical dynamics under the SOC concept is self-similar (power law) scaling of avalanche sizes (12 13 14 15 16 17–18). The close analogy between SOC and PE was noticed and explored by Bak and colleagues, the originators of the SOC concept, who developed models directly inspired by evolving biological systems and intended to describe their behavior (12, 15, 16, 18). In particular, the popular Bak–Sneppen model (15) explores how ecological connections between organisms (physical proximity in the model space) drive coevolution of the entire community. Extinction of the organisms with the lowest fitness disrupts the local environments and results in concomitant extinction of their closest neighbors. It has been shown that, after a short burn-in, such systems self-organize in a critical quasi-equilibrium interrupted by avalanches of extinction, with the power law distribution of avalanche sizes.

A distinct but related view of macroevolution is encapsulated in the concept of major transitions in evolution developed by Szathmáry and Maynard Smith (19 20–21). Under this concept, major evolutionary transitions, such as, for instance, emergence of multicellular organisms, involve emergence of new levels of selection (new Darwinian individuals), in this case selection affecting ensembles of multiple cells rather than individual cells. These evolutionary transitions resemble phase transitions in physics (22) and appear to occur rapidly, compared to the intervals of evolution within the same level of selection. The concept of evolutionary transitions can be generalized to apply to the emergence of any complex feature including those that do not amount to a major change in the level of biological organization (23).

We sought to assess the validity of evolutionary gradualism by mathematically investigating the simplest conceivable model of population evolution on a rugged fitness landscape (24). We show that, under the basic assumptions of a large population size and low mutation rate (weak-mutation limit), an evolving population spends most of the time in stasis, that is, percolating through near-neutral mutational networks around saddle points on the landscape. The intervals of stasis are punctuated by rapid transitions to new saddle points after fixation of beneficial mutations. Thus, contrary to the general perception of the weak-mutation limit as a paragon of gradualism (25), we find that the default evolutionary mode in this regime resembles PE while not requiring SOC or any other special conditions.

Results

Agent-Based Model of Competitive Exclusion.

We consider a well-mixed population of a large constant size N consisting of individuals, each with a specific genotype. To avoid dealing with the overwhelming complexity of the space of all genotypes, we work with a coarse-grained model that groups similar genotypes into “types.” The genotypes within the same type are considered to be homogeneous and densely connected by the mutation network. The only homogeneity assumption we need to make is that, within each type, the variations in fitness and available transitions to other classes due to mutations are negligible. We also assume that sizes of different types are comparable. The set of all types is denoted by $T$ .

The evolution of a population within the model involves reproduction and mutation. Reproduction of individuals occurs under the Moran model widely used in population genetics, that is, with rates proportional to their fitness, and is accompanied by removal of random individuals to keep N constant (26). Mutations are modeled by transitions in a mutational network E that might involve one or more elementary genetic mutations. The individual mutation rate λ is assumed to be low compared to the reproduction rates. The evolutionary regime depends on 1) the geometry of the graph ( $T$ , E), 2) the fitness function f, 3) the values of parameters N and λ, and 4) the initial configuration.

Let us now describe our basic model in more detail. We assume that the population size is a large number $N$ , constant in time. The set $T$ of all possible types is finite or countable. It can be viewed as a graph with adjacency matrix ${(E_{i j})}_{i, j \in T}$ . Two distinct types $i, j$ are connected by an edge if they differ by a mutation (at the scale of the model, a mutation is assumed to occur instantaneously and without intermediate steps). In that case, we set $E_{i j} = 1$ . Otherwise, $E_{i j} = 0$ .

Each type $i \in T$ is assigned a fitness value $f_{i} > 0$ which is identified with the reproduction rate. The numbers $f_{i}$ are assumed to be distinct and of the order of $1$ (more precisely, bounded), so essentially time is measured in reproductions. It is convenient to work with relative sizes $y_{i}$ of type populations (fractions) with respect to the total population size $N$ . We denote by $Δ$ the space of sequences ${(y_{i})}_{i \in T}$ such that $y_{i} \geq 0$ for all $i$ and $\sum_{i \in T} y_{i} = 1$ . Denoting the fraction of individuals of type $i \in T$ present in the population at time $t \in ℝ$ by $x_{i} (t)$ (taking values $0, N^{- 1}, 2 N^{- 1}, \dots$ ), we define random evolution of the vector ${(x_{i} (t))}_{i \in T} \in Δ$ as a continuous time pure jump $Δ$ -valued Markov process, by specifying the transition rates. A single individual of type $i \in T$ produces new individuals of the same type $i$ at the rate $f_{i}$ . Each reproduction is accompanied by removal of one individual that is randomly and uniformly chosen from the entire population. Thus, the total rate of reproduction of individuals of type $i$ is $N x_{i} f_{i}$ . Given that an individual of type $i$ is reproducing, the probability that the child individual will replace an individual of type $j$ is $x_{j}$ . Thus, the total rate of simultaneous change $x_{i} \to x_{i} + N^{- 1}$ and $x_{j} \to x_{j} - N^{- 1}$ is $N f_{i} x_{i} x_{j}$ . Let us now introduce mutations. We will assume that mutation rates are much lower than the reproduction rates. To model this, we introduce a small parameter $λ > 0$ . The rate of replacement of an individual of type $i \in I (x)$ , where

I (x) = \{i \in T : x_{i} > 0\}, x \in Δ,

by an individual of type

j

is given by

λ E_{i j} \in {0, λ}

. The total rate of such transitions occurring in a population is

N λ E_{i j} x_{i}

In what follows, we derive the PE-like evolutionary regime from several reasonable assumptions on the geometry of the graph, the fitness function, population size, mutation rates, and the initial state. Our results can be viewed as similar to those in previous work (27 28–29), where more mathematically sophisticated models were considered. However, our simple model allows for a more transparent analysis that is conducive to biological implications and we use it here to tie the PE concept to noisy dynamics near heteroclinic networks (30, 31) and emphasize the importance of saddle points on the landscape for the evolutionary process.

Evolution without Mutations in the Infinite Population Size Limit.

In this section, we examine the case where, in an infinite population, $λ = 0$ , that is, there are no mutations, and approximate the dynamics of our stochastic model by that of a deterministic ordinary differential equation (ODE):

{\dot{x}}_{i} = b_{i} (x), i \in T,

[1]

with the right-hand side given by

b_{i} (x) = x_{i} (f_{i} - \bar{f} (x)),

where

\bar{f} (x) = \underset{}{\sum_{j \in T}} f_{j} x_{j}

is the average fitness for the population state

x

. The system () is a well-known competitive exclusion system [see, e.g., (2.15) to (2.16) of ref. ] restricted to nonzero components of

x

. emerges due to the averaging effect and can be viewed as a law of large numbers for our model.

To state the results, we need to introduce some notations and definitions. We denote $I = I (x (0))$ for brevity and note that, given the absence of mutations, our stochastic model and ODE (1) are defined on the simplex $Δ_{I} = \{x \in ℝ_{+}^{I} : \sum_{i \in I} x_{i} = 1\}$ . This simplex is the convex hull of its vertices $e^{(i)}$ , $i \in I$ , corresponding to pure states where only one type is present:

e_{k}^{(i)} = {\begin{array}{l} 1, & i = k, \\ 0, & i \neq k . \end{array}

One of these vertices plays a special role. Let

i^{*}

be the type with maximum fitness

f^{*}

(within I), that is,

f^{*} = f_{i^{*}} = \max_{i \in I} f_{i}

. We will see that

e^{(i^{*})}

is an attractor for both deterministic dynamical system defined by (1) and for our stochastic model.

Solutions of the system (1) admit a concise analytic form (33):

x_{i} (t) = \frac{x_{i} (0) e^{f_{i} t}}{\sum_{j \in I} x_{j} (0) e^{f_{j} t}}, i \in I .

Since we assumed that the fitness values for each type are distinct, this formula implies straightforwardly that

x (t)

converges to

e^{(i^{*})}

exponentially fast, with exponential rate given by the selection coefficient

s = f^{*} - f^{* *},

[2]

where

f^{* *}

is the second-largest fitness value in

I .

On the one hand, our main approximation result given below implies that this attraction to the maximum fitness state

e^{(i^{*})}

under typical conditions holds also for the stochastic agent-based system; on the other hand, it quantifies possible deviations from this behavior.

For the approximation result, we need to define the discrepancy

D (t) = x (t) - Φ^{t} x (0),

[3]

where

x (t)

is the Markov process without mutations and for any

y

Φ^{t} y

is the solution of ODE () with the initial condition

y

, at time

t

. We are going to estimate the maximum discrepancy up to time

t

, that is,

D^{*} (t) = \sup_{u \in [0, t]} ∥ D (u) ∥

, where

∥ \cdot ∥

is the

L^{1}

norm in

ℝ^{I}

defined by

∥ x ∥ = \sum_{i \in I} | x_{i} | .

[4]

We assume that the number of types

| I |

is small compared to the population size, more precisely, there is

μ < 1 / 2

such that

| I | \leq N^{μ} .

[5]

Because this model does not include mutations, if a type

i

becomes extinct at time

t_{0}

, that is,

x_{i} (t_{0}) = 0

, then

x_{i} (t) = 0

for all

t \geq t_{0}

. We denote the event on which no type

i \in I

becomes extinct before time

t

B_{t} = {I (x (u)) = I

for all

u \in [0, t]}

. Events from a sequence

{(A_{N})}_{N \in ℕ}

are stretch-exponentially unlikely (SE-unlikely) if for some

C, γ > 0

P (A_{N}) \leq C e^{- N^{γ}}, N \in ℕ .

This is fast decay in

N,

just short of being truly exponentially fast. We are now ready to state our main result for the system without mutations and to examine the meaning of each of its parts.

Theorem 1.

Assume 5. Then:

1)There are constants $c, β > 0$ such that events $B_{c \ln N} \cap \{D^{*} (c \ln N) > N^{- β}\}$ are SE-unlikely.

2)Let $β$ be defined in Part 1 of the theorem. Then, for any $δ < β$ , there is a constant $C > 0$ such that, conditioned on the nonextinction of type $i^{*}$ , and up to an SE-unlikely event, $| x (C \ln N) - e^{(i^{*})} | \leq N^{- δ}$ .

3)There are constants $C', α > 0$ such that, if $| x (0) - e^{(i^{*})} | \leq N^{- δ}$ , then

P \{x (C' \ln N) = e^{(i^{*})}\} > 1 - N^{- α} .

4)There is a number $p > 0$ that does not depend on $N$ , such that the probability of nonextinction of type $i^{*}$ is bounded below by $p$ for all initial conditions $x (0)$ satisfying $x_{i^{*}} (0) > 0$ .

5)For any $δ \in (0,1)$ , if $x_{i^{*}} (0) > N^{- δ}$ , then, extinction of type $i^{*}$ is SE-unlikely.

Part 1 of the theorem shows that, up to time $c \ln N$ , if no type gets extinct, the stochastic process $x (t)$ follows the deterministic trajectory $Φ^{t} x (0)$ very closely, deviating from it at most by $N^{- β}$ . This happens with a probability very close to 1, exceptions being SE-unlikely.

Part 2 shows that, if type $i^{*}$ does not die out, then, with a high probability, by time $C \ln N$ , it will dominate the population whereas all other types will be almost extinct. The proof of the Theorem shows that $C = 1 / s$ , where s is the exponential rate of attraction to $e^{(i^{*})}$ given by the the selection coefficient and defined by (2).

Part 3 means that, after realization of the scenario described in Part 2 and an additional logarithmic time, $i^{*}$ will be the only surviving type.

Part 1 is conditioned on the nonextinction of any type, whereas Part 2 is conditioned on the nonextinction of type $i^{*}$ . If any type $i$ dies out, Part 1 still applies to the continuation of the process on the simplex $Δ_{I \ {i}}$ of a lower dimension. By contrast, for Part 2 to be meaningful, we need to provide a bound on the nonextinction of $i^{*}$ . This is done in Parts 4 and 5.

Part 4 states that there is a positive probability (independent of the population size) that the progeny of even a single individual of type $i^{*}$ will drive out all other types.

Part 5 states that, once the fraction of the individuals of type $i^{*}$ reaches a (small) threshold $N^{- δ}$ , then, it is almost certain that $i^{*}$ will dominate the population.

To summarize these results, the chance of extinction for the fittest type is nonnegligible only when there are very few individuals of this type, that is, when the initial state involves a recent mutation that produced a single individual of this type. Once the number of individuals reaches a certain modest threshold, the typical, effectively deterministic, behavior will follow the trajectory of Eq. 1 closely, eventually reaching the pure state of fixation where only individuals of type i* are present. The proof of Theorem 1 is given in the end of this section. Now, we turn to the analysis of the dynamics generated by ODE (1).

Heteroclinic Network.

The points $e^{(k)}$ are hyperbolic critical points (saddles) of various indices (the index of a saddle is the number of negative eigenvalues of the linearization of the vector field at that saddle). We can find these eigenvalues and associated eigenvectors explicitly. The linearization $(\partial_{j} b_{i} (e^{(k)}))$ of $b$ at $e^{(k)}$ is

\begin{aligned} \partial_{k} b_{k} (e^{(k)}) = & - f_{k}, \\ \partial_{i} b_{k} (e^{(k)}) = & - f_{i}, & i \neq k, \\ \partial_{i} b_{i} (e^{(k)}) = & f_{i} - f_{k}, & i \neq k, \\ \partial_{j} b_{i} (e^{(k)}) = & 0, & j \neq i, i \neq k . \end{aligned}

Therefore, for each

i \in I

such that

i \neq k

, there is an eigenvalue

f_{i} - f_{k}

(\partial_{j} b_{i} (e^{(k)}))

with an eigenvector

e^{(i)} - e^{(k)}

pointing along the simplex edge connecting

e^{(k)}

and

e^{(i)}

. These eigenvectors span the simplex

Δ_{I}

, so the additional eigenvalue

- f_{k}

with eigenvector

e^{(k)}

that is transversal to

Δ_{I}

can be ignored. To demonstrate explicitly that the vertex

e^{(k)}

is a saddle, we note that the eigendirections given by

e^{(i)} - e^{(k)}

are stable or unstable, depending on the sign of the associated eigenvalue, that is, on whether

f_{i} < f_{k}

f_{i} > f_{k}

. Moreover, there is a heteroclinic connection (a trajectory connecting two distinct saddle points) between

e^{(i)}

and

e^{(k)}

. This trajectory coincides with the simplex edge between

e^{(i)}

and

e^{(k)}

and corresponds to the presence of exactly two types

i, k

. The dynamics along this trajectory is described by the logistic equation

{\dot{x}}_{i} = (f_{i} - f_{k}) x_{i} (1 - x_{i})

(see Fig. 1 for the phase portrait).

$The phase portrait of the dynamical system (1). Four types, 1, 2, 3, and 4, are shown such that f1< f2<f3 <f4. The dynamics is defined on the simplex Δ{1,2,3,4} with vertices e(1),e(2),e(3),e(4), corresponding to pure states where the population consists entirely of individuals of one type. These vertices are critical points of the vector field b. The edges of the simplex are heteroclinic orbits connecting these critical points to each other. Several other orbits are also plotted as arrows. The vertex e(4) attracts every initial condition with nonzero fraction of individuals of the fittest type i∗=4.$

Fig. 1.

The phase portrait of the dynamical system (1). Four types, 1, 2, 3, and 4, are shown such that $f_{1}$ < $f_{2}$ < $f_{3}$ < $f_{4}$ . The dynamics is defined on the simplex $Δ_{{1,2,3,4}}$ with vertices $e^{(1)}, e^{(2)}, e^{(3)}, e^{(4)},$ corresponding to pure states where the population consists entirely of individuals of one type. These vertices are critical points of the vector field b. The edges of the simplex are heteroclinic orbits connecting these critical points to each other. Several other orbits are also plotted as arrows. The vertex $e^{(4)}$ attracts every initial condition with nonzero fraction of individuals of the fittest type $i^{*} = 4$ .

The key feature of this dynamics is a heteroclinic network formed by trajectories connecting saddle points to one another. The vertex $e^{(i^{*})}$ is a sink (a saddle with the maximal index) if considered in $Δ_{I}$ , but it can also be viewed as a saddle in simplices of higher dimensions based on coordinates (types) that include those with higher fitness than $f^{*}$ . The types with higher fitness will appear if we include mutations into the model.

Evolutionary Process with Mutations.

We now consider the full process with positive but small mutation rate $λ$ and recall that, for each type $i \in I (x)$ , the rate of mutation to type $j$ is given by $λ E_{i j}$ . We consider here only relatively late stages of evolution that are preceded by extensive evolutionary optimization so that the overwhelming majority of the mutations are either deleterious or at best neutral. More precisely, we assume that there is a constant M such that for each $i \in I (x)$ the total number of available fitness-increasing (beneficial) mutations, that is, vertices $j \in T$ such that $E_{i j} = 1$ and $f_{j} > f^{*}$ , is bounded by $M$ . Our first assumption on the magnitude of $λ$ is that

r (N) = λ N \ln N / s ≪ 1,

where

s

is the selection coefficient defined in 2.

According to Theorem 1 and the accompanying discussion, if the evolutionary process is conditioned on the survival of type $i^{*}$ , then, typically, it takes $C \ln N$ time for the process $x_{i^{*}}$ (t) to reach 1 (fixation), where $C = 1 / s$ . The probability of a beneficial mutation during this time interval is bounded by

1 - e^{- M N λ C \ln N} = 1 - e^{- M r (N)} \leq M r (N) .

[6]

Thus, the population is unlikely to produce a new beneficial mutation before it reaches the state of fixation where type

i^{*}

is the only surviving one. Once a new beneficial mutation occurs and, accordingly a new best-fit type emerges, it either gets extinct quickly or gets fixed in the population, in time of the order

\ln N / s

. The trajectory, driven by differential reproduction of random mutations, closely follows the heteroclinic connection, that is, the line connecting two vertices of the simplex

Δ

. The entire process can be described as follows: There is a moment when

i^{*}

is the only type present, after which it takes time of the order (kλN)⁻¹ to produce a new beneficial mutation, where k is the number of beneficial mutations that are available from

i^{*}

. Then, it takes time of the order

\ln N / s

for this fittest type to take over the entire population, after which the process repeats. At least for mutations with a relatively large

s

, on the order of 0.01 to 0.1, this time is much shorter than the time the population spends in the vicinity of a saddle point, “waiting” for the beneficial mutation to occur.

Now consider deleterious mutations. There are N individuals, and each produces a suboptimal (lower fitness) type with the rate λL, where L is the number of available deleterious mutations. Using the Poisson distribution, we obtain that, by the time $t$ , it is highly unlikely to produce more than $t$ NλL new suboptimal individuals. If $t = \ln N / s$ , then, this number is λLN lnN $/ s$ , so requiring

λ L \ln N / s ≪ 1,

[7]

we obtain

λ L N \ln N / s ≪ N

, that is, over the travel time between saddles the emerging individuals with deleterious mutations constitute an asymptotically negligible fraction of the entire population. Thus, the trajectory

x (t)

will be altered only by a term converging to

0

N \to \infty

The resulting picture is as follows: The evolving population spends most of the time in a “dynamic stasis” near saddle points. During this stage, a dynamic equilibrium exists under purifying selection: Deleterious mutations constantly produce individuals with fitness lower than the current maximum, and these individuals or their progeny die out. On time scale of (kλN)⁻¹, a new beneficial mutation will occur, and then either the new type will go extinct fast (in which case, the population has to wait for another beneficial mutation) or will get fixed such that, in time lnN $/ s$ , the new type (followed by a small, dynamic cloud of suboptimal types) will constitute the bulk of the population. The transition from one most common type to the next occurs along the heteroclinic trajectory coinciding with the edge of the infinite-dimensional simplex connecting the two vertices corresponding to monotypic populations. This iterative process of fast transitions between long stasis periods spent near saddle points is typical of noisy heteroclinic networks, as demonstrated in early, semiheuristic work (34, 35, 36) and later rigorously (30, 31). However, the two types of noisy contributions, from reproduction and mutation, play distinct roles here, so although the general punctuated character of the process that we describe here is the same as in the previous studies, their results do not apply to our case straightforwardly.

Because the process is random, deviations from this general description eventually will occur. SE-unlikely, extremely rare events can be ignored. However, the right-hand side of Eq. 6, albeit small, does not decay stretch-exponentially, and so, with a nonnegligible frequency, a new beneficial mutation would appear before the current fittest type takes over the entire population. The result will be clonal interference such that the current fittest type starts being replaced with the new one before reaching fixation.

Taking the Structure of the Landscape into Account.

In general, the structure of the landscape can be complicated. The available information on the structure of complex landscapes is limited, and there are few mathematical results. Several rigorous results based on random matrix theory have been obtained for centered Gaussian fields on Euclidean spheres of growing dimension with rotationally invariant covariances of polynomial type (37, 38). For those models, the average numbers of saddles of different indices at various levels of the landscape have been shown to grow exponentially with respect to the dimension of the model, and a variational characterization of the exponential rates has been obtained. Although formally limited to concrete models, these results indicate that there are many local maxima and many more saddle points in such complex landscapes. In the context of the evolutionary process, this indicates that the evolutionary path through a sequence of temporarily most fit types is likely to end up not in a global but in a local maximum. Consider now what transpires near a local fitness peak. Suppose the current most fit genotype differs in k₀ sites from the locally optimal genotype, and sequential beneficial mutations in these sites in an arbitrary order produce a succession of increasing fitness values. Ignoring shorter times of order ln N of transitioning between saddles and only taking into account the leading contributions (that is, the sum of the waiting times for the beneficial mutations), the time it takes to reach the peak is then of the order of ${(k_{0} λ N)}^{- 1} + {((k_{0} - 1) λ N)}^{- 1} + \dots + {(2 λ N)}^{- 1} + {(λ N)}^{- 1} \approx {(λ N)}^{- 1} \ln k_{0}$ (recall that our time units are comparable with reproduction rates). Once the peak is reached, it is extremely unlikely that the population moves anywhere else on the landscape. More specifically, the waiting time for the appearance of a new most fit genotype is exponentially large in N as follows from the metastability theory at the level of large deviations estimates.

Proof of Theorem 1:

To prove Part 1, our first goal is to represent the discrepancy $D (t)$ defined in Eq. 3 in a convenient way. We can write the solution $Φ^{t} x (0)$ of ODE (1) with initial value $x (0)$ as

{(Φ^{t} x (0))}_{i} - x_{i} (0) = \int_{0}^{t} b_{i} (Φ^{s} x (0)) d s, i \in I .

[8]

It is useful to represent

x (t)

in a similar form. To that end, we recall that every Markov process solves the martingale problem associated with its own generator. Therefore, introducing the projection function

π_{i} (x) = x_{i}

, we obtain that there is a martingale

M_{i}

such that

x_{i} (t) - x_{i} (0) = π_{i} (x (t)) - π_{i} (x (0)) = \int_{0}^{t} N π_{i} (x (s)) d s + M_{i} (t), i \in I,

[9]

where the generator

N h

is defined by

N h (x) = \lim_{t ↓ 0} \frac{E [h (x (t)) | x (0) = x] - h (x)}{t} .

For our pure jump process the generator is determined by transition rates:

N h (x) = N \sum_{\begin{matrix} i, j \in T \\ i \neq j \end{matrix}} f_{i} x_{i} x_{j} (h (σ^{i j} x) - h (x)),

where

σ^{i j} x

denotes the state obtained from state

x

by adding an individual of type

i

displacing an individual of type

j

{(σ^{i j} x)}_{k} = {\begin{array}{l} x_{k}, & k \neq i, j, \\ x_{i} + \frac{1}{N}, & k = i, \\ x_{j} - \frac{1}{N}, & k = j . \end{array}

We can compute directly

\begin{array}{rc} N π_{i} (x) & = N \sum_{j : j \neq i} f_{i} x_{i} x_{j} \frac{1}{N} + N \sum_{j : j \neq i} f_{j} x_{j} x_{i} (- \frac{1}{N}) = \sum_{j} (f_{i} - f_{j}) x_{i} x_{j} = b_{i} (x) . \end{array}

Plugging this into (), we obtain

x_{i} (t) - x_{i} (0) = \int_{0}^{t} b_{i} (x (s)) d s + M_{i} (t), i \in I .

[10]

Subtracting from (10), we obtain

D_{i} (t) = x_{i} (t) - {(Φ^{t} x (t))}_{i} = \int_{0}^{t} (b_{i} (x (s)) - b_{i} (Φ^{t} x (s))) d s + M_{i} (t), i \in I .

[11]

We will view

M (t) = {(M_{i} (t))}_{i \in I}

as a vector-valued martingale. To estimate the integral term, we recall the definition () and prove the following statement:

Lemma 1.

Let $F = \max_{i \in T} f_{i}$ . Then, for all $I \subset T$ , $∥ b (x) - b (y) ∥ \leq 3 F ∥ x - y ∥, x, y \in Δ_{I} .$

Proof:

We have

∥ b (x) - b (y) ∥ = \sum_{i} | b_{i} (x) - b_{i} (y) | = \sum_{i} | (f_{i} x_{i} - x_{i} \sum_{j} x_{j} f_{j}) - (f_{i} y_{i} + y_{i} \sum_{j} y_{j} f_{j}) |

\leq J_{1} (x, y) + J_{2} (x, y),

where

J_{1} (x, y) = | \sum_{i} f_{i} (x_{i} - y_{i}) | \leq F ∥ x - y ∥

and

\begin{array}{l} J_{2} (x, y) & \leq \sum_{i} | x_{i} \sum_{j} x_{j} f_{j} - y_{i} \sum_{j} y_{j} f_{j} | \\ \leq \sum_{i} | x_{i} (\sum_{j} x_{j} f_{j} - \sum_{j} y_{j} f_{j}) + (x_{i} - y_{i}) \sum_{j} y_{j} f_{j} | \\ \leq \sum_{i} x_{i} J_{1} (x, y) + \sum_{i} | x_{i} - y_{i} | F \leq J_{1} (x, y) + F ∥ x - y ∥_{1} \leq 2 F ∥ x - y ∥ . \end{array}

Combining three displays above, we complete the proof.

Taking the absolute value in Eq. 11, then taking the sum over $i \in I$ and applying Lemma 1, we obtain

∥ D (t) ∥ \leq 3 F \int_{0}^{t} ∥ D (s) ∥ d s + M^{*} (t),

where

M^{*} (t) = \sup_{s \in [0, t]} ∥ M (s) ∥

Using the Gronwall inequality, we obtain

∥ D (t) ∥ \leq M^{*} (t) e^{3 F t} .

[12]

To estimate

M^{*} (t)

, we first use to write for any

β > 0

P \{M^{*} (t) \geq N^{- β}\} \leq \sum_{i} P \{M_{i}^{*} (t) \geq N^{- β - μ}\} \leq N^{μ} \max_{i \in I} P \{M_{i}^{*} (t) \geq N^{- β - μ}\},

[13]

where

M_{i}^{*} (t) = \sup_{s \in [0, t]} | M_{i} (s) |

. Next, we will apply an exponential martingale inequality from Appendix B6 of ref. in the form given by van de Geer (lemma 2.1 of ref. ):

Lemma 2.

If jumps of a locally square integrable cadlag martingale ${(M (t))}_{t \geq 0}$ are uniformly bounded by a constant $K > 0$ , then

P \{\exists t : | M (t) | \geq A, {⟨ M ⟩}_{t} \leq B^{2}\} \leq 2 \exp [- \frac{A^{2}}{2 (A K + B^{2})}] .

Each

M_{i}

is a piecewise linear martingale with jumps of size

1 / N

[its jumps coincide with those of

x_{i} (t)

]. Since, in addition, the total jump rate is bounded by

N F

, we obtain that the predictable quadratic variation of

M_{i}

satisfies

{〈 M_{i} 〉}_{t} \leq t N F / N^{2} = t F / N

. Thus, we can apply Lemma 2 with

B^{2} = t F / N

K = 1 / N

, and

A = N^{- β - μ}

P \{M_{i}^{*} (t) \geq N^{- β - μ}\} \leq 2 \exp [- \frac{N^{- 2 (β + μ)}}{2 (N^{- (β + μ) - 1} + t F N^{- 1})}], \cdot i \in I .

Combining this with (13), choosing

β

so that

β + μ < 1 / 2

and using

t = c \ln N

, we can find constants

C, γ > 0

such that

P \{M^{*} (t) \geq N^{- β}\} \leq 2 N^{μ} \exp [- \frac{N^{- 2 (β + μ)}}{2 (N^{- (β + μ) - 1} + t F N^{- 1})}] \leq C e^{- N^{γ}} .

Using this in (12), we complete the proof of Part 1 of the theorem. To prove Part 2, we notice that according to Part 1, up to an SE-unlikely event, the stochastic process follows the deterministic trajectory

N^{- β}

closely up to time

τ_{e} \land c l n N

, where

τ_{e}

is the first moment when one of the types goes extinct. We can restart the process at

τ_{e} \land c \ln N

treating

x (τ_{e} \land c \ln N)

as a new starting point and apply the same estimate to the restarted process (in case

τ_{e} < c \ln N

, with fewer nonzero coordinates involved). Patching several ODE trajectories together in this way and noting that, conditioned on nonextinction of type

i^{*}

, the total time it takes to travel from any point

x \in Δ_{I}

with

x_{i^{*}} \geq N^{- 1}

to the neighborhood of

e^{(i^{*})}

of size

N^{- δ}

is bounded by

C \ln N

for some

C

, we obtain Part 2.

The remaining parts follow from an auxiliary statement. To state it, we define a jump Markov process $y (t)$ with values in ${0, N^{- 1}, 2 N^{- 1} \dots, 1}$ such that $y (0) = x (0)$ and $y (t)$ makes a jump from $x$ to $x + N^{- 1}$ with rate $N f^{*} x (1 - x)$ and to $x - N^{- 1}$ with rate $N f^{* *} x (1 - x)$ , where $f^{* *} < f^{*}$ is the second largest value among $f_{i,} i \in I$ .

Lemma 3.

1. The process $y (t)$ is stochastically dominated by $x_{i^{*}} (t)$ . 2. The process $y (t)$ considered only at times of jumps is an asymmetric random walk on ${0, N^{- 1}, 2 N^{- 1} \dots, 1}$ with absorption at $0$ and $N$ and probabilities of a step to the right and left being $p$ and $1 - p$ where $p \in (1 / 2,1)$ solves

\frac{p}{1 - p} = \frac{f^{*}}{f^{* *}} .

Proof:

The coordinate $x_{i^{*}}$ jumps to the right with rate $N f_{i} x_{i^{*}} (1 - x_{i^{*}})$ and to the left with rate

N x_{i^{*}} \sum_{j \neq i^{*}} f_{j} x_{j} \leq N x_{i^{*}} f^{* *} \sum_{j \neq i^{*}} x_{j} = N f^{* *} x_{i^{*}} (1 - x_{i^{*}}) .

So, the jump rates to the left for both processes coincide and the jump rates to the right for process

y (t)

do not exceed those for process

x_{i^{*}} (t)

, and Part 1 of the lemma follows. To prove Part 2, it suffices to note that the ratio of the jump right rate to the jump left rate for process

y (t)

is equal to

f^{*} / f^{* *}

everywhere (except the absorbing points

0

and

1

To prove Part 3, we can use this lemma and the fact that if $m \geq N / 2$ , then

N \frac{m}{N} \frac{m - N}{N} \geq \frac{1}{2} (m - N),

which implies that (except for an exponentially improbable event that

x_{i^{*}}

hits level

N / 2

before

1

), the time it takes for all non-

i^{*}

types to die out is stochastically dominated by the extinction time for the linear birth-and-death process with birth rate

λ_{k} = A k

and death rate

μ_{k} = B k

where

A = f^{* *} / 2 < B = f^{*} / 2

. The probabilty

p_{k} (t)

of extinction by time

t

starting with

k

individuals was probably first computed in ref. . There is a misprint in equation 78 in ref. but one can use equation 68 of that paper (for generating functions) to obtain

p_{k} (t) = {(\frac{B e^{(B - A) t} - B}{B e^{(B - A) t} - A})}^{k} = {(1 - \frac{B - A}{B e^{(B - A) t} - A})}^{k} .

Plugging

t = C' \ln N

and

k = N^{1 - δ}

into this formula we obtain

\begin{array}{rc} 1 - p_{N^{1 - δ}} (C' \ln N) & = 1 - {(1 - \frac{B - A}{B N^{C' (B - A)} - A})}^{N^{1 - δ}} \\ \sim \frac{(B - A) N^{1 - δ}}{B N^{C' (B - A)} - A} \sim \frac{B - A}{B} N^{1 - δ - C' (B - A)}, \end{array}

and since

α = C' (B - A) - 1 + δ > 0

if we choose

C'

large enough, the desired result follows.

The last two parts of Theorem 1 follow from Lemma 3, and similar well-known statements for asymmetric random walks.

Discussion

Despite some disagreements regarding its extent, fossil record analysis suggests that PE is important in organismal evolution (7, 8, 10), which is, therefore, in general, not gradualist. Here we examine mathematically a simple population-genetic model and show that the default regime of population evolution under basic, realistic assumptions, namely, large effective population size, low mutation rate, and rarity of beneficial mutations phenomenologically resembles PE. It has to be stressed that this model is entirely within the classical framework of population genetics which also includes estimates of mutation fixation times and the waiting times between fixation events (42, 43). We reformulate it here, in order to take advantage of the mathematical toolkit of heteroclinic network analysis that provides for a rigorous treatment.

We show that, in the weak-mutation limit, large populations spend most of their time in “dynamic stasis,” that is, exercising short-range random walks within their local neutral networks in the vicinity of saddle points on the fitness landscape, without shifting to a new distinct state. The stasis periods are punctuated by rapid transitions between saddle points upon emergence of new beneficial mutations; these transitions appear effectively instantaneous compared to the duration of stasis, even when they evolve through more than one elementary mutation event (Fig. 2). Eventually, the population might reach a local fitness peak where no beneficial mutations are available. This would lead to indefinite stasis as long as the fitness landscape does not change and the population size stays large (drift to a different peak is exponentially rare in N, that is, impractical for large N).

Fig. 2.

Evolution under PE on a fitness landscape dominated by saddles: stasis around saddle points punctuated by fast adaptive transitions. Planar shapes depict distinct classes of genotypes. The color scale shows a range of fitness values. Gray “ramp” strips show available transitions between the genotype classes (k transitions leading to classes with higher fitness and L transitions leading to classes with lower fitness, $k ≪ L$ ). The two blue circles indicate the original and the current states of the population; blue arrows show succession of genotypes within the same class, occurring within the effectively neutral network during the “dynamic stasis” phase; red arrows indicate fast adaptive transitions from a lower-fitness genotype to one with a higher fitness.

Two conditions determine the behavior described by this model: 1) low overall mutation rate (dominated by deleterious mutations), [Eq. (7)], $λ L ≪ s / \ln N$ and 2) an even lower rate of beneficial mutations, which results in the difference in scale between the waiting time ${(λ k N)}^{- 1}$ and the saddle-to-saddle transition time $\ln N / s$ , that is, $λ k N ≪ s / \ln N$ . Comparison of the expressions for these conditions suggests that, for the PE-like regime to be pronounced, in the vicinity of all saddle points where the population is at equilibrium, deleterious mutations should outnumber the beneficial mutations by at least a factor of $N$ . This is a large but not unrealistic difference in the case of “highly adapted” organisms, that is, in situations that are most common in the extant biosphere, where the pool of straightforward optimizations that presumably were available at the earliest stages of the evolution of life had been exhausted long ago. For example, with population and genomic parameters characteristic of animals, namely, N ∼10⁵ and ∼10⁷ amino acid-encoding sites in the genome, the local mutational neighborhood in the sequence space consists of 19 × 10⁷ mutations. Assuming that about half of these mutations are deleterious and noting that the number of beneficial mutations should be less by a factor of 10⁵, there must be 1 < k < 1,000 beneficial mutations available, apparently, a realistic value. It should be emphasized that, once a beneficial mutation occurs it can open the door for additional ones due to pervasive epistasis, so that a transition might involve multiple adaptive mutational events.

The condition on the overall mutation rate ( $λ L ≪ s / \ln N$ ) is more difficult to assert because both $λ$ and $L$ depend on the clustering of the whole sequence space into a coarse-grained network of distinct types. Note, however, that, as the first approximation, $λ$ is bounded by the sequence-level mutation rate $μ$ (only some of the sequence-level mutations lead to transitions between distinct types) and $L$ is bounded by the genome size $G$ (the number of available sequence-level single-position mutations is on the order of the genome size, although only some of these mutations have detectable deleterious effect). Thus, $λ L < μ G$ , where $μ G$ is the expected number of sequence-level mutations per genome per generation. It has been shown that the values of $μ G$ tend to stay of the order of $1 / N$ under “normal” conditions (44, 45); therefore,

λ L < μ G \sim 1 / N ≪ s / \ln N .

The condition for a mutation to be subject to selection in an asexual population is Ns >>1 (). According to the above inequality, the weak-mutation regime applies when

N s ≫ \ln N

. For realistic population sizes, N < 10⁹, lnN < 20, so this regime holds for a broad range of conditions.

Thus, our model suggests that the PE-like regime is common and is likely to be the default in the evolution of natural populations. The probable exceptions include stress-induced mutagenesis (46), whereby the mutation rate can rise by orders of magnitude, locally blooming microbial populations that might violate the $k N ≪ L$ condition, and abrupt changes in the fitness landscape that might temporarily increase the number of immediately beneficial mutations $k$ . All of these situations, however, are likely to be transient.

Theoretically, PE has been linked to SOC as the underlying mechanism (12, 15). However, we show here that a PE-like regime is readily observed in extremely simple models of population evolution that do not involve any criticality. The major conclusion from this analysis is that PE-like evolution rather than gradualism is the fundamental character of sufficiently large populations in the weak-mutation limit which is, arguably, the most common evolutionary regime across the entire diversity of life. The parameter values that lead to this regime appear to hold for evolving populations of all organisms, including viruses, under “normal” conditions. Situations can emerge in the course of evolution when the PE regime breaks through disruption of the stasis phase. This could be the case in very small populations that rapidly evolve via drift or in cases of a dramatically increased mutation rate, such as stress-induced mutagenesis, and especially when these two conditions combine (46 47–48). In many cases, disruption of stasis will lead to extinction but, on occasion, a population could move to a different part of the landscape, potentially, the basin of attraction of a higher peak. The evolution of cancers, at least at advanced stages, does not appear to include stasis either, due to the high rate of nearly neutral and deleterious mutations and low effective population size (46). Furthermore, the PE-like regime is characteristic of “normal” evolution of well-adapted populations in which the fraction of beneficial mutations is small. If many, perhaps the majority, of the mutations are beneficial, there will be no stasis but rather a succession of rapid transitions in a fast adaptive evolution regime. Conceivably, this was the mode of evolution of primordial replicators at precellular stages of evolution.

One of the most fundamental—and most difficult—problems in biology is the origin of major biological innovations (more or less synonymous to macroevolution). In modern evolutionary biology, Darwin’s central idea of survival of the fittest transformed into the concept of fitness landscape with numerous peaks, where each stable form occupies one of the peaks (24, 49). Then, the fundamental problem arises: If a population has reached a local peak further adaptive evolution is possible only via a stage of temporary decrease of fitness. How can this happen? A common answer is based on Wright’s concept of random genetic drift: The smaller the effective population size N_e (or simply N, for a well-mixed population) the greater the probability of random drift through (not excessively deep) valleys in the fitness landscape (49 50–51). This notion implies that evolutionary transitions occur through narrow population bottlenecks. As formalized in our previous work, the evolutionary “innovation potential” is inversely proportional to N_e (22). There are, however, multiple indications that drift is unlikely to be the only mode of evolutionary innovation and that novelty often arises in large populations thanks to their high mutational diversity (52 53 54–55). Nevertheless, it remains unclear, within the tenets of classical population genetics, how a large population can cross a valley on the landscape. One obvious way to overcome this conundrum is to assume that the landscape changes in time due to environmental changes, so that peaks could become saddle points, and vice versa, and a population might find itself in the basin of attraction of a new fitness peak (56, 57). The analysis presented here suggests a greater innovation potential of large populations than usually assumed, stemming from the fact that a typical landscape in a multidimensional space contains many more saddle points than peaks. On the one hand, this intuitively obvious claim follows from the observation that, for any two peaks, the path connecting the peaks and maximizing the minimum height must pass through a saddle point. On the other hand, it is justified by precise computations of exponential (with respect to the model dimension) growth rates of the expected numbers of saddle points of various indices (including peaks) for random Gaussian landscapes under certain restrictions on covariance (37, 38). Thus, typical fitness landscapes are likely to allow numerous transitions and extensive, innovative evolution without the need for valley crossing, as also argued previously from the analysis of “holey” fitness landscapes (24). In biological terms, it seems to be impossible to maximize fitness in all numerous directions (the number of these being at least on the order of the genome size), and therefore the probability of beneficial mutations is (almost) never zero, however small it might be (in general, this pertains not only to single point mutations but also to beneficial epistatic combinations of mutations as well as large-scale genomic changes, such as gene gain, loss, and duplication). In other words, the landscape is dominated by saddle points that are far more common than peaks, so that there is almost always an upward path which an evolving population will follow provided it is large enough to afford a long wait in saddles without risking extinction due to fluctuations.

Results similar to ours have been reported in the mathematical biology literature (27 28–29). Specifically, it has been proven that a trait substitution sequence process (sequential transition from one dominant trait to another) occurs in the limit of large population size and small beneficial mutation rate. Here we employ a very simple model to demonstrate the fundamental character of the concept of PE, to tie it to the noisy dynamics near heteroclinic networks (30, 31) and to stress the key role of saddle points, in contrast to the widespread perception of peaks as the central structural elements of fitness landscapes.

To conclude, the results presented here show that PE-like evolution is not only characteristic of speciation or evolutionary transitions but rather is the default mode of evolution under weak-mutation limit which is the most common evolutionary regime (25). In our previous work, we have identified conditions under which saltational evolution becomes feasible, under the strong-mutation limit (48). Here we show that, even for evolution in the weak-mutation limit that is generally perceived as gradual (25), PE is the default regime. Even during periods of stasis in phenotypic evolution, the underlying microevolutionary process appears to be punctuated.

Acknowledgements

Y.I.W. and E.V.K. are supported by the Intramural Research Program of the NIH. Y.B. is partially supported by NSF Grant DMS-1811444. M.I.K. was supported by Spinoza Prize funds.

The authors declare no competing interest.

Data Availability.

There are no data underlying this work.

References

Darwin

, On the Origin of Species (A.F. Murray, London, 1859).

Huxley

T. H.

, Darwin on the origin of species. Westminster Rev., 541–570 (1860).

Goldschmidt

R. B.

, The Material Basis of Evolution (Yale University Press, New Haven, CT, 1940).

Simpson

G. G.

, Tempo and Mode in Evolution (Columbia University Press, New York, 1983).

Eldredge

, Gould

S. J.

, “Punctuated equilibria: An alternative to phyletic gradualism” in Models in Paleobiology, Schopf

T. J. M.

, Ed. (Freeman Cooper, San Francisco, 1972), pp. 193–223.

Gould

S. J.

, Eldredge

, Punctuated equilibrium: The tempo and mode of evolution reconsidered. Paleobiology 3, 115–151 (1977).

Gould

S. J.

, Eldredge

, Punctuated equilibrium comes of age. Nature 366, 223–227 (1993).

Eldredge

, Gould

S. J.

, On punctuated equilibria. Science 276, 338–341 (1997).

Gould

S. J.

, Tempo and mode in the macroevolutionary reconstruction of Darwinism. Proc. Natl. Acad. Sci. U.S.A. 91, 6764–6771 (1994).

Gould

S. J.

, The Structure of Evolutionary Theory (Harvard University Press, Cambridge, MA, 2002).

Eldredge

, The dynamics of evolutionary stasis. Paleobiology 31, 133–145 (2005).

Bak

, How Nature Works: The Science of Self-Organized Criticality (Springer, New York, 1996).

Bak

, Tang

, Wiesenfeld

, Self-organized criticality: An explanation of the 1/f noise. Phys. Rev. Lett. 59, 381–384 (1987).

Bak

, Tang

, Wiesenfeld

, Self-organized criticality. Phys. Rev. A Gen. Phys. 38, 364–374 (1988).

Bak

, Sneppen

, Punctuated equilibrium and criticality in a simple model of evolution. Phys. Rev. Lett. 71, 4083–4086 (1993).

Maslov

, Paczuski

, Bak

, Avalanches and 1/f noise in evolution and growth models. Phys. Rev. Lett. 73, 2162–2165 (1994).

Maslov

, Zhang

Y. C.

, Exactly solved model of self-organized criticality. Phys. Rev. Lett. 75, 1550–1553 (1995).

Bak

, Paczuski

, Complexity, contingency, and criticality. Proc. Natl. Acad. Sci. U.S.A. 92, 6689–6696 (1995).

Szathmáry

, Smith

J. M.

, The major evolutionary transitions. Nature 374, 227–232 (1995).

Maynard Smith

, Szathmary

, The Major Transitions in Evolution (Oxford University Press, Oxford, 1997).

Szathmáry

, Toward major evolutionary transitions theory 2.0. Proc. Natl. Acad. Sci. U.S.A. 112, 10104–10111 (2015).

Katsnelson

M. I.

, Wolf

Y. I.

, Koonin

E. V.

, Towards physical principles of biological evolution. Phys. Scr., 93043001 (2018).

Wolf

Y. I.

, Katsnelson

M. I.

, Koonin

E. V.

, Physical foundations of biological complexity. Proc. Natl. Acad. Sci. U.S.A. 115, E8678–E8687 (2018).

Gavrilets

, Fitness Landscapes and the Origin of Species (Princeton University Press, Princeton, NJ, 2004).

Gillespie

J. H.

, The Causes of Molecular Evolution (Oxford University Press, Oxford, 1994).

Moran

P. A.

, Random processes in genetics. Proc. Philos. Soc. Math. and Phys. Sci. 54, 60–71 (1958).

Champagnat

, A microscopic interpretation for adaptive dynamics trait substitution sequence models. Stochastic Process. Appl. 116, 1127–1160 (2006).

Champagnat

, Meleard

, Polymorphic evolution sequence and evolutionary branching. Probab. Theory Relat. Fields 151, 45–94 (2011).

Kraut

, Bovier

, From adaptive dynamics to adaptive walks. J. Math. Biol. 79, 1699–1747 (2019).

Bakhtin

, Small noise limit for diffusions near heteroclinic networks. Dyn. Syst. 25, 413–431 (2010).

Bakhtin

, Noisy heteroclinic networks. Probab. Theory Relat. Fields 150, 1–42 (2011).

Nowak

M. A.

, Evolutionary Dynamics: Exploring the Equations of Life (Belknap Press, Cambridge, MA, 2006).

Karev

G. P.

, How to explore replicator equations? arXiv:0812.4295v1 (22 December 2008).

Stone

, Holmes

, Random perturbation of heteroclinic attractors. SIAM J. Appl. Math. 50, 726–743 (1990).

Stone

, Armbruster

, Noise and O(1) amplitude effects on heteroclinic cycles. Chaos 9, 499–506 (1999).

Armbruster

, Stone

, Kirk

, Noisy heteroclinic networks. Chaos 13, 71–86 (2003).

Auffinger

, Ben Arous

, Complexity of random smooth functions on the high-dimensional sphere. Ann. Probab. 41, 4214–4247 (2013).

Ben Arous

, Mei

, Montanari

, Nica

, The landscape of the spiked tensor model. Commun. Pure Appl. Math. 72, 2282–2330 (2019).

Shorack

G. R.

, Wellner

J. A.

, Empirical Processes with Applications to Statistics (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2009).

Van de Geer

, Exponential inequalities for martingales, with application to maximum likelihood extimation for counting processes. Ann. Stat. 23, 1779–1801 (1995).

Bartholomay

A. F.

, On the linear birth and death processes of biology as Markoff chains. Bull. Math. Biophys. 20, 97–118 (1958).

Ewens

W. J.

, Population Genetics (Methuen, New York, 1969).

Crow

J. F.

, Kimura

, Introduction to Population Genetics Theory (Harper and Row, New York, 1970).

Lynch

, Evolution of the mutation rate. Trends Genet. 26, 345–352 (2010).

Lynch

, Genetic drift, selection and the evolution of the mutation rate. Nat. Rev. Genet. 17, 704–714 (2016).

Fitzgerald

D. M.

, Hastings

P. J.

, Rosenberg

S. M.

, Stress-induced mutagenesis: Implications in cancer and drug resistance. Annu. Rev. Cancer Biol. 1, 119–140 (2017).

Ram

, Hadany

, Evolution of stress-induced mutagenesis in the presence of horizontal gene transfer. Am. Nat. 194, 73–89 (2019).

Katsnelson

M. I.

, Wolf

Y. I.

, Koonin

E. V.

, On the feasibility of saltational evolution. Proc. Natl. Acad. Sci. U.S.A. 116, 21068–21075 (2019).

Wright

, Adaptation and Selection: Genetics, Paleontology and Evolution (Princeton University Press, Princeton, NJ, 1949).

Lynch

, The Origins of Genome Archiecture (Sinauer Associates, Sunderland, MA, 2007).

Lynch

, Conery

J. S.

, The origins of genome complexity. Science 302, 1401–1404 (2003).

Masel

, Cryptic genetic variation is enriched for potential adaptations. Genetics 172, 1985–1991 (2006).

Rajon

, Masel

, Compensatory evolution and the origins of innovations. Genetics 193, 1209–1220 (2013).

Lynch

, Abegg

, The rate of establishment of complex adaptations. Mol. Biol. Evol. 27, 1404–1414 (2010).

Lynch

, Phylogenetic divergence of cell biological features. eLife 7, e34820 (2018).

Gavrilets

, Vose

, Dynamic patterns of adaptive radiation. Proc. Natl. Acad. Sci. U.S.A. 102, 18040–18045 (2005).

Mustonen

, Lässig

, From fitness landscapes to seascapes: Non-equilibrium dynamics of selection and adaptation. Trends Genet. 25, 111–119 (2009).