Survival of the flattest in the quasispecies model

Abstract

Viruses present an amazing genetic variability. An ensemble of infecting viruses, also called a viral quasispecies, is a cloud of mutants centered around a specific genotype. The simplest model of evolution, whose equilibrium state is described by the quasispecies equation, is the Moran–Kingman model. For the sharp-peak landscape, we perform several exact computations and derive several exact formulas. We also obtain an exact formula for the quasispecies distribution, involving a series and the mean fitness. A very simple formula for the mean Hamming distance is derived, which is exact and does not require a specific asymptotic expansion (such as sending the length of the macromolecules to $\infty$ or the mutation probability to 0). With the help of these formulas, we present an original proof for the well-known phenomenon of the error threshold. We recover the limiting quasispecies distribution in the long-chain regime. We try also to extend these formulas to a general fitness landscape. We obtain an equation involving the covariance of the fitness and the Hamming class number in the quasispecies distribution. Going beyond the sharp-peak landscape, we consider fitness landscapes having finitely many peaks and a plateau-type landscape. Finally, within this framework, we prove rigorously the possible occurrence of the survival of the flattest, a phenomenon which was previously discovered by Wilke et al. (Nature 412, 2001) and which has been investigated in several works (see e.g. Codoñer et al. (PLOS Pathogens 2, 2006), Franklin et al. (Artificial Life 25, 2019), Sardanyés et al. (J. Theoret. Biol. 250, 2008), and Tejero et al. (BMC Evolutionary Biol. 11, 2011)).

1. Introduction

Some viruses are very hard to fight because of their genetic diversity: they mutate constantly and this allows them to escape the attacks of the immune system of their host. At the same time, the genotypes of an infecting population of viruses are strongly correlated; they look like a cloud of mutants centered around a specific genotype. This population structure, called a quasispecies, can be observed experimentally for viruses (in vivo studies have been conducted for HIV [17] and the hepatitis C virus [12]). Nowadays, deep sequencing techniques make it possible to collect data on the structure of viral quasispecies. Yet the various biological mechanisms involved in the creation and evolution of quasispecies are far from being understood (see [8] for a recent review or the book [9] for a more comprehensive account). The design of simple mathematical models of quasispecies may be helpful in developing efficient antiviral therapies. In principle, a quasispecies can occur in any biological population whose evolution is driven by mutation and selection. In the next subsection, we present a very simple model for the evolution of a population under mutation and selection. We then write the equations describing the equilibrium of this model, thereby recovering directly the quasispecies equations associated with the equilibrium of Eigen’s model. All the results presented afterwards concern the solutions of these equations.

1.1. The quasispecies equations

We consider a population of haploid individuals evolving under the conjugate effects of mutation and selection. We suppose that the population has reached an equilibrium around a specific well-adapted individual, called the wild type, and denoted by $w^*$ . Individuals reproduce, yet the reproduction mechanism is error-prone, and mutations occur constantly. These mutations drive the genotypes away from $w^*$ . Yet $w^*$ has a selective advantage because it reproduces faster.

We would like to characterize this kind of equilibrium mathematically. We denote by E the set of possible genotypes, which we assume to be finite. Generic elements of E are denoted by the letters u, v. The Darwinian fitness of an individual having genotype u is denoted by f(u); it can be thought of as its mean number of offspring. In addition, mutations occur in each reproduction cycle: an individual of type u might appear as the result of mutations from offspring of other types. Let us denote by M(v, u) the probability that the offspring of an individual of type v is of type u. For simplicity, we will assume that all the coefficients of the mutation matrix M are positive. Of course, we have

\begin{align*}\forall\, v\in E{\qquad}\sum_{u\in E}M(v,u)=1.\end{align*}

We introduce next the linear model, one of the simplest models for the evolution of a population with selection and mutation. We suppose that the successive generations do not overlap, and we denote by $N_n(u)$ the number of individuals of type u in generation n. The linear model assumes that an individual of type v produces a number of offspring proportional to its fitness f(v), and that a proportion M(v, u) of the offspring mutates and becomes of type u; thus $N_{n+1}(u)$ is given by the formula

\begin{align*}\forall\,u\in E{\qquad} N_{n+1}(u)\,=\,\sum_{v\in E}N_n(v)f(v)M(v,u).\end{align*}

The trouble with this formula is that the sum is not necessarily an integer. A natural way to circumvent this problem is to develop stochastic population models, in such a way that the above formula describes the evolution of the mean number of individuals. The archetype of this kind of model is the Galton–Watson branching process. If we introduce in addition a constraint on the total size of the population, then we arrive at the classical Wright–Fisher model. Yet the randomness adds an extra layer of complexity; stochastic models are considerably harder to study. Another, simpler, possibility is to consider the proportion of each type of individual in the population, instead of the number, as Moran and Kingman did in the late seventies [16, 18]. Let us denote by $x_n(u)$ the proportion of individuals of type u in generation n. The model proposed by Moran is given by

\begin{equation*}\forall\, u\in E{\qquad} x_{n+1}(u)\,=\,\frac{\sum_{v\in E}x_n(v)f(v)M(v,u)}{\sum_{v\in E}x_n(v)f(v)}.\end{equation*}

The denominator has been chosen to ensure that $\sum_{u\in E}x_{n+1}(u)=1$ . In fact, this model consists in iterating a deterministic map on the function encoding the proportion of each type present in the population. The possible equilibria of the model correspond to the fixed points of this map, that is, the solutions of the following equations:

(1.1) \begin{equation}\forall\, u\in E{\qquad} x(u)\,\sum_{v\in E}x(v)f(v)\,=\,\,\sum_{v\in E}x(v)f(v)M(v,u)\,\end{equation}

subject to the constraint

(1.2) \begin{equation}\forall\, u\in E{\qquad} x(u)\geq 0\,,{\qquad}\sum_{u\in E}x(u)\,=\,1.\end{equation}

We call these equations the quasispecies equations. They also characterize the equilibrium in the model originally developed by Eigen [10], which led to quasispecies theory. Some of our results are derived for the general quasispecies equations (1.1). To go further, we will focus on a particular choice of the set of genotypes E and of the mutation matrix M. For both practical and historical reasons, we make the same choice as Eigen and Schuster [11], which leads to the sharp-peak landscape.

Genotypes. We consider the different genotypes to be sequences of length $\ell\geq1$ over an alphabet $\mathcal{A}$ of cardinality $\kappa\geq 1$ . Standard examples are $\mathcal{A}=\{\,0,1\,\}$ , $\kappa=2$ , $E=\lbrace\,0,1\,\rbrace^\ell$ , for binary sequences, and $\mathcal{A}=\lbrace\,A, C, G, T\,\rbrace$ , $\kappa=4$ , $E=\lbrace\,A,C,G,T\,\rbrace^\ell$ , to model DNA molecules. The ith component of a genotype u is denoted by u(i). The Hamming distance d counts the number of differences between two sequences:

\begin{align*}\forall\,u,v\in E{\qquad}d(u,v)\,=\,\text{card}\,\big\lbrace\,1\leq i\leq \ell\;:\;u(i)\neq v(i)\,\big\rbrace.\end{align*}

When we wish to work with the simplest model, we use binary sequences. When we seek to state a more general result, we use sequences over a finite alphabet.

Mutations. We suppose that mutations happen independently along the sequence with probability $q\in\,]0,1[\,$ , and that, upon mutation, a new letter of the alphabet is chosen with equal probability. For $u,v\in E$ , the mutation probability M(u, v) is given by

(1.3) \begin{equation}M(u,v)\,=\,\Big(\frac{q}{\kappa-1}\Big)^{d(u,v)}(1-q)^{\ell-d(u,v)}.\end{equation}

The final ingredient needed to completely specify the model is the fitness function. We will consider different types of fitness functions in the following subsections. For the initial models, for which the fitness function is quite simple, we work with the genotype space $E=\mathcal{A}^\ell$ where $\mathcal{A}$ is a finite alphabet. For the more complicated fitness landscapes considered at the end, which contain a kind of plateau, we work with the genotype space of binary sequences $E=\lbrace\,0,1\,\rbrace^\ell$ .

1.2. The sharp-peak landscape

Ideally, we would like to have explicit formulas for x in terms of f and M. Unfortunately, there is little hope of obtaining such explicit formulas in the general case. So we consider the simplest non-neutral fitness function which comes to mind, the sharp-peak landscape: there is a privileged genotype, $w^*$ , referred to as the wild type, which has a higher fitness than the rest. We work with the genotype space $E=\mathcal{A}^\ell$ where $\mathcal{A}$ is a finite alphabet of cardinality $\kappa$ . Let $\sigma>1$ , and let the fitness function $f_{\text{SP}}$ be given by

(1.4) \begin{equation}\forall\,u\in E{\qquad}f_{\text{SP}}(u)\,=\,\begin{cases}\quad\sigma&\quad\text{if}\quad u=w^*\,,\\[5pt] \quad1&\quad\text{if}\quad u\neq w^*\,.\end{cases}\end{equation}

This is the fitness function that Eigen studied in detail in his article [10]. Despite its apparent simplicity, this model leads to interesting mathematical questions, and it provides new insight into the genetic structure of a population. Eigen and Schuster [11] discussed the sharp-peak landscape (which is fully presented in Section 2.2) with the help of approximation techniques and in a specific asymptotic regime. We present here a new approach, which is more elementary. The sharp-peak landscape is studied in detail in Section 2. We perform several exact computations and derive several exact formulas. In particular, we obtain a remarkable equation for the mean fitness $\lambda$ at equilibrium, which was previously discovered by Bratus et al.; see the formulas $5.5$ and $5.6$ in [3]. We state next a probabilistic version of this formula. Let $S_\ell$ be a random variable with distribution the binomial law with parameters $\ell$ and $(\kappa-1)/\kappa$ . Let $\mathbb{E}$ denote the expectation; then the mean fitness $\lambda$ at equilibrium satisfies

(1.5) \begin{equation}\frac{1}{\sigma-1}\,=\, \mathbb{E}\Bigg(\frac{1}{\frac{\lambda}{\rho^{S_\ell}}-1}\Bigg),\end{equation}

where we have set

(1.6) \begin{equation}\rho\,=\,1-\frac{\kappa q}{\kappa-1}.\end{equation}

Bratus et al. exploited this equation to derive rigorous bounds on the mean fitness, but this formula has wider applications. Furthermore, in Section 2.4 we obtain an exact formula for the quasispecies distribution, involving a series and the mean fitness. A very simple formula for the mean Hamming distance is derived, which is exact and which does not require a specific asymptotic expansion, such as sending the length of the macromolecules to $\infty$ or the mutation probability to 0 (recall that the Hamming distance between two sequences is the number of sites where they differ). The mean Hamming distance $\overline{\mathcal{Q}}$ in the quasispecies distribution can be interpreted as the mean number of ones in a genotype at equilibrium, and the formula is

(1.7) \begin{equation}\overline{\mathcal{Q}} \,=\,\frac{\ell q\, \lambda} {\lambda-\rho},\end{equation}

where $\lambda$ is the mean fitness of the population (here we assume that the wild type is the genotype $(0,\dots,0)$ and the Hamming distances are computed with respect to the wild type).

We stress that all the formulas obtained in Section 2 are exact and hold for a fixed value of the length $\ell$ . In the next subsection and in Section 3, we perform asymptotic expansions of these formulas in the long-chain regime.

1.3. The error threshold

We shall demonstrate that the quasispecies equations for the sharp-peak landscape undergo a phenomenon similar to a phase transition. This is a classical and well-known result. However, we adopt an original approach for the proof; indeed, we will show that it follows easily from Equation (1.5). Furthermore, with the help of the exact formulas derived in Section 2.4, we recover the limiting quasispecies distribution in the long-chain regime in Section 3.3.

We consider the asymptotic regime, called the long-chain regime, where

(1.8) \begin{equation}\ell\to+\infty\,,{\qquad}q\to 0\,,{\qquad}\ell q\to a\in\,[0,+\infty].\end{equation}

The parameter a is in fact the asymptotic mean number of observed mutations per individual in each reproduction cycle. We denote by $\lambda$ the mean fitness at equilibrium and by $x(w^*)$ the proportion of the wild type $w^*$ in the population at equilibrium.

Theorem 1. (Error threshold.) For the sharp-peak landscape, we have the following dichotomy:

• • If $a>\ln \sigma$ , then $\lambda\rightarrow 1$ , whence $x(w^*)\rightarrow 0$ .

• • If $a<\ln \sigma$ , then $\lambda\rightarrow \sigma e^{-a}>1$ , and

  \begin{equation*}x(w^*)\quad\to\quad\frac{\sigma e^{-a}-1}{\sigma-1}.\end{equation*}

Thus the limit of the proportion $x(w^*)$ of the wild type present in the population is positive only if $\sigma e^{-a}>1$ , or equivalently, if asymptotically $q<\ell^{-1}\ln\sigma$ . The value

\begin{align*}q^*\,=\,\frac{\ln\sigma}{\ell}\end{align*}

is the error threshold, described by Eigen as the critical mutation rate above which the wild type $w^*$ disappears from the population. Whenever the occurrence of mutations is negligible compared to the sequence length, that is, for $a=0$ , we have $x(w^*)=1$ and the only genotype present in the population is the wild type $w^*$ . In the presence of a small number of mutations, that is, for $a>0$ , genetic diversity is constantly reintroduced in the population and there is a positive proportion of genotypes that differ from $w^*$ . Yet most of the genotypes are very close to the wild type $w^*$ . In fact, the population looks like a cloud of mutants centered around the wild type. As the mutation rate q is raised, the proportion of the wild type $w^*$ present in the population decreases. When q reaches the error threshold, the wild type disappears from the population. More precisely, the proportion of the wild type becomes comparable to the proportions of the other genotypes. This is the error catastrophe: the selective advantage of the wild type $w^*$ is annihilated by the mutations.

This kind of equilibrium was discovered within the framework of Eigen’s model and was called a quasispecies [10, 11], as opposed to a species, which refers to a homogeneous solution in chemistry. In fact, Eigen’s model is a model for the evolution of a population of macromolecules governed by a system of chemical reactions, and the laws of kinetics yield a differential system of equations whose equilibrium is described by the quasispecies equation. This system was historically analyzed with approximation and expansion techniques, which in the asymptotic regime (1.8) led to the discovery of the error threshold.

Although the original goal of Eigen was to understand the first stages of life on Earth, the notion of quasispecies and the error threshold had a profound impact on the understanding of molecular evolution [8]. Indeed, many living organisms seem to satisfy approximately the scaling relation

\begin{align*}\text{mutation rate}\quad\sim\quad\frac{1}{\text{genome length}}.\end{align*}

Unfortunately, for complex creatures, it is very difficult to estimate the mutation rate, which is usually extremely small. Viruses, however, have a rather high mutation rate, and the orders of magnitude of their genome length and their mutation rate are compatible with this scaling law. Moreover, some RNA viruses, such as HIV, evolve with a high mutation rate which seems to be close to an error threshold [8]. Why is that so?

In order to survive, a virus must achieve two goals. First, its genetic information must be preserved from one generation to another; hence its mutation rate has to be below the error threshold. Second, it has to escape the attacks of the immune system of its host, and to do so, it should mutate as fast as possible. The most efficient strategy is therefore to adjust the mutation rate to the left of the error threshold: this will achieve simultaneously a huge genetic variability, and the preservation of the pertinent genetic information across generations. Some promising antiviral strategies, called lethal mutagenesis, consist in using mutagenic drugs which slightly increase the mutation rate of the virus in order to induce an error catastrophe [1, 7, 8].

Most of the rigorous mathematical analysis of Eigen’s model deals with the sharp-peak landscape, or with landscapes presenting a lot of symmetries (see [2] for the analysis of a permutation-invariant fitness landscape and [4] for a review of the literature).

1.4. Extension to a general landscape

We try to extend the formulas obtained in Section 2 to a general fitness landscape. Here we work with the genotype space $E=\mathcal{A}^\ell$ where $\mathcal{A}$ is a finite alphabet of cardinality $\kappa$ . We make the following hypothesis on the fitness function f.

Hypothesis 1. The fitness function f is larger than or equal to 1, and it is not identically equal to 1.

In Theorem 2, we obtain an interesting equation involving the covariance of the fitness and the sum of an arbitrary function g over the letters of the sequences in the quasispecies distribution. This formula is a considerable generalization of the formula (1.7).

Theorem 2. Let $\mathcal{A}$ be a finite alphabet of cardinality $\kappa$ and let $E=\mathcal{A}^\ell$ . Let g be a function defined on $\mathcal{A}$ with values in $\mathbb{R}$ , and let $G\;:\;E\to\mathbb{R}$ be defined by

\begin{align*}\forall u\in E{\qquad} G(u)\,=\,\sum_{i=1}^\ell g(u(i)).\end{align*}

For any fitness function f over the genotype set E which is larger than or equal to 1, and which is not identically equal to 1, we have

(1.9) \begin{equation} \overline{f}\,\overline{G}\,=\,\rho\,\overline{fG}+\ell(1-\rho)\overline{f}\,\overline{g},\end{equation}

where $\rho$ is given in (1.6) and

\begin{gather*} \overline{fG}\,=\,\sum_{u\in E}f(u)\,G(u)\,x(u)\,,\quad \overline{g}\,=\,\frac{1}{\kappa}\sum_{a\in\mathcal{A}}g(a)\,,\cr\overline{f}\,=\,\sum_{u\in E}f(u)x(u)\,,\quad\overline{G}\,=\,\sum_{u\in E}G(u)x(u)\,.\end{gather*}

Let us look at the specific case of the sharp-peak landscape. We take $\kappa=2$ , $E=\lbrace\,0,1\,\rbrace^\ell$ , and the fitness function defined in the formula (1.4). For the function g, we take the identity function on $\lbrace\,0,1\,\rbrace$ , so that, for any $u\in E$ , the value G(u) is simply the Hamming class of u (that is, the number of ones present in u). In particular, we have $\overline{fG}\,=\,\overline{G}$ and $\overline{g}=1/2$ , so that the formula (1.9) becomes

(1.10) \begin{equation} \overline{f}\,\overline{G}\,=\,\rho\,\overline{G}+\frac{\ell}{2}(1-\rho)\overline{f},\end{equation}

which is precisely the formula (1.7) that we obtained previously in this context. We can also consider the neutral case. If the fitness function is constant, then the quasispecies distribution is in fact the uniform distribution over the genotype space, so that $\overline{G}=\ell/2$ , and the above relation (1.10) still holds. This is not too surprising, because it amounts to the case where the fitness $\sigma$ of the wild type $w^*$ is equal to 1.

The proof of Theorem 2 is given in Section 4, where we also perform several useful computations valid for a general state space E. The various formulas obtained in Section 4 shed new light on the classical phenomenon of the error threshold and the notion of quasispecies, which are discussed in Section 3.2. Finally, in Section 4, we outline a robust strategy for analyzing more general fitness landscapes. This strategy will be used to study the fitness landscapes introduced in the following subsections.

1.5. Finitely many peaks

Now that we have thoroughly understood the case of the sharp-peak landscape, we are ready to consider more complex landscapes. Here we work with the genotype space $E=\mathcal{A}^\ell$ where $\mathcal{A}$ is a finite alphabet of cardinality $\kappa$ . We make the following hypothesis on the fitness function f.

Hypothesis 2. There exist a fixed integer k and k fixed values $\sigma_1,\dots,\sigma_k$ such that the following holds: for any values of $\ell,q$ , there exist k peaks $w_1^*,\cdots,w_k^*$ in E (whose locations may vary with $\ell$ and q) such that the fitness function $f_{\textit{FP}}$ is given by

(1.11) \begin{equation}\forall\,u\in E{\qquad}f_{\textit{FP}} (u)\,=\,\begin{cases}\quad\sigma_i&\quad\textit{if}\quad u=w^*_i\,,\quad 1\leq i\leq k\,,\\[5pt] \quad1&\quad\textit{otherwise}\,.\end{cases}\end{equation}

We consider again the long-chain regime

\begin{equation*}\ell\to+\infty\,,{\qquad}q\to 0\,,{\qquad}\ell q\to a\in\,[0,+\infty].\end{equation*}

To simplify the proofs, we suppose that $a<+\infty$ . Indeed, the case $a=+\infty$ should be considered separately. We shall prove the following theorem.

Theorem 3. (Error threshold for finitely many peaks.) Let us set

\begin{align*}\sigma\,=\,\max_{1\leq i\leq k} \sigma_i.\end{align*}

Let $\lambda_{\textit{FP}}(\ell,q)$ be the mean fitness for the model associated with the fitness function $f_{\textit{FP}}$ and the parameters $\ell,q$ . We have the following dichotomy:

• • If $a>\ln\sigma$ , then $\lambda_{\textit{FP}}(\ell,q)\rightarrow 1$ .

• • If $a<\ln \sigma$ , then $\lambda_{\textit{FP}}(\ell,q)\rightarrow \sigma e^{-a}>1$ .

We see that the situation is very similar to the case of the sharp-peak landscape. The critical parameter depends now on the fitness of the highest peak. The conclusion is that, when there are finitely many peaks, the lower peaks do not really interfere with the highest peak. Two proofs of Theorem 3 are presented in Section 5.

1.6. Plateau

We have seen that a finite number of small peaks do not interact significantly with the highest peak. We consider here a fitness landscape which contains a kind of plateau, that is, a large subset of genotypes having a fitness $\sigma>1$ . More precisely, we consider the genotype space of binary sequences $E=\lbrace\,0,1\,\rbrace^\ell$ and the usual mutation kernel M. To simplify the formulas, we suppose that the length $\ell$ is even. The Hamming class number H(u) of $u\in\lbrace\,0,1\,\rbrace^\ell$ is the number of ones present in u, that is,

\begin{align*}H(u)\,=\,\sum_{1\leq i\leq\ell}u(i).\end{align*}

We consider the fitness function $f_{\text{PL}}$ defined as follows.

Hypothesis 3. The fitness function $f_{\textit{PL}}$ is given by

(1.12) \begin{equation}\forall\,u\in\lbrace\,0,1\,\rbrace^\ell{\qquad}f_{\textit{PL}} (u)\,=\,\begin{cases}\quad\sigma&\quad\textit{if }H(u)=\ell/2\,,\\[5pt] \quad1&\quad\textit{otherwise}\,.\end{cases}\end{equation}

The plateau associated with the fitness function $f_{\text{PL}} (u)$ consists of the sequences having $\ell/2$ zeroes and $\ell/2$ ones. In particular, it contains sequences which are very far away in the Hamming distance. This might go against the intuitive notion of a plateau, especially in low-dimensional spaces, where we would usually take for a plateau a closed neighborhood around a fixed connected set. However, our goal is to find a subset of $\lbrace\,0,1\,\rbrace^\ell$ which is very stable with respect to mutations, and the plateau of $f_{\text{PL}} (u)$ is the simplest example of such a set. We could also modify this set so that it is closer to our intuition of a plateau, by adding to it all the sequences which are at Hamming distance less than $\sqrt{\ell}$ from it. All the results presented here and in the subsequent section would hold for this larger set as well (but the proofs would become more complicated).

We consider again the long-chain regime with a finite:

\begin{equation*}\ell\to+\infty\,,{\qquad}q\to 0\,,{\qquad}\ell q\to a\in\,[0,+\infty[.\end{equation*}

Theorem 4. Let $\lambda_{\textit{PL}}(\ell,q)$ be the mean fitness for the model associated with the fitness function $f_{\textit{PL}}$ and the parameters $\ell,q$ . For any $a>0$ , we have

\begin{align*}\lim_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }\lambda_{\textit{PL}}(\ell,q)\,=\,\lambda(a,\sigma)\,,\end{align*}

where $\lambda(a,\sigma)$ is the unique real number $\lambda>1$ such that

(1.13) \begin{equation}\frac{1}{\sigma-1}\,=\,\sum_{n\geq 1}\frac{e^{-an}}{\lambda^n}\sum_{k\geq 0}\frac{(an)^{2k}}{2^{2k}(k!)^2}.\end{equation}

The situation for the fitness landscape $f_{\text{PL}}$ with the huge plateau is very different from the sharp-peak landscape. Indeed, the error threshold phenomenon described in Theorem 1 does not occur, and there is no phase transition in the long-chain regime. It turns out that, if we shift the plateau away from $\ell/2$ , then an error threshold phenomenon reappears. We consider the fitness function $f_{\text{PL}}$ defined as follows.

Hypothesis 4. Let $\sigma>1$ and $\alpha\in [0,1]$ . The fitness function $f_{\textit{PL}}$ is given by

(1.14) \begin{equation}\forall\,u\in\lbrace\,0,1\,\rbrace^\ell{\qquad}f_{\textit{PL}} (u)\,=\,\begin{cases}\quad\sigma&\quad\textit{if }H(u)=\lfloor \alpha\ell\rfloor\,,\\[5pt] \quad1&\quad\textit{otherwise}\,,\end{cases}\end{equation}

where $\lfloor \alpha\ell\rfloor$ is the integer part of $\alpha\ell$ .

The cases $\alpha=0$ and $\alpha=1$ correspond in fact to the sharp-peak landscape, which has already been studied. We shall therefore exclude these cases.

Theorem 5. Suppose that $\alpha\not\in\{\,0,1/2,1\,\}$ . Let $\lambda_{\textit{PL}}(\ell,q)$ be the mean fitness for the model associated with the fitness function $f_{\textit{PL}}$ and the parameters $\ell,q$ . There exist two positive values $a^1_c(\sigma,\alpha)\leq a^2_c(\sigma,\alpha)$ such that the following holds: if $a\geq a^2_c(\sigma,\alpha)$ , then

\begin{align*}\lim_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }\lambda_{\textit{PL}}(\ell,q)\,=\,1.\end{align*}

If $a<a^1_c(\sigma,\alpha)$ , then

\begin{align*}\liminf_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }\lambda_{\textit{PL}}(\ell,q)\,>\,1.\end{align*}

Of course, it would be more satisfactory to know that $a^1_c(\sigma,\alpha)=a^2_c(\sigma,\alpha)$ . This would be a consequence of the monotonicity of the function $\phi_\alpha(a)$ defined in (6.14); unfortunately this fact does not seem so obvious.

The proofs of Theorems 4 and 5 are presented in Section 6. The reason underlying these two results is the following. In the neutral fitness landscape, the most likely Hamming class is already the class $\ell/2$ . By installing a plateau precisely on this Hamming class, we reinforce its stability. When we shift the plateau away from the Hamming class $\ell/2$ , we create competition between two Hamming classes, and this leads to the occurrence of a phase transition.

1.7. Survival of the flattest

We finally prove rigorously within our framework the possible occurrence of the survival of the flattest, a phenomenon previously discovered by Wilke et al. [21], which has been investigated in several works (see e.g. [6, 14, 19, 20]). We now state the precise mathematical result we obtain. Here we will combine the sharp-peak landscape with a plateau. We consider the genotype space of binary sequences $E=\lbrace\,0,1\,\rbrace^\ell$ and the usual mutation kernel M, with the following fitness function. To simplify the formulas, we suppose that the length $\ell$ is even.

Hypothesis 5. Let $\delta,\sigma\geq 1$ and let $f_{\textit{SP/PL}}$ be the fitness function given by

(1.15) \begin{equation}\forall\,u\in\lbrace\,0,1\,\rbrace^\ell{\qquad}f_{\textit{SP/PL}}(u)\,=\,\begin{cases}\quad\delta&\quad\textit{if }u=0\cdots 0\,,\\[5pt] \quad\sigma&\quad\textit{if }u\textit{ contains }\frac{\ell}{2}\textit{ zeroes and }\frac{\ell}{2}\textit{ ones}\,,\\[5pt] \quad1&\quad\textit{otherwise}\,.\end{cases}\end{equation}

Theorem 6. Let $a>0$ , and let $\lambda(a,\sigma)$ be the unique real number $\lambda>1$ such that

(1.16) \begin{equation}\frac{1}{\sigma-1}\,=\,\sum_{n\geq 1}\frac{e^{-an}}{\lambda^n}\sum_{k\geq 0}\frac{(an)^{2k}}{2^{2k}(k!)^2}.\end{equation}

Let $\lambda_{\textit{SP/PL}}(\delta,\sigma,\ell,q)$ be the mean fitness for the model associated with the fitness function $f_{\textit{SP/PL}}$ and the parameters $\ell,q$ . We have the following convergence in the long-chain regime:

\begin{align*}\lim_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }\lambda_{\textit{SP/PL}}(\delta,\sigma,\ell,q)\,=\,\max\big(\delta e^{-a},\lambda(a,\sigma)\big).\end{align*}

Let us denote by y(0) the fraction of the sequence $0\cdots 0$ in the population at equilibrium, and by $y(\ell /2)$ the fraction of the sequences having $\ell/2$ zeroes and $\ell/2$ ones in the population at equilibrium.

• • If $\lambda(a,\sigma)>\delta e^{-a}$ , then $y(0)\to 0$ , $y(\ell/2)\to \frac{\lambda(a,\sigma)-1}{\sigma-1}$ .

• • If $\delta e^{-a}>\lambda(a,\sigma)$ , then $y(\ell/2)\to 0$ , $y(0)\to \frac{\delta e^{-a}-1}{\delta-1}$ .

A plateau is always more stable than a peak of the same height; that is, for any $a>0$ , any $\sigma>1$ , we have $\lambda(a,\sigma)>\sigma e^{-a}$ . However, for a fixed value $a>0$ , there exist positive values $\sigma<\delta$ such that $\lambda(a,\sigma)>\delta e^{-a}$ . Indeed, let us compute the right-hand side of Equation (1.16) with $\lambda$ replaced by $\delta e^{-a}$ :

\begin{align*}\sum_{n\geq 1}\frac{e^{-an}}{\big(\delta e^{-a}\big)^n}\sum_{k\geq 0}\frac{(an)^{2k}}{2^{2k}(k!)^2}&\,=\,\sum_{n\geq 1}\frac{1}{\delta^n}\sum_{k\geq 0}\frac{(an)^{2k}}{2^{2k}(k!)^2}\cr &\,\geq \,\sum_{n\geq 1}\frac{1}{\delta^n}\Big(1+ \frac{(an)^{2}}{4}\Big)\,=\,\frac{1}{\delta-1}+\frac{a^{2}}{4}\frac{\delta(\delta+1)}{(\delta-1)^3}.\end{align*}

So, for any value of $\delta>\sigma$ such that

\begin{align*}\frac{1}{\sigma-1}\,<\,\frac{1}{\delta-1}+\frac{a^{2}}{4}\frac{\delta(\delta+1)}{(\delta-1)^3}\,,\end{align*}

we will indeed have $\lambda(a,\sigma)>\delta e^{-a}$ . For such values, the quasispecies will be located on the plateau, despite the fact that the height $\sigma$ of the plateau is strictly less than the height of the peak. In this situation, we witness the survival of the flattest.

Theorem 6 is proved in Section 7. The route to proving Theorems 3 and 6 is quite long. It rests on the exact formulas proved in Section 4.2 and the general strategy developed in Section 4.

2. The sharp-peak landscape

Our first goal is to develop an exact formula for the quasispecies on the sharp-peak landscape. All our subsequent study of the sharp-peak landscape will rely on this formula.

2.1. A representation formula for the quasispecies

The computations performed in this subsection are in fact valid without any assumption on the geometry of the genotype space E. We merely suppose that E is finite and that the coefficients of the matrix M are all positive. Suppose that $(x(u), u\in E)$ is a solution of the system (1.1) satisfying the constraint (1.2). Let $\lambda$ be the mean fitness of this solution, defined by

(2.1) \begin{equation}\lambda\,=\,\sum_{v\in E}x(v)f(v)\,=\,\sigma x(w^*)+1-x(w^*).\end{equation}

Obviously, the mean fitness $\lambda$ satisfies $1\leq\lambda\leq\sigma$ and $\lambda=1$ if and only if $x(w^*)=0$ . By hypothesis, all the coefficients of the mutation matrix M are positive, and it follows from Equation (1.1) that

(2.2) \begin{equation}\forall\, u\in E{\qquad} x(u)\,\sigma\,\geq\,\,\sum_{v\in E}x(v)f(v)M(v,u)\,>\,0.\end{equation}

Therefore $x(u)>0$ for all $u\in E$ . Equation (1.1) for a generic genotype u reads as follows:

(2.3) \begin{equation}\lambda x(u)\,=\,{\sigma} x(w^*)M(w^*,u)+\sum_{v\neq w^*}x(v)M(v,u).\end{equation}

We rewrite this equation as

(2.4) \begin{equation}x(u)\,=\,\frac{\sigma-1}{\lambda} x(w^*)M(w^*,u)+\frac{1}{\lambda}\sum_{v\in E}x(v)M(v,u).\end{equation}

Using this very formula for x(v), we replace x(v) in the sum, thereby getting

\begin{multline*}x(u)\,=\,\frac{\sigma-1}{\lambda} x(w^*)M(w^*,u)+\frac{\sigma-1}{\lambda^2}\sum_{v\in E} x(w^*) M(w^*,v) M(v,u)\cr+ \frac{1}{\lambda^2} \sum_{v,v'\in E}x(v')M(v',v)M(v,u)\,.\end{multline*}

We iterate this process. After $N-1$ steps, we get

\begin{multline*}x(u)\,=\,\sum_{1\leq n\leq N}\frac{\sigma-1}{\lambda^n}\sum_{v_1,\dots,v_{n-1}\in E}\kern-7ptx(w^*)M(w^*,v_1)\cdots M(v_{n-1},u)\cr\,+\,\frac{1}{\lambda^N}\sum_{v_1,\dots,v_{N}\in E} \kern-7ptx(v_1)M(v_1,v_2)\cdots M(v_{N},u)\,.\end{multline*}

The sums of products involving the matrix M have a natural probabilistic interpretation. Let $(X_n)_{n\in\mathbb N}$ be the Markov chain on E with transition matrix M. The previous formula for x(u) can now be rewritten as

\begin{multline*}x(u)\,=\,\sum_{1\leq n\leq N}\frac{\sigma-1}{\lambda^n}x(w^*)\mathbb{P}\big(X_{n}=u\,\big|\,X_0=w^*\big)\cr\,+\,\frac{1}{\lambda^N}\sum_{v_1\in E}x(v_1)\mathbb{P}\big(X_{N}=u\,\big|\,X_0=v_1\big)\,.\end{multline*}

By (2.1) and (2.2), the mean fitness $\lambda$ is strictly larger than 1; therefore the last term goes to 0 as N goes to $\infty$ . Sending N to $\infty$ , and using the fact that

(2.5) \begin{equation}\lambda-1\,=\,(\sigma-1)x(w^*),\end{equation}

we get the following result.

Proposition 1. Any solution $(x(u),u\in E)$ of the quasispecies equations (1.1) for the sharp-peak function $f_{SP}$ defined in (1.4) satisfies

(2.6) \begin{equation}x(u)\,=\,\sum_{n\geq 1}\frac{\lambda-1}{\lambda^n}\mathbb{P}\big(X_{n}=u\,\big|\,X_0=w^*\big),\end{equation}

where $(X_n)_{n\geq0}$ is the Markov chain on E with transition matrix M.

Let us take $u=w^*$ in the formula (2.6). Using (2.5) again, we obtain

(2.7) \begin{equation}\frac{1}{\sigma-1}\,=\,\sum_{n\geq 1}\frac{1}{\lambda^n}\mathbb{P}\big(X_{n}=w^*\,\big|\,X_0=w^*\big).\end{equation}

The right-hand side of this equation is a continuous decreasing function of $\lambda$ , which is equal to $+\infty$ when $\lambda=1$ and is less than or equal to $1/(\sigma-1)$ when $\lambda=\sigma$ . Therefore there exists exactly one value of $\lambda$ in $[1,\sigma]$ which satisfies this equation. Our next goal is to obtain an explicit expression for $P\big(X_{n}=w^*\,\big|\,X_0=w^*\big)$ and to perform the summation of the series appearing in (2.7). To do so, we shall make additional assumptions on the geometry of the fitness landscape.

2.2. Mutation dynamics

We now turn our attention to the genotype space $E=\mathcal{A}^\ell$ where $\mathcal{A}$ is a finite alphabet of cardinality $\kappa$ , and we consider the mutation scheme (1.3). In this model, the mutations occur independently at each site. An important consequence of this structural assumption is that the components $(X_n(i),1\leq i\leq\ell)$ of the Markov chain $X_n$ are themselves independent Markov chains with state space $\mathcal{A}$ and transition matrix

\begin{align*}\begin{pmatrix} 1-q\;\;\;\; & \frac{q}{\kappa-1}\;\;\;\; &\cdots\;\;\;\; & \frac{q}{\kappa-1}\cr \frac{q}{\kappa-1}\;\;\;\; & 1-q\;\;\;\; &\ddots\;\;\;\; & \frac{q}{\kappa-1}\cr\vdots\;\;\;\; &\ddots\;\;\;\; & \ddots\;\;\;\; &\vdots\cr \frac{q}{\kappa-1}\;\;\;\;& \cdots\;\;\;\; & \frac{q}{\kappa-1}\;\;\;\; & 1-q\cr\end{pmatrix}.\end{align*}

The non-diagonal terms in this matrix are all equal to $q/(\kappa-1)$ . Since we want to compute $\mathbb{P}\big(X_{n}=w^*\,\big|\,X_0=w^*\big)$ , we shall lump together the letters which differ from the wild type $w^*$ , and we shall do this for each component. More precisely, we define a process $(V_n)_{n\geq 0}$ by setting

\begin{align*}\forall i\in\{\,1,\dots,\ell\,\}{\qquad} V_n(i)\,=\,\begin{cases}0 &\text{if $X_n(i)=w^*(i)$}\,,\\[5pt] 1 &\text{if $X_n(i)\neq w^*(i)$}\,. \end{cases}\end{align*}

The binary word $V_n$ indicates the sites where $X_n$ and $w^*$ differ. In particular, $V_n$ is a deterministic function of $X_n$ . Thanks to the specific form of the transition matrix, it turns out that $V_n(i)$ is still a Markov chain. This is a particular case of the lumping theorem of Kemeny and Snell (see Theorem $6.3.2$ of [15]). To see why, consider for instance the first two steps of the first component $V_n(1)$ , and let us make the following computation:

\begin{multline*}\mathbb{P}(V_2(1)=0,\,V_1(1)=1\,\big|\,V_0(1)=0\big)\cr\,=\,\sum_{x_1\in \mathcal{A}\setminus \{\,w^*(1)\,\} }\mathbb{P}(X_2(1)=w^*(1),\,X_1(1)=x_1\,\big|\,V_0(1)=0\big)\cr\,=\,\sum_{x_1\in \mathcal{A}\setminus \{\,w^*(1)\,\} }\mathbb{P}(X_2(1)=w^*(1)\,|\,X_1(1)=x_1)\cr\times\mathbb{P}( X_1(1)=x_1 \,|\,X_0(1)=w^*(1))\,=\, (\kappa-1)\times\frac{q}{\kappa-1}\times \frac{q}{\kappa-1}\,=\, {q}\frac{q}{\kappa-1}\,.\end{multline*}

The simplification in the final result comes from the fact that the probabilities of finding or losing a letter from the wild type $w^*$ are the same for all the letters in $\mathcal{A}\setminus \{\,w^*(1)\,\}$ . A similar computation can be made for a finite number of steps. We conclude that $V_n(i)$ is the two-state Markov chain that we define and study next. Let $(E_n)_{n\geq 0}$ be the Markov chain with state space $\{0,1\}$ and transition matrix

\begin{align*}T\,=\,\begin{pmatrix} 1-q\;\;\;\; & q \cr \frac{q}{\kappa-1}\;\;\;\; & 1- \frac{q}{\kappa-1}\cr\end{pmatrix}.\end{align*}

The eigenvalues of T are 1 and

(2.8) \begin{equation}\rho\,=\,1- \frac{\kappa q}{\kappa-1}.\end{equation}

We will study some asymptotic regimes in which q is sent to 0, so we can suppose that $q\leq 1-1/\kappa$ (otherwise $\rho$ might be negative, which would create additional complications). It is a standard exercise to compute the powers of T:

\begin{align*}\forall n\geq 1{\qquad} T^n\,=\,\begin{pmatrix} \frac{1}{\kappa} + \frac{\kappa-1}{\kappa}\rho^n\;\;\;\; & \frac{\kappa-1}{\kappa} \big(1-\rho^n\big) \\[10pt] \frac{1}{\kappa}\big(1-\rho^n \big)\;\;\;\; & \frac{\kappa-1}{\kappa}+ \frac{1}{\kappa}\rho^n \cr\end{pmatrix}.\end{align*}

Here is a simple, illuminating way to realize the Markov chain $(E_n)_{n\geq 0}$ and to understand the expression for the nth power $T^n$ . Let $(\varepsilon_n)_{n\geq 1}$ be a sequence of independent and identically distributed Bernoulli random variables with parameter $\rho$ . Suppose that $E_{n-1}=e\in \{\,0,1\,\}$ . If $\varepsilon_n=1$ , then we set $E_n=e$ . If $\varepsilon_n=0$ , then we choose a letter uniformly over $\mathcal{A}$ , independently of the past history up to time n, and we set $E_n=0$ if the chosen letter is the one of the wild type and $E_n=1$ otherwise. Now, the event $E_n=E_0$ can occur in two different ways. Either $\varepsilon_1=\cdots=\varepsilon_n=1$ , or one of $\varepsilon_1,\dots,\varepsilon_n$ is zero, in which case $E_n=0$ with probability $1/\kappa$ and $E_n=1$ with probability $(\kappa-1)/\kappa$ ; thus

\begin{align*}\mathbb{P}(E_n=0\,|\,E_0=0)\,&=\,\rho^n+\big(1-\rho^n\big)\frac{1}{\kappa}\,,\cr\mathbb{P}(E_n=1\,|\,E_0=1)\,&=\,\rho^n+\big(1-\rho^n\big)\frac{\kappa-1}{\kappa}\,,\end{align*}

and we recover the expression for the diagonal coefficients of $T^n$ . In words, the status $E_n$ at step n is the same as at time 0 if no mutation has occurred, or if the last mutation results in the same letter as the wild type (case $E_0=0$ ) or in a different letter (case $E_0=1$ ). Similarly, the event $E_n\neq E_0$ can occur only if one of $\varepsilon_1,\dots,\varepsilon_n$ is zero, and the last mutation event yields the adequate letter; thus

\begin{align*}\mathbb{P}(E_n=1\,|\,E_0=0)\,&=\,\big(1-\rho^n\big)\frac{\kappa-1}{\kappa}\,,\cr\mathbb{P}(E_n=0\,|\,E_0=1)\,&=\,\big(1-\rho^n\big)\frac{1}{\kappa}.\end{align*}

Now the probability $\mathbb{P}\big(X_{n}=w^*\,\big|\,X_0=w^*\big)$ can be rewritten with the help of $V_n$ and $E_n$ as

(2.9) \begin{align}\mathbb{P}\big(X_{n}=w^*&\,\big|\,X_0=w^*\big)\,=\,\mathbb{P}\big(X_{n}(i)=w^*(i), \,1\leq i\leq\ell\,\big|\,X_0=w^*\big)\cr &\,=\,\mathbb{P}\big(V_{n}(i)=0, \,1\leq i\leq\ell\,\big|\,X_0=w^*\big)\cr &\,=\,{\prod_{1\leq i\leq\ell}} \mathbb{P}\big(V_{n}(i)=0\,\big|\,X_0=w^*\big)\cr &\,=\,\Big(\mathbb{P}\big(E_{n}=0\,\big|\,E_0=0\big)\Big)^\ell\,=\,\Big( \frac{1}{\kappa} + \frac{\kappa-1}{\kappa}\rho^n \Big)^\ell.\end{align}

2.3. The equation for the mean fitness

We now have all the material necessary to obtain Equation (1.5). Substituting the expression obtained in (2.9) into the formula (2.7), we get

(2.10) \begin{equation}\frac{1}{\sigma-1}\,=\,\sum_{n\geq 1}\frac{1}{\lambda^n}\Big( \frac{1}{\kappa} + \frac{\kappa-1}{\kappa}\rho^n \Big)^\ell.\end{equation}

We expand the power in $\ell$ , exchange the summations, and re-sum the series as follows:

\begin{align*}\frac{1}{\sigma-1}&\,=\,\sum_{n\geq 1}\frac{1}{\lambda^n}\sum_{0\leq k\leq \ell}\binom{\ell}{k}\Big( \frac{1}{\kappa}\Big)^{\ell-k}\Big( \frac{\kappa-1}{\kappa}\rho^n \Big)^{k}\cr&\,=\, \frac{1}{\kappa^\ell}\sum_{0\leq k\leq \ell}\binom{\ell}{k}(\kappa-1)^k\sum_{n\geq 1} \frac{ \rho^{kn}}{\lambda^n}\cr& \,=\, \frac{1}{\kappa^\ell}\sum_{0\leq k\leq \ell}\binom{\ell}{k}(\kappa-1)^k\frac{1}{\frac{\lambda}{\rho^{k}}-1}.\end{align*}

Thus the equation satisfied by the mean fitness $\lambda$ reads

(2.11) \begin{equation}\frac{1}{\sigma-1}\,=\, \frac{1}{\kappa^\ell}\sum_{0\leq k\leq \ell}\binom{\ell}{k}(\kappa-1)^k\frac{1}{\frac{\lambda}{\rho^{k}}-1}.\end{equation}

This equation was discovered by Bratus et al. (see the formulas (5.5) and (5.6) in [3]), who exploited it to derive rigorous bounds on the mean fitness. This formula calls naturally for a probabilistic interpretation. Let $S_\ell$ be a random variable with distribution the binomial law with parameters $\ell$ and $(\kappa-1)/\kappa$ . Denoting the expectation by $\mathbb{E}$ , we can rewrite Equation (2.11) as

(2.12) \begin{equation}\frac{1}{\sigma-1}\,=\, \mathbb{E}\left(\frac{1}{\frac{\lambda}{\rho^{S_\ell}}-1}\right).\end{equation}

Our analysis of the quasispecies equations on the sharp-peak landscape rests on the identity (2.10) or its equivalent form (2.12). Indeed, if we manage to estimate the mean fitness $\lambda$ , then we will also have an estimate for the proportion $x(w^*)$ of the wild type present in the population. Moreover, once we know $\lambda$ or $x(w^*)$ , the proportions of the other types are completely determined by the formula (2.6).

2.4. An exact formula for the quasispecies

Loosely speaking, a quasispecies is a cloud of mutants centered around a specific genotype. The structure of this cloud depends on the parameters $\sigma$ and q. In fact, the proportion of the wild type $w^*$ decreases as the mutation rate q increases, and it becomes comparable to the proportions of the other genotypes when q reaches the error threshold $q^*$ . This fascinating phenomenon will be established rigorously in Section 3.2. When q is a little below $q^*$ , we observe a cloud of mutants centered around $w^*$ which contains a very small proportion of the wild type $w^*$ ; the vast majority of the genotypes present in the population will differ from $w^*$ . Yet the genetic information carried out by $w^*$ is still present in the population and determines its structure. This paradoxical situation has led several biologists to argue that selection operates at the level of the quasispecies, and not at the level of individuals [8]. Thus an important goal is to understand better the statistical composition of the cloud of mutants. In order to do so, we classify the genotypes according to their Hamming distance from the wild type $w^*$ . For $k\geq 1$ , we define the Hamming class k as the set $\mathcal{H}_k$ of the genotypes which differ from $w^*$ at exactly k indices, i.e.,

(2.13) \begin{equation}\mathcal{H}_k\,=\,\big\{\,u\in E \;:\;d(u,w^*)=k\,\big\}.\end{equation}

We shall exploit further the formula (2.6) in order to derive an exact formula for the proportion ${\mathcal{Q}(k)}$ of the genotypes belonging to the Hamming class k. For $k\geq 1$ , we have, thanks to the formula (2.6),

(2.14) \begin{align}{\mathcal{Q}}(k)&\,=\,\sum_{u\in\mathcal{H}_k}\sum_{n\geq 1}\frac{\lambda-1}{\lambda^n}\mathbb{P}\big(X_{n}=u\,\big|\,X_0=w^*\big)\,\cr&\,=\,\sum_{n\geq 1}\frac{\lambda-1}{\lambda^n}\sum_{u\in\mathcal{H}_k}\mathbb{P}\big(X_{n}=u\,\big|\,X_0=w^*\big)\,\cr&\,=\,\sum_{n\geq 1}\frac{\lambda-1}{\lambda^n}\mathbb{P}\big(X_{n} \in\mathcal{H}_k\,\big|\,X_0=w^*\big).\end{align}

We have already noticed that the components of $X_n$ , $(X_n(i),1\leq i\leq\ell)$ , are independent Markov chains. Therefore, using the notation of Section 2.2, we have

(2.15) \begin{align}\mathbb{P}\big(X_{n} \in\mathcal{H}_k&\,\big|\,X_0=w^*\big)\,=\,\mathbb{P}\Big(\sum_{1\leq i\leq\ell}V_{n}(i)=k\,\big|\,X_0=w^*\Big)\cr&\,=\,\binom{\ell}{k}\Big(\mathbb{P}\big(E_{n}=1\,\big|\,E_0=0\big)\Big)^k\Big(\mathbb{P}\big(E_{n}=0\,\big|\,E_0=0\big)\Big)^{\ell-k}\cr&\,=\,\binom{\ell}{k}\Bigg( \big(1-\rho^n\big)\frac{\kappa-1}{\kappa} \Bigg)^k\Bigg(\rho^n+\big(1-\rho^n\big)\frac{1}{\kappa}\Bigg)^{\ell-k}.\end{align}

Recall that $\rho=1-{\kappa q}/({\kappa-1})$ . Plugging (2.15) into (2.14), we get an exact formula for the quasispecies distribution, in terms of a series involving $\lambda$ and $\rho$ :

(2.16) \begin{equation}{\mathcal{Q}}(k)\,=\,\sum_{n\geq 1}\frac{\lambda-1}{\lambda^n}\binom{\ell}{k}\Bigg( \big(1-\rho^n\big)\frac{\kappa-1}{\kappa} \Bigg)^k\Bigg(\rho^n+\big(1-\rho^n\big)\frac{1}{\kappa}\Bigg)^{\ell-k}.\end{equation}

We can even compute the sum of the series, by developing the powers, and obtain a finite algebraic formula:

(2.17) \begin{equation}{\mathcal{Q}}(k)\,=\,(\lambda-1)\binom{\ell}{k}\frac{(\kappa-1)^k}{\kappa^\ell}\sum_{i=0}^k\sum_{j=0}^{\ell-k}\binom{k}{i}\binom{\ell-k}{j}\frac{(-1)^i(\kappa-1)^j\rho^{i+j}}{\lambda-\rho^{i+j}}.\end{equation}

Of course this formula is complicated, yet it expresses completely the dependence of ${\mathcal{Q}}(k)$ on $\lambda$ and q, and it shows the complexity of the sharp-peak landscape. Surprisingly, the mean and the variance of the Hamming distance of the quasispecies distribution have very simple expressions. Indeed, the mean of the Hamming distance of the distribution $\overline{\mathcal{Q}}$ is given by

(2.18) \begin{equation}\overline{\mathcal{Q}} =\sum_{k=0}^\ell k \,{\mathcal{Q}}(k) =\frac{\ell q\, \lambda}{\lambda-\rho}.\end{equation}

To compute this mean, we rely on the formula (2.16) to write

(2.19) \begin{align}\overline{\mathcal{Q}}&\, = \,\sum_{k=0}^\ell k \,{\mathcal{Q}}(k)\cr &\,= \,\sum_{n\geq 1}\frac{\lambda-1}{\lambda^n}\sum_{k=0}^\ell k\binom{\ell}{k}\Bigg( \big(1-\rho^n\big)\frac{\kappa-1}{\kappa} \Bigg)^k\Bigg(\rho^n+\big(1-\rho^n\big)\frac{1}{\kappa}\Bigg)^{\ell-k}.\end{align}

In the inner sum we recognize the expectation of a binomial law, whence

\begin{align*}\overline{\mathcal{Q}} \,=\, \sum_{n\geq 1}\frac{ \lambda-1 }{\lambda^n} \ell \big(1-\rho^n\big)\frac{\kappa-1}{\kappa} .\end{align*}

By summing the geometric series and using the definition of $\rho$ , we obtain the formula (2.18). In fact, all the successive moments of the quasispecies distribution can be computed in this way.

3. The error threshold

We continue to consider the case of the sharp-peak landscape described in Section 1.2. In Section 2, we performed exact computations with the length $\ell$ fixed. Here we shall perform asymptotic expansions in the long-chain regime. We will revisit the classical error threshold phenomenon and compute the limiting quasispecies distribution.

3.1. An upper bound on $\lambda$

The goal of this subsection is to prepare the ground for the technical proof of the error threshold result. We start from the identity (2.12), and we wish to obtain an upper bound on $\lambda$ . From the weak law of large numbers, we know that the most likely values for $S_\ell$ are those around its mean $(\kappa-1)\ell/\kappa$ , so we pick a positive number $\varepsilon$ and split the expectation as follows:

(3.1) \begin{align}\frac{1}{\sigma-1}&=\sum_{0\leq k\leq \ell}\mathbb{P}(S_\ell=k)\frac{1}{\frac{\lambda}{\rho^{k}}-1}\nonumber \\[5pt]&=\sum_{k:|k-(\kappa-1)\ell/\kappa|>\varepsilon\ell}\mathbb{P}(S_\ell=k)\frac{1}{\frac{\lambda}{\rho^{k}}-1}+\sum_{k:|k-(\kappa-1)\ell/\kappa|\leq\varepsilon\ell}\mathbb{P}(S_\ell=k)\frac{1}{\frac{\lambda}{\rho^{k}}-1}\,.\end{align}

Taking advantage of the fact that $\rho\leq 1$ , we bound $\rho^k$ by 1 in the first sum and by $\rho^{(\kappa-1)\ell/\kappa-\varepsilon\ell}$ in the second; this leads to the inequality

(3.2) \begin{equation}\frac{1}{\sigma-1}\,\leq\,\frac{1}{ {\lambda}-1}\mathbb{P}\Big(\Big|S_\ell-\frac{\kappa-1}{\kappa}\ell\Big|>\varepsilon\ell\Big)+\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)\ell/\kappa-\varepsilon\ell} }-1}.\end{equation}

We apply the Chebyshev inequality to the binomial random variable $S_\ell$ to obtain the classical estimate

(3.3) \begin{equation}\mathbb{P}\Big(\Big|S_\ell-\frac{\kappa-1}{\kappa}\ell\Big|>\varepsilon\ell\Big)\,\leq\, \frac{\text{variance}(S_\ell)}{\varepsilon^2\ell^2}\,=\, \frac{\kappa-1}{\kappa^2\varepsilon^2\ell}.\end{equation}

Plugging the Chebyshev inequality into (3.2), we conclude that

(3.4) \begin{equation}\forall\varepsilon>0{\qquad}\frac{1}{\sigma-1}\,\leq\,\frac{\kappa-1}{\kappa^2\varepsilon^2\ell(\lambda-1)}+\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)\ell/\kappa-\varepsilon\ell} }-1}.\end{equation}

Let us pause for one moment to look at this inequality. Obviously, as $\lambda$ goes to $\infty$ , the right-hand side goes to 0, and the inequality cannot hold. In other words, for this inequality to hold, $\lambda$ must not be too large. Thus, from this inequality, we should be able to derive an upper bound on $\lambda$ . To get the best possible upper bound, we could rewrite it as a polynomial inequality of degree two in $\lambda$ and compute the associated roots, but this leads to messy expressions. Instead, after a few trials (hidden from the reader), we obtain the following inequality.

Lemma 1. The mean fitness $\lambda$ satisfies

(3.5) \begin{equation}\lambda\,\leq\,\max\Bigg(1+\frac{4\sigma}{\ell^{1/4}},\sigma \Big(1- \frac{\kappa q}{\kappa-1}\Big)^{{(\kappa-1)\ell}/{\kappa}-\ell^{2/3}}\Big(1+\frac{1}{\ell^{1/13}}\Big)\Bigg).\end{equation}

Proof. Suppose that $\lambda$ is larger than the upper bound in (3.5). We check that this is not compatible with the inequality (3.4) with the choice $\varepsilon={\ell}^{-1/3}$ . Indeed, we have on the one hand

\begin{align*}\lambda>1+\frac{4\sigma}{\ell^{1/4}}{\qquad}\Longrightarrow{\qquad}\frac{1}{\ell^{1/3}(\lambda-1)}\,\leq\,\frac{1}{4\sigma\ell^{1/12}}.\end{align*}

On the other hand,

\begin{multline*}\lambda>\sigma \Big(1- \frac{\kappa q}{\kappa-1}\Big)^{{(\kappa-1)\ell}/{\kappa}-\ell^{2/3}}\Big(1+\frac{1}{\ell^{1/13}}\Big){\qquad}\Longrightarrow{\qquad}\cr\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)\ell/\kappa-\ell^{2/3}} }-1}\,\leq\,\frac{1}{\sigma\Big(1+\frac{1}{\ell^{1/13}}\Big)-1}\,.\end{multline*}

Applying the mean value theorem to the function $x\mapsto 1/(\sigma-1+x)$ on the interval $[0,\sigma/ {\ell^{1/13}}]$ , we have

(3.6) \begin{equation}\frac{1}{ \sigma\Big(1+\frac{1}{\ell^{1/13}}\Big)-1}\,\leq\,\frac{1}{ \sigma -1}-\frac{1}{ ( 2\sigma -1)^2}\frac{\sigma} {\ell^{1/13}}\,\leq\,\frac{1}{ \sigma -1}-\frac{1} {4\sigma\ell^{1/13}}.\end{equation}

Combining the previous inequalities, we obtain that

\begin{equation*}\frac{\kappa-1}{\kappa^2\ell^{1/3}(\lambda-1)}\,+\,\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)\ell/\kappa-\ell^{2/3}} }-1}\,\leq\,\frac{1}{ \sigma -1}+\frac{1} {4\sigma}\Big(\frac{1}{\ell^{1/12}}-\frac{1}{\ell^{1/13}}\Big)\,<\,\frac{1}{ \sigma -1},\end{equation*}

which stands in contradiction to the inequality (3.4) with $\varepsilon={\ell}^{-1/3}$ .

3.2. Revisiting the error threshold

Because of its importance, we shall give two proofs of Theorem 1. The first proof is a soft proof based on a compactness argument, whose starting point is the initial quasispecies equation. The second proof is a more technical proof, whose starting point is the Equation (2.10) discovered by Bratus et al. Naturally, the second proof is more informative: it yields a control on the speed of convergence and opens the way to performing an asymptotic expansion of $\lambda$ with respect to $\ell$ and q.

From Equation (2.1), we can express $x(w^*)$ as

(3.7) \begin{equation}x(w^*)\,=\,\frac{\lambda-1}{\sigma-1}.\end{equation}

Therefore it is enough to study the asymptotic behavior of $\lambda$ . We shall prove that

(3.8) \begin{equation}\lambda\quad\longrightarrow\quad \max\big(1, \,\sigma e^{-a}\big).\end{equation}

Soft proof. In the long-chain regime, we have

(3.9) \begin{equation}\sigma \big(1-q\big)^\ell\quad\longrightarrow\quad \sigma e^{-a}.\end{equation}

The quasispecies equation (1.1) associated with $w^*$ reads

(3.10) \begin{equation}\lambda x(w^*)\,=\,\sigma x(w^*)M(w^*,w^*)\,+\,\sum_{v\in E\setminus \{w^*\}}x(v)M(v,w^*).\end{equation}

We know that $x(w^*)>0$ ; therefore we deduce from (3.10) that

(3.11) \begin{equation}\lambda\,\geq\,\sigma M(w^*,w^*)\,=\,\sigma (1-q)^\ell.\end{equation}

Using the lower bound (3.11) and the fact that $\lambda\geq 1$ , we see that

(3.12) \begin{equation}\liminf\,\,\lambda\,\geq\,\max\big(1, \,\sigma e^{-a}\big).\end{equation}

Suppose that $\lambda$ does not converge towards $\max\big(1, \,\sigma e^{-a}\big)$ . In this case, there exist $\varepsilon>0$ and two sequences of parameters $(\ell_n)_{n\geq 0}$ , $(q_n)_{n\geq 0}$ such that

\begin{equation*}\ell_n\to+\infty\,,{\qquad}q_n\to 0\,,{\qquad}\ell_n q_n\to a,\end{equation*}

(3.13) \begin{equation}\forall n\geq 0{\qquad}\lambda(\ell_n,q_n)\,\geq\,\max\big(1, \,\sigma e^{-a}\big)+\varepsilon.\end{equation}

Using (3.7) and (3.13), we obtain a lower bound on $x(w^*)$ :

(3.14) \begin{equation}x(w^*)\,\geq\,\frac{\varepsilon}{\sigma-1}.\end{equation}

We bound from above the quasispecies equation (3.10) associated with $w^*$ , as follows:

(3.15) \begin{equation}\lambda x(w^*)\,\leq\,\sigma x(w^*)(1-q)^\ell+\,\big(1-x(w^*)\big)q,\end{equation}

where we have used the fact that $M(v,w^*)\leq q$ whenever $v\neq w^*$ . Using (3.14) and (3.15) together, we obtain that

\begin{equation*}\forall n\geq 0{\qquad}\lambda(\ell_n,q_n)\,\leq\,\sigma (1-q_n)^{\ell_n}+\frac{\sigma-1}{\varepsilon}q_n.\end{equation*}

Sending n to $\infty$ , we see that this inequality is not coherent with the inequality (3.13). Therefore it must be the case that $\lambda$ converges towards $\max\big(1, \,\sigma e^{-a}\big)$ .

Technical proof. This proof is entirely based on the equation (2.10) discovered by Bratus et al. (see the formulas $5.5$ and $5.6$ in [3]) and its probabilistic interpretation (2.12). Let us start with a simple yet interesting inequality on $\lambda$ . The map

\begin{align*}\phi\;:\;x\in [0,+\infty[\,\,\mapsto\,\frac{1}{\frac{\lambda}{\rho^{x}}-1}\,\end{align*}

is convex (recall that we suppose that $\rho>0$ ). The equation (2.10) (or rather (2.12)) can be rewritten with the help of the function $\phi$ as

(3.16) \begin{equation}\frac{1}{\sigma-1}\,=\, \mathbb{E}\left(\phi( {S_\ell}) \right),\end{equation}

where $S_\ell$ is a random variable with distribution the binomial law with parameters $\ell$ and $(\kappa-1)/\kappa$ . By the classical Jensen inequality, we therefore have

(3.17) \begin{equation}\frac{1}{\sigma-1}\,\geq\,\phi\big(\mathbb{E}( {S_\ell})\big)\,=\, \phi\Big(\frac{\kappa-1}{\kappa}\ell\Big).\end{equation}

This yields the inequality

(3.18) \begin{equation}\lambda\,\geq\,\sigma\rho^{ {(\kappa-1)\ell}/{\kappa}}\,=\,\sigma\Big(1- \frac{\kappa q}{\kappa-1}\Big)^{ {(\kappa-1)\ell}/{\kappa}}.\end{equation}

Notice that this inequality is not better than the lower bound (3.11) used in the soft proof. However, it can certainly be improved with additional work. In the long-chain regime, we have

(3.19) \begin{equation}\sigma \Big(1- \frac{\kappa q}{\kappa-1}\Big)^{{(\kappa-1)\ell}/{\kappa}}\quad\longrightarrow\quad \sigma e^{-a}.\end{equation}

We combine the lower bound (3.18) and the upper bound (3.5) (which were both derived from the equation (2.10)) with the limit stated in (3.19) to obtain the desired conclusion.

3.3. The limiting quasispecies distribution

In Section 2.4, we computed an exact expression for the mean Hamming distance under the quasispecies distribution, which was

(3.20) \begin{equation}\overline{\mathcal{Q}} \,=\,\sum_{k=0}^\ell k \,{\mathcal{Q}({\sigma},a)}(k) \,=\,\frac{\ell q\, \lambda}{\lambda-\rho}\,.\end{equation}

This formula is particularly illuminating for the error threshold phenomenon. We consider again the long-chain regime

(3.21) \begin{equation}\ell\to+\infty\,,{\qquad}q\to 0\,,{\qquad}\ell q\to a\in\,[0,+\infty].\end{equation}

We have the following dichotomy:

• • If $a>\ln \sigma$ , then $\lambda \rightarrow 1$ , $\lambda-\rho\to 0$ , whence $\overline{\mathcal{Q}}\rightarrow +\infty$ .

• • If $a<\ln \sigma$ , then $\lambda\rightarrow \sigma e^{-a}>1$ , and

  (3.22) \begin{equation}\overline{\mathcal{Q}}\quad\longrightarrow\quad\frac{a\sigma e^{-a}}{\sigma e^{-a}-1}\,.\end{equation}

We can go even further; indeed we can perform the expansion of the formula (2.15) in the regime (3.21). From this we obtain

\begin{align*}\mathbb{P}\big(X_{n}\in\mathcal{H}_k\,\big|\,X_0=w^*\big)\,\sim\,\frac{\ell^k}{k!}(nq)^k(1-nq)^\ell\,\sim\,\frac{(na)^k}{k!}e^{-na}\,;\end{align*}

thus the asymptotic distribution of $X_n$ is the Poisson distribution with parameter na. Combined with the formula (3.8), this yields that, for any fixed $n\geq 1$ ,

\begin{align*}\frac{1}{\lambda^n}\mathbb{P}\big(X_{n}\in\mathcal{H}_k\,\big|\,X_0=w^*\big)\,\sim\,\frac{(na)^k}{\sigma^nk!}\,.\end{align*}

Through brute substitution of this expansion into the formula (2.14), and again applying (3.8), we finally get

\begin{align*}{\mathcal{Q}}(k)\,\sim\,(\sigma e^{-a}-1)\frac{a^k}{k!}\sum_{n\geq 1}\frac{n^k}{\sigma^n}.\end{align*}

Some steps in this computation must be made rigorous, but this can be done by a routine application of the dominated convergence theorem. In fact, this formula was obtained in [5] while studying the quasispecies arising in a stochastic model for the evolution of a finite population. Thanks to the formula, we can easily plot the quasispecies distribution (Figures 1 and 2) and study its dependence on the parameters $\sigma,a$ . In the quasispecies literature, this was previously done by numerically integrating the differential system (see for instance [9, Figure 3, p. 9]). In addition, we can now perform simple theoretical computations with the help of this distribution. For instance, the mutation rate at which the proportion of the Hamming class 1 becomes larger than the proportion of the wild type $w^*$ is $a=1-1/{\sigma}$ .

Figure 1. Fractions of the different types as functions of a for ${\sigma}=5$ .

Figure 2. Fractions of the different types as functions of a for ${\sigma}=10^6$ .

4. Extension to a general landscape

Our starting point is the quasispecies equation (1.1), rewritten in the following form:

(4.1) \begin{equation}\forall\, u\in E{\qquad}\lambda x(u) \,=\,\,\sum_{v\in E}x(v)f(v)M(v,u)\,\end{equation}

subject to the constraint

(4.2) \begin{equation}\forall\, u\in E{\qquad} x(u)\geq 0\,,{\qquad}\sum_{u\in E}x(u)\,=\,1,\end{equation}

where

\begin{align*}\lambda\,=\,\sum_{v\in E}x(v)f(v)\end{align*}

is the mean fitness, or equivalently the Perron–Frobenius eigenvalue of the matrix $(f(v)M(v,u))_{u,v\in E}$ .

4.1. Reformulation of the quasispecies equations

The computations performed in this subsection are in fact valid without any assumption on the geometry of the genotype space E. Thanks to the hypothesis on the fitness function, we have that $\lambda\geq 1$ . Moreover, since the matrix M has only positive entries and since the fitness function is not identically equal to 1, the Perron–Frobenius eigenvalue is strictly larger than 1, i.e., $\lambda>1$ . We rewrite the equation (4.1) as

(4.3) \begin{equation}\forall\, u\in E{\qquad}\lambda x(u)-\,\sum_{v\in E}x(v)M(v,u)\,\,=\,\,\sum_{v\in E}x(v)\big(f(v)-1\big)M(v,u).\end{equation}

We denote by $x^t$ the row vector which is the transpose of the column vector x, and we introduce the square matrix whose diagonal coefficients are the fitness values $(f(u),u\in E)$ :

(4.4) \begin{equation} F\,=\,\begin{pmatrix}\ddots\;\;\;\; & 0\;\;\;\; & 0\cr0\;\;\;\; & f(u)\;\;\;\; & 0\cr0\;\;\;\; & 0\;\;\;\; & \ddots\cr\end{pmatrix}.\end{equation}

After dividing by $\lambda$ (recall that $\lambda> 1$ ), we can rewrite Equation (4.3) in matrix form as

(4.5) \begin{equation} x^t\Big(I-\frac{1}{\lambda}M\Big)\,=\,x^t\frac{1}{\lambda}\big(F-I\big)M.\end{equation}

Now the matrix M is a stochastic matrix, so its spectral radius is equal to 1. Since $\lambda>1$ , the matrix $I-\frac{1}{\lambda}M$ is invertible, and its inverse is given by the geometric series

\begin{align*} \Big(\text{I}-\frac{M}{\lambda}\Big)^{-1}\, =\,\, \sum_{n\geq 0} \frac{M^n}{\lambda^n}\,. \end{align*}

Plugging this identity into Equation (4.5), we obtain

(4.6) \begin{equation} x^t\,=\,x^t\big(F-I\big) \sum_{n\geq 1} \frac{M^n}{\lambda^n}.\end{equation}

Let us introduce the set $W^*$ of the genotypes where the fitness function is strictly larger than 1:

\begin{align*}W^*\,=\,\big\{\,u\in E\;:\;f(u)>1\,\big\}.\end{align*}

We call the set $W^*$ the set of wild types. We define the matrix $N_\lambda$ , indexed by the elements of $W^*$ , as follows:

(4.7) \begin{equation} \forall v,w\in W^*{\qquad} N_\lambda(v,w)\,=\,\Bigg(\big(F-I\big) \sum_{n\geq 1} \frac{M^n}{\lambda^n}\Bigg)(v,w).\end{equation}

With this notation, the linear system (4.6) of size $\kappa^\ell$ can be split into the system of size $|W^*|$ given by

(4.8) \begin{equation} \forall v\in W^*{\qquad} x(v)\,=\,\sum_{w\in W^*}x(w)\,N_\lambda(w,v)\,\end{equation}

and the $\kappa^\ell-|W^*|$ remaining equations

(4.9) \begin{equation} \forall v\in E\setminus W^*{\qquad} x(v)\,=\,\sum_{w\in W^*}x(w)(f(w)-1)\, \sum_{n\geq 1} \frac{M^n}{\lambda^n}(w,v).\end{equation}

The system (4.8) expresses that the vector $(x(v),\,v\in W^*)$ is a left Perron–Frobenius eigenvector of the matrix $N_\lambda$ , and that the Perron–Frobenius eigenvalue of $N_\lambda$ is equal to 1. The equations (4.9) show that the remaining coordinates $(x(v),\,v\in E\setminus W^*)$ are completely determined by $(x(v),\,v\in W^*)$ . These considerations suggest that the quasispecies equation can in principle be solved through the following procedure.

General strategy. We start by isolating the set $W^*$ of the wild types, and we form the matrix $(N_\lambda(w,v))_{v,w\in W^*}$ , defined in the formula (4.7), where $\lambda> 1$ is considered as a parameter. Now, the Perron–Frobenius eigenvalue of $N_\lambda$ is a decreasing function of $\lambda$ which tends to 0 as $\lambda$ goes to $+\infty$ and to $+\infty$ as $\lambda$ goes to 1. We choose for $\lambda$ the unique value such that the Perron–Frobenius eigenvalue of $N_\lambda$ is equal to 1. Once this value is fixed, we solve the system (4.8) on $W^*$ . The remaining coordinates of x on $E\setminus W^*$ are determined by the equations (4.9).

A general lower bound on $\lambda$ . We close this section with a general lower bound on $\lambda$ . The quasispecies equations imply that

\begin{align*}\forall u\in E{\qquad}\lambda x(u)\,\geq\,x(u)\,f(u)\,M(u,u)\,=\,x(u)f(u)\,(1-q)^\ell.\end{align*}

Since $x(u)>0$ for any $u\in E$ , this readily implies that

(4.10) \begin{equation} \lambda \,\geq\,\big(\max_Ef\big)\,(1-q)^\ell.\end{equation}

4.2. Computation of $ \boldsymbol\sum_{\boldsymbol{n}\boldsymbol\geq \textbf{1}} {\boldsymbol{M}^{\boldsymbol{n}}}/{\boldsymbol\lambda^{\boldsymbol{n}}}$

From now on, we work with the genotype space $E=\mathcal{A}^\ell$ where $\mathcal{A}$ is a finite alphabet of cardinality $\kappa$ , and we consider the mutation scheme (1.3). We start by computing $M^n$ . We use the notation and the technique of Section 2.2, and we generalize the formula (2.9) as follows:

(4.11) \begin{align}M^n(u,v)&\,=\,\mathbb{P}\big(X_{n}=v\,\big|\,X_0=u\big)\,=\,\mathbb{P}\big(X_{n}(i)=v(i), \,1\leq i\leq\ell\,\big|\,X_0=u\big)\cr&\,=\,{\prod_{1\leq i\leq\ell}}\mathbb{P}\big(X_{n}(i)=v(i)\,\big|\,X_0(i)=u(i)\big).\end{align}

Now, we have

(4.12) \begin{equation}\mathbb{P}\big(X_{n}(i)=v(i)\,\big|\,X_0(i)=u(i)\big)\,=\,\begin{cases} \frac{1}{\kappa}+ \frac{\kappa-1}{\kappa}\rho^n\phantom{\frac{1}{\Big|}}& \text{if }u(i)=v(i)\,,\cr \frac{1}{\kappa}\big(1-\rho^n \big)& \text{if }u(i)\neq v(i)\,,\cr\end{cases}\end{equation}

where

\begin{align*}\rho\,= \,1-\frac{\kappa q}{\kappa-1}.\end{align*}

Coming back to the formula (4.11), we conclude that

(4.13) \begin{equation}M^n(u,v)\,=\, \Big(\frac{1}{\kappa} + \frac{\kappa-1}{\kappa}\rho^n\Big)^{\ell-d(u,v)} \Big(\frac{1}{\kappa}\big(1-\rho^n\big) \Big)^{d(u,v)}.\end{equation}

We wish now to compute

\begin{align*} \sum_{n\geq 1} \frac{M^n}{\lambda^n}\,. \end{align*}

Our first strategy consists of expanding the formula (4.13). We end up with several geometric series, from which we obtain a closed finite formula. So, starting from (4.13) and setting $d=d(u,v)$ to simplify the notation, we have

(4.14) \begin{equation}M^n(u,v)\,=\,\sum_{i=0}^{\ell-d}\sum_{j=0}^{d}\binom{\ell-d}{i} \Big(\frac{1}{\kappa}\Big)^{\ell-d-i} \Big(\frac{\kappa-1}{\kappa}\Big)^i\rho^{ni}\binom{d}{j} \frac{1}{\kappa^d}(-1)^j\rho^{nj}.\end{equation}

We can now sum the geometric series to get

(4.15) \begin{equation}\sum_{n\geq 1} \frac{M^n}{\lambda^n}(u,v)\,=\,\sum_{i=0}^{\ell-d}\sum_{j=0}^{d}\binom{\ell-d}{i} \Big(\frac{1}{\kappa}\Big)^{\ell-d-i} \Big(\frac{\kappa-1}{\kappa}\Big)^i\binom{d}{j} \frac{1}{\kappa^d} \frac{(-1)^j}{ \frac{\lambda}{\rho^{i+j}}-1}.\end{equation}

4.3. Probabilistic interpretation

We provide furthermore a probabilistic interpretation of the apparently complex formula (4.15). Let $S_{\ell-d}$ be a binomial random variable with parameters $\ell-d,1-1/\kappa$ . We have

(4.16) \begin{equation}\sum_{n\geq 1} \frac{M^n}{\lambda^n}(u,v)\,=\,\mathbb{E}\Bigg(\frac{1}{\kappa^d}\sum_{j=0}^{d}\binom{d}{j} \frac{(-1)^j}{\frac{\lambda}{\rho^{S_{\ell-d}+j}}-1}\Bigg).\end{equation}

The above formula is quite nice. Its drawback is that it contains negative terms, so it is not obvious that the global result is non-negative. We present next a third formula, which avoids this problem.

4.4. Third formula

The trick is to expand the factor $(1-\rho^n)^d$ in the formula (4.13) in a more complicated way, as follows:

(4.17) \begin{equation}M^n(u,v)\,=\,\sum_{i=0}^{\ell-d}\binom{\ell-d}{i} \Big(\frac{1}{\kappa}\Big)^{\ell-d-i} \Big(\frac{\kappa-1}{\kappa}\Big)^i\rho^{ni} \frac{1}{\kappa^d}(1-\rho)^d\Bigg(\sum_{j=0}^{n-1}\rho^j\Bigg)^{d}.\end{equation}

We wish to compute the series $\sum_{n\geq 1} {M^n}/{\lambda^n}$ . From now on, we deal with the case $d\geq 1$ (for $d=0$ , we use the expression obtained in (4.15), where all the terms are positive). Focusing on the two terms which depend on n in the formula (4.17), we end up with the series

(4.18) \begin{align}\sum_{n\geq 1} \frac{ \rho^{ni} }{\lambda^n}\Bigg(\sum_{j=0}^{n-1}\rho^j\Bigg)^{d}&\,=\,\sum_{n\geq 1} \frac{ \rho^{ni} }{\lambda^n}\sum_{j_1=0}^{n-1} \cdots\sum_{j_d=0}^{n-1}\rho^{j_1+\cdots+j_d}\cr&\,=\,\sum_{j_1,\dots,j_d\geq 0}\,\,\sum_{n> \max\!(j_1,\dots,j_d)}\frac{ \rho^{ni} }{\lambda^n}\rho^{j_1+\cdots+j_d}\cr&\,=\,\sum_{j_1,\dots,j_d\geq 0}\,\,\frac{{\Big(\frac{ \rho^{i} }{\lambda}\Big)^{\max\!(j_1,\dots,j_d)+1}\kern-11pt}}{ 1-\frac{ \rho^{i} }{\lambda}}\,\,\rho^{j_1+\cdots+j_d}.\end{align}

Notice that this formula is not quite legitimate for $d=0$ , unless we adopt the convention that for $d=0$ there is only one term in the sum, corresponding to $j_1=\cdots=j_d=0$ . Therefore it seems safer to perform the above computation only when $d\geq 1$ . Let us set $\gamma= { \rho^{i} }/{\lambda}$ and focus on the sum

(4.19) \begin{equation}\sum_{j_1,\dots,j_d\geq 0}\,\,{\gamma^{\max\!(j_1,\dots,j_d)}}\,\,\rho^{j_1+\cdots+j_d}.\end{equation}

We shall decompose the sum according to the number r of distinct indices among $j_1,\cdots,j_d$ , their different values $i_1,\dots,i_r$ , and the number of times $n_1,\dots,n_r$ each value appears. We get

(4.20) \begin{equation}\sum_{r=1}^d\sum_{0\leq i_1<\cdots< i_r}\,\,\sum_{\genfrac{}{}{0pt}{1}{n_1,\dots,n_r\geq 1}{n_1+\cdots+n_r=d}}\sum_{j_1,\dots,j_d\in\mathcal J}{\gamma^{ \max\!(j_1,\dots,j_d)} }\,\,\rho^{j_1+\cdots+j_d},\end{equation}

where the innermost sum extends over the set $\mathcal J$ defined by

\begin{equation*}{\mathcal J}\,=\,\Big\{\,j_1,\dots,j_d\geq 0\;:\;{\text{ for }}1\leq k\leq d,\,i_k\text{ appears }{n_k\text{ times in }} {j_1,\dots,j_d}\,\Big\} .\end{equation*}

The number of terms appearing in the last summation corresponds to the number of placements of d balls in r cells with occupancy numbers given by $n_1,\dots,n_r$ (see for instance the classical book of Feller [13, Section II.5]), that is, $d!/(n_1!\cdots n_r!)$ . For each of these placements, we have

\begin{align*}{ \max\!(j_1,\dots,j_d)} \,=\,i_r\,,{\qquad}{j_1+\cdots+j_d}\,=\,n_1i_1+\cdots +n_ri_r\,;\end{align*}

therefore the formula (4.20) becomes

(4.21) \begin{equation}\sum_{r=1}^d\sum_{0\leq i_1<\cdots< i_r}\,\,\sum_{\genfrac{}{}{0pt}{1}{n_1,\dots,n_r\geq 1}{n_1+\cdots+n_r=d}}\frac{d!}{n_1!\cdots n_r!}\gamma^{ i_r}\,\,\rho^{ n_1i_1+\cdots +n_ri_r }.\end{equation}

We can perform the summation over $i_r$ , which yields

(4.22) \begin{multline}\sum_{0\leq i_1<\cdots< i_r}\,\,\gamma^{ i_r}\,\,\rho^{ n_1i_1+\cdots +n_ri_r }\,=\,\cr\frac{\gamma\rho^{n_r}} {1-\gamma\rho^{n_r}}\sum_{0\leq i_1<\cdots< i_{r-1}}\kern-3pt\big(\gamma\rho^{n_r}\big)^{ i_{r-1}}\,\,\rho^{n_1i_1+\cdots +n_{r-1}i_{r-1}}\,.\end{multline}

The number of indices $i_1,\dots,i_{r-1}$ has decreased by 1, and $\gamma$ has been replaced by $\gamma\rho^{n_r}$ . Proceeding in the same way, we perform successively the summations over $i_{r-1},\dots,i_1$ until we obtain

\begin{multline*}\sum_{0\leq i_1<\cdots< i_r}\,\,\gamma^{ i_r}\,\,\rho^{ n_1i_1+\cdots +n_ri_r }\,=\,\cr\frac{\gamma\rho^{n_r}} {1-\gamma\rho^{n_r}}\frac{\gamma\rho^{n_r+n_{r-1}}} {1-\gamma\rho^{n_r+n_{r-1}}}\cdots\frac{\gamma\rho^{n_r+\cdots+n_{2}}} {1-\gamma\rho^{n_r+\cdots+n_{2}}}\sum_{0\leq i_1}\,\,\big(\gamma\rho^{ n_r+\cdots +n_2 }\big)^{i_1}\rho^{ n_1i_1}\,.\end{multline*}

Recalling that $n_1+\cdots+n_r=d$ , we can compute the last sum

\begin{align*}\sum_{0\leq i_1}\,\,\big(\gamma\rho^{ n_r+\cdots +n_2 }\big)^{i_1}\rho^{ n_1i_1}\,=\,\frac{1} {1-\gamma\rho^{d}}\end{align*}

and conclude that

(4.23) \begin{multline}\sum_{0\leq i_1<\cdots< i_r}\,\,\gamma^{ i_r}\,\,\rho^{ n_1i_1+\cdots +n_ri_r }\,=\,\cr\frac{\gamma\rho^{n_r}} {1-\gamma\rho^{n_r}}\frac{\gamma\rho^{n_r+n_{r-1}}} {1-\gamma\rho^{n_r+n_{r-1}}}\cdots\frac{\gamma\rho^{n_r+\cdots+n_{2}}} {1-\gamma\rho^{n_r+\cdots+n_{2}}}\frac{1} {1-\gamma\rho^{d}}\,.\end{multline}

Putting together the formulas (4.21) and (4.23), we obtain a finite formula for the sum (4.19):

(4.24) \begin{multline}\sum_{j_1,\dots,j_d\geq 0}\,\,{\gamma^{\max\!(j_1,\dots,j_d)}}\,\,\rho^{j_1+\cdots+j_d}\,=\,\cr\sum_{r=1}^d\sum_{\genfrac{}{}{0pt}{1}{n_1,\dots,n_r\geq 1}{n_1+\cdots+n_r=d}}\frac{d!}{n_1!\cdots n_r!}\frac{1}{\gamma\rho^{d}}\frac{\rho^{rn_r+(r-1)n_{r-1}+\cdots+n_{1}}}{\prod_{k=1}^r \Big({\frac{1}{\gamma}-\rho^{n_r+\cdots+n_{r-k+1}}}\Big)}\,.\end{multline}

Plugging this formula into (4.18) and (4.17) yields that, for $d=d(u,v)\geq 1$ ,

(4.25) \begin{multline}\sum_{n\geq 1} \frac{M^n}{\lambda^n}(u,v)\,=\,\sum_{i=0}^{\ell-d}\binom{\ell-d}{i} \Big(\frac{1}{\kappa}\Big)^{\ell-i} \Big(\frac{\kappa-1}{\kappa}\Big)^i \Big(\frac{1-\rho}{\rho}\Big)^d\frac{ 1 }{1-\frac{\rho^i}{\lambda}}\cr\times\sum_{r=1}^d\sum_{\genfrac{}{}{0pt}{1}{n_1,\dots,n_r\geq 1}{n_1+\cdots+n_r=d}}\frac{d!}{n_1!\cdots n_r!}\frac{\rho^{rn_r+(r-1)n_{r-1}+\cdots+n_{1}}}{\prod_{k=1}^r \Big({\frac{\lambda}{\rho^i}-\rho^{n_r+\cdots+n_{r-k+1}}}\Big)}\,.\end{multline}

We can further rewrite this formula using the expectation with respect to a binomial random variable $S_{\ell-d}$ with parameters $\ell-d$ and $1-1/\kappa$ . Indeed, we have

(4.26) \begin{multline}\sum_{n\geq 1} \frac{M^n}{\lambda^n}(u,v)\,=\, \frac{1}{\kappa^d} \Big(\frac{1-\rho}{\rho}\Big)^d\times\cr\mathbb{E}\Bigg(\frac{ 1 }{1-\frac{\rho^{S_{\ell-d}}}{\lambda}}\sum_{r=1}^d\sum_{\genfrac{}{}{0pt}{1}{n_1,\dots,n_r\geq 1}{n_1+\cdots+n_r=d}}\frac{d!}{n_1!\cdots n_r!}\frac{\rho^{rn_r+(r-1)n_{r-1}+\cdots+n_{1}}}{\prod_{k=1}^r\Big(\frac{\lambda}{\rho^{S_{\ell-d}}}-\rho^{n_r+\cdots+n_{r-k+1}}\Big)}\Bigg)\,.\end{multline}

Each of the previous formulas is of interest on its own and may be useful in certain contexts. For instance, with the formula (4.15), we see directly that the coefficients of the matrix $N_\lambda$ are rational functions of the parameter $\lambda$ , and that there exists a unique choice for $\lambda$ such that the Perron–Frobenius eigenvalue of $N_\lambda$ becomes equal to 1. The analysis of the error threshold conducted in Section 3.2 rests entirely on the probabilistic representation presented in the formula (4.16) in the specific case where $u=v=w^*$ and $d=0$ . Finally, the third formula (4.26) will be useful for analyzing the asymptotic behavior of the non-diagonal entries of the matrix $N_\lambda$ for the landscape with finitely many peaks, in Section 5.

4.5. Proof of Theorem 2

Let $G\;:\;E\to\mathbb{R}$ be an additive functional, that is, a function given by

\begin{align*}\forall u\in E{\qquad} G(u)\,=\,\sum_{i=1}^\ell g(u(i))\,,\end{align*}

where g is a function defined on $\mathcal{A}$ with values in $\mathbb{R}$ . We consider an arbitrary fitness function satisfying Hypothesis 1, and we denote by $(x(u), u\in E)$ the solution of the quasispecies equation associated with f. Our goal here is to relate the mean value $\overline{G}$ of G in the quasispecies $(x(u), u\in E)$ , defined by

\begin{align*}\overline{G}\,=\,\sum_{u\in E}G(u)x(u)\,,\end{align*}

to the mean fitness $\lambda$ or $\overline{f}$ , defined by

\begin{align*}\overline{f}\,=\,\sum_{u\in E}f(u)x(u).\end{align*}

Our starting point is the formula (4.6), which yields

(4.27) \begin{align}\overline{G} &\,=\,\sum_{u\in E}G(u)\sum_{v\in E} x(v)\big(f(v)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M^n(v,u) \cr&\,=\,\sum_{v\in E} x(v)\big(f(v)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}\sum_{u\in E} M^n(v,u)\, G(u) .\end{align}

We compute next the innermost sum. By the definition of G,

(4.28) \begin{equation} \sum_{u\in E} M^n(v,u)\, G(u)\,=\,\sum_{i=1}^\ell\sum_{u\in E} M^n(v,u)\, g(u(i)).\end{equation}

Moreover, the formulas (4.11) and (4.12) yield that, for $1\leq i\leq\ell$ ,

(4.29) \begin{align}\sum_{u\in E} M^n(v,u)\, g(u(i))&\,=\,\Big( \frac{1}{\kappa} + \frac{\kappa-1}{\kappa}\rho^n \Big)g(v(i))+\Big( \frac{1}{\kappa}\big(1-\rho^n \big) \Big)\sum_{a\in\mathcal{A}\setminus\{v(i)\}}g(a)\cr&\,=\, \rho^n g(v(i)) + \big(1-\rho^n \big) \overline{g}\,,\kern20pt\end{align}

where $\overline{g}$ is the mean value of g over the alphabet $\mathcal{A}$ , i.e.,

(4.30) \begin{equation} \overline{g}\,=\,\frac{1}{\kappa}\sum_{a\in\mathcal{A}}g(a).\end{equation}

Plugging the formula (4.29) into the formula (4.28), we obtain

(4.31) \begin{equation} \sum_{u\in E} M^n(v,u)\, G(u)\,=\, \rho^n G(v) + \ell\big(1-\rho^n \big) \overline{g}.\end{equation}

Inserting this last formula into (4.27), we obtain

(4.32) \begin{align}\overline{G} &\,=\,\sum_{v\in E} x(v)\big(f(v)-1\big)\sum_{n\geq 1} \frac{1}{\lambda^n}\Big( \rho^n G(v) + \ell\big(1-\rho^n \big) \overline{g}\Big) \cr&\,=\,\sum_{v\in E} x(v)\big(f(v)-1\big)\Big( \frac{\rho}{\lambda-\rho} G(v) + \frac{\ell\overline{g}}{\lambda-1}- \frac{\ell\overline{g}\rho}{\lambda-\rho} \Big) \cr&\,=\, \frac{\rho}{\lambda-\rho} \big( \overline{fG}-\overline{G}\big) + \big( \overline{f}-1\big) {\ell\overline{g}} \frac{\lambda(1-\rho)}{(\lambda-1)(\lambda-\rho)} \,,\end{align}

where we have introduced the notation

(4.33) \begin{equation} \overline{fG}\,=\,\sum_{u\in E}f(u)\,G(u)\,x(u).\end{equation}

Recalling that $\overline{f}=\lambda$ , we have proved the curious formula stated in Theorem 2.

5. Finitely many peaks

As we did for the sharp-peak landscape, we shall give two proofs of Theorem 3: a soft one, whose starting point is the initial quasispecies equation, and a technical one, based on the strategy explained at the end of Section 4. Again, the soft proof might look easier, but the technical proof is more informative and could potentially yield more accurate estimates. We define first the set $W^*$ of the wild types as

\begin{align*}W^*\,=\,\big\{\,w_i\;:\;1\leq i\leq k\,\big\}.\end{align*}

Throughout the proofs, we write simply $\lambda(\ell,q)$ or even $\lambda$ instead of $\lambda_{\text{FP}}(\ell,q)$ .

5.1. Soft proof

We start as in the case of the sharp-peak landscape. In the long-chain regime, we have that $\sigma\smash{ \big(1-q\big)^\ell}\to\sigma e^{-a}$ . Using the lower bound (4.10) and the fact that $\lambda\geq 1$ , we see that

(5.1) \begin{equation}\liminf\,\,\lambda\,\geq\,\max\big(1, \,\sigma e^{-a}\big).\end{equation}

Suppose that $\lambda$ does not converge towards $\max\big(1, \,\sigma e^{-a}\big)$ . In this case, there exist $\varepsilon>0$ and two sequences of parameters $(\ell_n)_{n\geq 0}$ , $(q_n)_{n\geq 0}$ such that

\begin{equation*}\ell_n\to+\infty\,,{\qquad}q_n\to 0\,,{\qquad}\ell_n q_n\to a,\end{equation*}

(5.2) \begin{equation}\forall n\geq 0{\qquad}\lambda(\ell_n,q_n)\,\geq\,\max\big(1, \,\sigma e^{-a}\big)+\varepsilon.\end{equation}

For $1\leq i\leq k$ , we bound from above the quasispecies equation (1.1) associated with $w^*_i$ , as follows:

(5.3) \begin{equation}\lambda x(w^*_i)\,\leq\,\sum_{v\in {W^*}}\sigma x(v)M(v,w^*_i)\,+\,\sum_{v\in E\setminus {W^*}}x(v)M(v,w^*_i).\end{equation}

We sum Equation (5.3) over $i\in\{\,1,\dots,k\,\}$ . Setting

\begin{align*}&\qquad x(W^*)\,=\,\sum_{1\leq i\leq k}x(w^*_i)\,,\\[5pt] &\forall v\in E{\qquad} M(v,W^*)\,=\,\sum_{1\leq i\leq k}M(v,w^*_i)\,,\end{align*}

we obtain

(5.4) \begin{equation}\lambda x(W^*)\,\leq\,\sum_{v\in {W^*}}\sigma x(v)M(v,W^*)\,+\,\sum_{v\in E\setminus {W^*}}x(v)M(v,W^*)\,.\end{equation}

Next, using the fact that $M(v,w^*_i)\leq q$ whenever $v\neq w^*_i$ , we obtain that

\begin{align*}\forall v\in W^*{\qquad} M(v,W^*)\,\leq\,M(v,v)+kq\,=\,(1-q)^\ell +kq\,,\end{align*}

\begin{align*}\forall v\in E\setminus W^*{\qquad} M(v,W^*)\,\leq\,kq.\end{align*}

Inserting these inequalities into (5.4), we obtain that

(5.5) \begin{align}\lambda x(W^*)&\,\leq\,\sigma x(W^*)\big((1-q)^\ell+kq\big)+\,\big(1-x(W^*)\big)kq\cr&\,\leq\,\sigma x(W^*)(1-q)^\ell+\sigma kq.\end{align}

In addition, we have that

(5.6) \begin{equation}\lambda\,\leq\, \sigma x(W^*)+1-x(W^*)\,\leq\, \sigma x(W^*)+1.\end{equation}

Equation (5.6) and the condition (5.2) yield that

(5.7) \begin{equation}\frac{1}{x(W^*)}\,\leq\,\frac{ \sigma}{\lambda-1}\,\leq\,\frac{ \sigma}{\varepsilon}.\end{equation}

Using (5.5) and (5.7) together, we obtain that

\begin{equation*}\forall n\geq 0{\qquad}\lambda(\ell_n,q_n)\,\leq\,\sigma (1-q_n)^{\ell_n}+\frac{\sigma}{\varepsilon}\sigma kq_n.\end{equation*}

Sending n to $\infty$ , we see that this inequality is not coherent with the inequality (5.2). Therefore it must be the case that $\lambda$ converges towards $\max\big(1, \,\sigma e^{-a}\big)$ .

5.2. Technical proof

For this proof, we shall implement the strategy explained at the end of Section 4. So we consider the matrix $\big(N_\lambda(i,j)\big)_{1\leq i,j\leq k}$ defined by

(5.8) \begin{equation} \forall i,j\in \{\,1,\dots,k\,\}{\qquad} N_\lambda(i,j)\,=\,\Bigg(\big(F-I\big) \sum_{n\geq 1} \frac{M^n}{\lambda^n}\Bigg)(w^*_i,w^*_j),\end{equation}

where the matrix F was introduced in the formula (4.4). The advantage is that the size of the matrix $N_\lambda$ is $k\times k$ and this size does not vary with $\ell$ and q. The diagonal elements of the matrix are given by

\begin{align*}N_\lambda(i,i)\,=\,(\sigma_i-1) \sum_{n\geq 1} \frac{1}{\lambda^n}M^n(w^*_i,w^*_i)\,,\quad 1\leq i\leq k\,. \end{align*}

In fact, the above series does not depend on i. Using the formula (4.16) with $d=0$ , we have

(5.9) \begin{equation} \forall i\in\{\,1,\dots,k\,\}{\qquad} N_\lambda(i,i)\,=\,(\sigma_i-1)\,\mathbb{E}\Bigg( \frac{1}{ \frac{\lambda}{\rho^{S_{\ell}}}-1}\Bigg),\end{equation}

where $S_{\ell}$ is a binomial random variable with parameters $\ell,1-1/\kappa$ . The parameter $\lambda$ is adjusted so that the Perron–Frobenius eigenvalue of the matrix $N_\lambda$ is equal to 1. This Perron–Frobenius eigenvalue is always larger than or equal to any diagonal element; thus we must have $N_\lambda(i,i)\,\leq \,1$ . Proceeding exactly as in the technical proof of Theorem 1, we obtain with the help of Jensen’s inequality that

(5.10) \begin{equation} \forall i\in\{\,1,\dots,k\,\}{\qquad}\lambda\,\geq\,\sigma_i\rho^{{(\kappa-1)\ell}/{\kappa}}.\end{equation}

Taking the maximum over $i\in\{\,1,\dots,k\,\}$ , we get

(5.11) \begin{equation}\lambda\,\geq\,\sigma\rho^{{(\kappa-1)\ell}/{\kappa}}.\end{equation}

Next we seek an upper bound on the value of $\lambda$ . The Perron–Frobenius eigenvalue of the matrix $N_\lambda$ is equal to 1, and all its entries are non-negative; thus

\begin{align*}1\,\leq\,\max_{1\leq i\leq k}\sum_{1\leq j\leq k}N_{\lambda}(i,j).\end{align*}

In particular, there exists an index $i\in\{\,1,\dots,k\,\}$ such that

(5.12) \begin{equation}1\,\leq\,\sum_{1\leq j\leq k}N_{\lambda}(i,j).\end{equation}

Let us study the non-diagonal elements of the matrix $N_\lambda$ . For $i\neq j$ , we have

\begin{align*}N_\lambda(i,j)\,=\,(\sigma_i-1) \sum_{n\geq 1} \frac{1}{\lambda^n}M^n(w^*_i,w^*_j)\,. \end{align*}

In fact, the series depends only on the Hamming distance between $w^*_i$ and $w^*_j$ , as we can see from the formula (4.13). The same formula shows also that $M^n(w^*_i,w^*_j)$ is a non-increasing function of $d(w^*_i,w^*_j)$ . Therefore, using the formula (4.26) with $d=1$ , we have the inequality

(5.13) \begin{equation} N_\lambda(i,j)\,\leq\,(\sigma_i-1)\frac{\lambda(1-\rho)}{\kappa}\,\mathbb{E}\Bigg(\frac{ 1}{{\rho^{S_{\ell-1}}}\Big( \frac{\lambda}{\rho^{S_{\ell-1}}}-1\Big)\Big({ \frac{\lambda}{\rho^{S_{\ell-1}}}-\rho}\Big)}\Bigg),\end{equation}

where $S_{\ell-1}$ is a binomial random variable with parameters $\ell-1,1-1/\kappa$ . Moreover we know that $1\leq \lambda\leq\sigma$ , so

(5.14) \begin{equation} N_\lambda(i,j)\,\leq\,(\sigma_i-1)\sigma\frac{(1-\rho)}{\kappa}\,\mathbb{E}\Bigg(\frac{ 1}{\Big( \frac{\lambda}{\rho^{S_{\ell-1}}}-1\Big)\Big({ \lambda-{\rho^{S_{\ell-1}+1}}}\Big)}\Bigg).\end{equation}

Plugging (5.9) and (5.14) into the inequality (5.12), and recalling that $\sigma\geq \sigma_i$ , we get

(5.15) \begin{equation} \frac{1}{\sigma-1}\,\leq\,\mathbb{E}\Bigg(\frac{1}{ \frac{\lambda}{\rho^{S_{\ell}}}-1}\Bigg)\, +\,k\sigma \frac{(1-\rho)}{\kappa}\,\mathbb{E}\Bigg(\frac{ 1}{\Big( \frac{\lambda}{\rho^{S_{\ell-1}}}-1\Big)\Big({ \lambda-{\rho^{S_{\ell-1}+1}}}\Big)}\Bigg).\end{equation}

We proceed in a way similar to what we did for the sharp-peak landscape. Namely, we pick a positive number $\varepsilon$ and split the expectations according to the values of $S_\ell$ and $S_{\ell-1}$ , as follows:

(5.16) \begin{multline}\frac{1}{\sigma-1}\,\leq\,\frac{1}{ {\lambda}-1}\mathbb{P}\Big(\Big|S_\ell-\frac{\kappa-1}{\kappa}\ell\Big|>\varepsilon\ell\Big)+\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)\ell/\kappa-\varepsilon\ell} }-1}\cr+ \frac{k\sigma(1-\rho)}{\kappa ( {\lambda}-1) ( {\lambda}-\rho)}\mathbb{P}\Big(\Big|S_{\ell-1}-\frac{\kappa-1}{\kappa}(\ell-1)\Big|>\varepsilon(\ell-1)\Big)\cr+ \frac{k\sigma(1-\rho)}{\kappa}\frac{1}{ \Big( \frac{\lambda}{ \rho^{(\kappa-1)(\ell-1)/\kappa-\varepsilon\ell} }-1\Big) \Big( {\lambda}-{ \rho^{(\kappa-1)(\ell-1)/\kappa-\varepsilon\ell} }\Big)}\,.\end{multline}

Using the estimate given by Chebyshev’s inequality (3.3), we obtain

(5.17) \begin{multline}\frac{1}{\sigma-1}\,\leq\,\frac{1}{ {\lambda}-1}\frac{\kappa-1}{\kappa^2\varepsilon^2\ell}+\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)\ell/\kappa-\varepsilon\ell} }-1}+ \frac{k\sigma(1-\rho)}{\kappa ( {\lambda}-1) ( {\lambda}-\rho)}\frac{\kappa-1}{\kappa^2\varepsilon^2(\ell-1)}\cr+ \frac{k\sigma(1-\rho)}{\kappa}\frac{1}{ \Big( \frac{\lambda}{ \rho^{(\kappa-1)(\ell-1)/\kappa-\varepsilon\ell} }-1\Big) \Big( {\lambda}-{ \rho^{(\kappa-1)(\ell-1)/\kappa-\varepsilon\ell} }\Big)}\,.\end{multline}

Finally, regrouping the first and the third term together, as well as the second and the fourth term, and using that $1-\rho\leq \lambda-\rho$ and $\rho\leq 1$ , we get

(5.18) \begin{multline}\frac{1}{\sigma-1}\,\leq\, \frac{\kappa+k\sigma}{\kappa ({\lambda}-1) }\frac{\kappa-1}{\kappa^2\varepsilon^2(\ell-1)}\cr+\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)(\ell-1)/\kappa-\varepsilon\ell} }-1}\Bigg(1+\frac{ \frac{k\sigma(1-\rho)}{\kappa} }{ \Big( {\lambda}-{ \rho^{(\kappa-1)(\ell-1)/\kappa-\varepsilon\ell} }\Big)} \Bigg)\,.\end{multline}

By making an appropriate choice for $\varepsilon$ , we shall now deduce an upper bound on $\lambda$ from this inequality.

Lemma 2. There exists a positive constant c such that, in the long-chain regime, the mean fitness $\lambda$ satisfies

(5.19) \begin{equation}\lambda\,\leq\,\max\bigg(1+\frac{c}{\ell^{1/4}},\sigma \Big(1- \frac{\kappa q}{\kappa-1}\Big)^{{(\kappa-1)(\ell-1)}/{\kappa}-\ell^{2/3}}\Big(1+\frac{1}{\ell^{1/13}}\Big)\bigg).\end{equation}

Before proving Lemma 2, we explain the end of the technical proof. The lower bound (5.11) on $\lambda$ and the upper bound (5.19) on $\lambda$ yield the desired conclusion, which terminates the technical proof. We have also obtained estimates on the speed of convergence of $\lambda$ , which could be improved with some additional work.

Proof of lemma 2. Suppose that $\lambda$ is larger than the upper bound in (5.19). We check that this is not compatible with the inequality (5.18) with the choice $\varepsilon={\ell}^{-1/3}$ . Indeed, we have on the one hand

\begin{align*}\lambda>1+\frac{c}{\ell^{1/4}}{\qquad}\Longrightarrow{\qquad} \frac{\kappa+k\sigma}{\kappa ({\lambda}-1) }\frac{\kappa-1}{\kappa^2 \ell^{-2/3} (\ell-1)}\,=\,O\Big( \frac{1}{\ell^{1/12}}\Big).\end{align*}

On the other hand,

\begin{multline*}\lambda>\sigma \Big(1- \frac{\kappa q}{\kappa-1}\Big)^{{(\kappa-1)(\ell-1)}/{\kappa}-\ell^{2/3}}\Big(1+\frac{1}{\ell^{1/13}}\Big){\qquad}\Longrightarrow{\qquad}\cr\exists c'>0{\qquad}\frac{1}{ \frac{\lambda}{\rho^{(\kappa-1)(\ell-1)/\kappa-\ell^{2/3}} }-1}\,\leq\,\frac{1}{ \sigma -1}-\frac{c'}{\ell^{1/13} }+o\Big( \frac{1}{\ell^{1/13} } \Big)\,.\end{multline*}

We have also that

\begin{equation*}\frac{ \frac{k\sigma(1-\rho)}{\kappa} }{ \Big( {\lambda}-{ \rho^{(\kappa-1)(\ell-1)/\kappa-\ell^{2/3}} }\Big)} \,=\,O\Big( \frac{1}{\ell}\Big)\end{equation*}

(notice that we use here the fact that a is finite). Combining the previous inequalities, we find that the right-hand side of the inequality (5.18) has the following expansion:

\begin{align*}O\Big( \frac{1}{\ell^{1/12}}\Big)+\frac{1}{ \sigma -1}-\frac{c'}{\ell^{1/13} }+o\Big( \frac{1}{\ell^{1/13} } \Big)+O\Big( \frac{1}{\ell}\Big)\,.\end{align*}

This quantity becomes strictly less than $1/(\sigma-1)$ for $\ell$ large enough, and this is not compatible with the inequality (5.18).

6. Plateau

This section is devoted to the proof of Theorems 4 and 5. We first show that the value $\lambda(a,\sigma)$ is well-defined. Then we perform a lumping procedure and study the asymptotics of the reduced equations. This yields the convergence of $\lambda$ towards $\lambda(a,\sigma)$ . The other claims of the theorem are proved by similar arguments.

6.1. Definition of $\lambda(a,\sigma)$

We start by proving that Equation (1.13) admits a unique solution. Let us define the function $\phi_n$ by

(6.1) \begin{equation}\phi_n(a)\,=\,e^{-an}\sum_{k\geq 0}\frac{(an)^{2k}}{2^{2k}(k!)^2}.\end{equation}

For any $a>0$ , we have

\begin{align*} \phi_n(a)\,=\, e^{-an}\sum_{k\geq 0}\Bigg( \frac{1}{k!} \Big(\frac{an}{2}\Big)^k\Bigg)^2\,\leq\, e^{-an}\Bigg(\sum_{k\geq 0} \frac{1}{k!} \Big(\frac{an}{2}\Big)^k\Bigg)^2\,=\,1\,.\end{align*}

It follows that the map

\begin{align*}\lambda\mapsto\sum_{n\geq 1}\frac{{\phi_n(a)}}{\lambda^n}\end{align*}

is continuous and decreasing on $]1,+\infty[$ . Moreover it converges to 0 when $\lambda$ goes to $\infty$ and to $+\infty$ when $\lambda$ goes to 1, because

\begin{align*}&\sum_{n\geq 1}{{\phi_n(a)}}\,\geq\,\sum_{n\geq 1}e^{-an}\sum_{k\geq 0}\frac{(an)^{2k}}{(2k+1)!}\\[5pt] &\qquad =\sum_{n\geq 1}\frac{e^{-an}}{an}\sum_{k\geq 0}\frac{(an)^{2k+1}}{(2k+1)!}\,=\,\sum_{n\geq 1}\frac{e^{-an}}{an}\text{sh}(an)\,=\,+\infty\,.\end{align*}

Thus there exists a unique value $\lambda(a,\sigma)$ satisfying Equation (1.13). We prove next the main claim of the theorem. Throughout the proof, we write simply $\lambda(\ell,q)$ or even $\lambda$ instead of $\lambda_{\text{PL}}(\ell,q)$ .

6.2. Lumping

We are dealing with a fitness function which depends only on the Hamming classes. So we lump together the genotypes according to their Hamming classes and introduce the new variables y(h), $0\leq h\leq\ell$ , given by

\begin{align*}\forall h\in\{\,0,\dots,\ell\,\}{\qquad} y(h)\,=\,\sum_{u:H(u)=h}x(u).\end{align*}

Using the formula (4.6), we have

(6.2) \begin{multline}\forall h\in\{\,0,\dots,\ell\,\}{\qquad} y(h)\,=\,\sum_{u:H(u)=h}\sum_{v} x(v)\big(f(v)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M^n(v,u)\,\cr\,=\,\sum_{v}x(v)\big(f(v)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}\sum_{u:H(u)=h}M^n(v,u)\,.\end{multline}

Now, the point is that the innermost sum depends only on the Hamming class of v. Indeed, let us fix $b\in\{\,0,\dots,\ell\,\}$ and let v be such that $H(v)=b$ . After n mutation steps starting from v, we obtain a genotype $X_n$ whose components are $\ell$ independent Bernoulli random variables. Among them, exactly b have parameter $(1+\rho^n)/2$ , and $\ell-b$ have parameter $(1-\rho^n)/2$ (this is a consequence of the computations performed in Section 2.2). Therefore the Hamming class of $X_n$ is distributed as the sum of two independent binomial random variables with respective parameters $(b,(1+\rho^n)/2)$ and $(\ell-b,(1-\rho^n)/2)$ . We conclude that

(6.3) \begin{multline}\forall\, b,c\in\{\,0,\dots,\ell\,\}\,,\quad\forall v\in\lbrace\,0,1\,\rbrace^\ell\,,{\qquad}\cr H(v)=b\quad\Longrightarrow\quad\sum_{u:H(u)=c}M^n(v,u)\,=\,M^n_H(b,c)\,,\end{multline}

where the matrix $M^n_H(b,c)$ is given by

(6.4) \begin{equation}M^n_H(b,c)\,=\,\mathbb{P}\bigg( \text{Bin}\Big(b, \frac{1+\rho^n}{2}\Big)\,+\, \text{Bin}\Big(\ell-b, \frac{1-\rho^n}{2}\Big)\,=\,c \bigg),\end{equation}

$\text{Bin}(n,p)$ being a generic binomial random variable with parameters (n, p), and the two binomial random variables appearing above being independent. From the definition (1.12), we see that the fitness $f_{\text{PL}}$ depends only on the Hamming class, i.e., there exists a function f such that $f_{\text{PL}}(v)=f(H(v))$ for any v. With a slight abuse of notation, we denote this function also by $f_{\text{PL}}$ . Rearranging the sum in the formula (6.2) according to the Hamming class of v and using (6.3), we get

(6.5) \begin{align}\forall h\in\{\,0,\dots,\ell\,\}{\qquad}y(h)&\,=\,\sum_{0\leq k\leq\ell}\sum_{v:H(v)=k}x(v)\big(f_{\text{PL}}(v)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n(k,h)\,\cr&\,=\,\sum_{0\leq k\leq\ell}y(k)\big(f_{\text{PL}}(k)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n(k,h).\end{align}

Equation (6.5) is in fact valid as long as the fitness function depends only on the Hamming classes. In the case of the plateau landscape, the equation for $h=\ell/2$ yields

(6.6) \begin{equation}\frac{1}{\sigma-1}\,=\, \sum_{n\geq 1} \frac{1}{\lambda^n}M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big).\end{equation}

This is the counterpart of Equation (2.7) for the sharp-peak landscape.

6.3. Asymptotics of $\boldsymbol{M}^{\boldsymbol{n}}_{\boldsymbol{H}}({\boldsymbol\ell}/{\textbf{2}},\boldsymbol\ell/\textbf{2})$

We shall rely on Equation (6.6) to analyze the asymptotic behavior of $\lambda$ in the long-chain regime. Let us begin with the term $M^n_H({\ell}/{2},\ell/2)$ .

Lemma 3. In the long-chain regime, we have the convergence

(6.7) \begin{equation}\forall n\geq 1{\qquad}\lim_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\,=\,\phi_n(a),\end{equation}

where the sequence of functions $\phi_n$ is defined in (6.1).

Proof. From the formula (6.4) with $b=c=\ell/2$ , we have

\begin{align*}M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)&\,=\,\mathbb{P}\bigg( \text{Bin}\Big(\frac {\ell}{2}, \frac{1+\rho^n}{2}\Big)\,+\, \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big)\,=\,\frac {\ell}{2} \bigg)\cr&\,=\,\mathbb{P}\bigg( \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big)\,=\, \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big) \bigg)\,,\end{align*}

with the understanding that the two binomial random variables appearing above are independent. Moreover, for $a>0$ we have

\begin{align*}\frac {\ell}{2} \frac{1-\rho^n}{2}\,=\,\frac {\ell}{2} \frac{1-(1-2q)^n}{2}\,\sim\,\frac{\ell}{2} nq\,\sim\,\frac{na}{2}\,;\end{align*}

hence we are in the regime where these binomial laws converge towards the Poisson distribution $\mathcal{P}(na/2)$ with parameter $na/2$ . If $a=0$ , the limit distribution is a Dirac mass at 0. Let us fix $K\geq 1$ . We have

\begin{equation*}M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\,\geq\,\sum_{0\leq k\leq K}\mathbb{P}\bigg( \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big) \,=\,k\bigg)^2,\end{equation*}

whence, passing to the limit,

\begin{equation*}\liminf_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\,\geq\,e^{-an}\sum_{0\leq k\leq K}\frac{(an)^{2k}}{2^{2k}(k!)^2}.\end{equation*}

Sending K to $\infty$ , we obtain that

(6.8) \begin{equation}\liminf_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\,\geq\,e^{-an}\sum_{k\geq 0}\frac{(an)^{2k}}{2^{2k}(k!)^2}\,=\,\phi_n(a).\end{equation}

Let us look for the reverse inequality. We again fix $K\geq 1$ and write

\begin{equation*}M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\,\leq\,\sum_{0\leq k\leq K}\mathbb{P}\bigg( \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big) \,=\,k\bigg)^2\!\!+\,\mathbb{P}\bigg( \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big) \,>\,K\bigg).\end{equation*}

We bound the last term with the help of Markov’s inequality:

\begin{equation*}\mathbb{P}\bigg( \text{Bin}\Big(\frac {\ell}{2}, \frac{1-\rho^n}{2}\Big) \,>\,K\bigg)\,\leq\, \frac{1}{K} \frac {\ell}{2} \frac{1-\rho^n}{2}\,\leq\, \frac{an}{K} ,\end{equation*}

where the last inequality holds asymptotically. From the previous two inequalities, we conclude that

\begin{equation*}\limsup_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\,\leq\,e^{-an}\sum_{0\leq k\leq K}\frac{(an)^{2k}}{2^{2k}(k!)^2}\,+\, \frac{an}{K}.\end{equation*}

Finally, we send K to $\infty$ to complete the proof of the lemma.

6.4. Convergence of $\lambda$

Let $\lambda^*$ be an accumulation point of $\lambda$ . By definition, there exist two sequences of parameters $(\ell_m)_{m\geq 0}$ , $(q_m)_{m\geq 0}$ such that

\begin{equation*}\ell_m\to+\infty\,,{\qquad}q_m\to 0\,,{\qquad}\ell_m q_m\to a\,,{\qquad}\lambda(\ell_m, q_m)\to \lambda^*.\end{equation*}

By Lemma 3 and Fatou’s lemma, we have

(6.9) \begin{align}\sum_{n\geq 1} \frac{ \phi_n(a) }{(\lambda^*)^n}&\,=\, \sum_{n\geq 1} \lim_{m\to\infty} \frac{1}{\lambda(\ell_m, q_m)^n}M^n_H\Big(\frac {\ell_m}{2},\frac {\ell_m}{2}\Big)\cr&\,\leq\, \lim_{m\to\infty} \sum_{n\geq 1} \frac{1}{\lambda(\ell_m, q_m)^n}M^n_H\Big(\frac {\ell_m}{2},\frac {\ell_m}{2}\Big)\,=\, \frac{1}{\sigma-1}\,,\end{align}

where the last equality comes from (6.6). This inequality and the very definition of $\lambda(a,\sigma)$ (see Equation (1.13)) imply that $\lambda^*\geq\lambda(a,\sigma)>1$ . In particular, there exist $\varepsilon>0$ and $M\geq 1$ such that

(6.10) \begin{equation}\forall m\geq M{\qquad}\lambda(\ell_m,q_m)\,\geq\,1+\varepsilon.\end{equation}

This last inequality ensures that the series

\begin{align*}\sum_{n\geq 1} \frac{1}{\lambda(\ell_m, q_m)^n}M^n_H\Big(\frac {\ell_m}{2},\frac {\ell_m}{2}\Big)\end{align*}

converges uniformly with respect to m, so that we can interchange the limits and the infinite sum to get

(6.11) \begin{align}\frac{1}{\sigma-1}&\,=\, \lim_{m\to\infty} \sum_{n\geq 1} \frac{1}{\lambda(\ell_m, q_m)^n}M^n_H\Big(\frac {\ell_m}{2},\frac {\ell_m}{2}\Big)\cr &\,=\, \sum_{n\geq 1} \lim_{m\to\infty} \frac{1}{\lambda(\ell_m, q_m)^n}M^n_H\Big(\frac {\ell_m}{2},\frac {\ell_m}{2}\Big)\,=\, \sum_{n\geq 1} \frac{\phi_n(a) }{(\lambda^*)^n}.\end{align}

Thus we see that $\lambda^*$ has to be the solution $\lambda(a,\sigma)$ of Equation (1.13). In conclusion, the only possible accumulation point for $\lambda$ in the long-chain regime is $\lambda(a,\sigma)$ . Therefore $\lambda$ converges towards $\lambda(a,\sigma)$ as stated, and this completes the proof of Theorem 4.

6.5. Asymptotics of $\boldsymbol{M}^{\boldsymbol{n}}_{\boldsymbol{H}}\big(\lfloor {\boldsymbol\alpha\boldsymbol\ell}\rfloor, \lfloor {\boldsymbol\alpha\boldsymbol\ell}\rfloor \big)$

The next lemma is the missing ingredient needed to complete the proof of Theorem 5.

Lemma 4. In the long-chain regime, we have the convergence

(6.12) \begin{equation}\forall n\geq 1{\qquad}\lim_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }M^n_H\big(\lfloor {\alpha\ell}\rfloor, \lfloor {\alpha\ell}\rfloor \big)\,=\,\phi_{n,\alpha}(a),\end{equation}

where the sequence of functions $\phi_{n,\alpha}$ is defined by

(6.13) \begin{equation}\phi_{n,\alpha}(a)\,=\,e^{-an}\sum_{k\geq 0}\frac{(\sqrt{4\alpha(1-\alpha)}an)^{2k}}{2^{2k}(k!)^2}.\end{equation}

Proof. From the formula (6.4) with $b=c=\lfloor\alpha\ell\rfloor$ , we have

\begin{align*}M^n_H\big(\lfloor {\alpha\ell}\rfloor, \lfloor {\alpha\ell}\rfloor \big)&\,=\,\mathbb{P}\Big( \text{Bin}\Big(\lfloor {\alpha\ell}\big\rfloor, \frac{1+\rho^n}{2}\Big)\,+\,\text{Bin}\Big({\ell}-\lfloor {\alpha\ell}\rfloor, \frac{1-\rho^n}{2}\Big)\,=\,\lfloor {\alpha\ell}\rfloor \Big)\cr&\,=\,\mathbb{P}\Big( \text{Bin}\Big({\ell}-\lfloor {\alpha\ell}\rfloor, \frac{1-\rho^n}{2}\Big)\,=\, \text{Bin}\Big(\lfloor {\alpha\ell}\rfloor, \frac{1-\rho^n}{2}\Big) \Big)\,,\end{align*}

with the understanding that the two binomial random variables appearing above are independent. Moreover, we have

\begin{align*}{\ell} \frac{1-\rho^n}{2}\,=\,{\ell} \frac{1-(1-2q)^n}{2}\,\sim\,{\ell} nq\,\sim\,{na}\,,\end{align*}

so we are in the regime where the first binomial law converges towards the Poisson distribution with parameter $(1-\alpha)na$ while the second binomial law converges towards the Poisson distribution with parameter $\alpha na$ . We proceed as in the proof of Lemma 3 to conclude that

\begin{equation*}\forall n\geq 1{\qquad}\lim_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }M^n_H\big(\lfloor {\alpha\ell}\rfloor, \lfloor {\alpha\ell}\rfloor \big)\,=\,e^{-an}\sum_{k\geq 0}\frac{((1-\alpha)an)^{k}}{k!}\frac{(\alpha an)^{k}}{k!},\end{equation*}

and we rewrite the right-hand quantity as in the formula (6.13).

6.6. Completion of the proof of Theorem 5

Let $\phi_{\alpha}(a)$ be the function defined by

(6.14) \begin{equation}\phi_{\alpha}(a)\,=\,\sum_{n\geq 1}\phi_{n,\alpha}(a)\,=\,\sum_{n\geq 1}e^{-an}\sum_{k\geq 0}\frac{(\sqrt{4\alpha(1-\alpha)}an)^{2k}}{2^{2k}(k!)^2}.\end{equation}

In order to study the behavior of $\phi_\alpha(a)$ as a goes to 0 or $\infty$ , we take advantage of the simple inequalities $(2k)!\leq {2^{2k}(k!)^2}\leq (2k+1)!$ to bound $\phi_{n,\alpha}(a)$ from above and from below as follows:

\begin{align*}e^{-an}\frac{\sinh\big(\sqrt{4\alpha(1-\alpha)}an\big)}{ \sqrt{4\alpha(1-\alpha)}an }\,\leq\,\phi_{n,\alpha}(a)\,\leq\,e^{-an}\cosh\big( \sqrt{4\alpha(1-\alpha)}an\big).\end{align*}

From these inequalities, we deduce that, when $\alpha\neq 1/2$ , we have

(6.15) \begin{equation}\lim_{a\to 0} \phi_{\alpha}(a)\,=\,+\infty\,,{\qquad}\lim_{a\to +\infty} \phi_{\alpha}(a)\,=\,0\,. \end{equation}

Moreover, for any $a_0>0$ , we have

\begin{equation*}\forall a\geq a_0{\qquad}\phi_{n,\alpha}(a)\,\leq\,\exp\big({-a_0n}(1-\sqrt{4\alpha(1-\alpha)})\big).\end{equation*}

Thus, when $\alpha\neq 1/2$ , the series $\sum_{n\geq 1}\phi_{n,\alpha}(a)$ is uniformly convergent over $[a_0,+\infty[$ for any $a_0>0$ , and therefore the function $\phi_\alpha$ is continuous on $]0,+\infty[$ . We consider the equation

\begin{align*}\frac{1}{\sigma-1}\,=\,\phi_\alpha(a)\,,\end{align*}

and we denote by $a^1_c(\sigma,\alpha)$ (respectively $a^2_c(\sigma,\alpha)$ ) the smallest (respectively the largest) solution to this equation in $]0,+\infty[$ . These two real numbers are well-defined, thanks to the limits (6.15) and the continuity of $\phi_\alpha$ on $]0,+\infty[$ .

We use a strategy similar to that of the proof of Theorem 4. Suppose that $a\geq a_c^2(\sigma,\alpha)$ . The equation for $h=\lfloor\alpha\ell\rfloor$ in the system (6.5) yields

(6.16) \begin{equation}\frac{1}{\sigma-1}\,=\, \sum_{n\geq 1} \frac{1}{\lambda^n}M^n_H\big(\lfloor {\alpha\ell}\rfloor, \lfloor {\alpha\ell}\rfloor \big).\end{equation}

Let $\lambda^*$ be an accumulation point of $\lambda$ . If $\lambda^*>1$ , then, proceeding as in the proof of Theorem 4, we pass to the limit along a subsequence in Equation (6.16) to get

(6.17) \begin{equation}\frac{1}{\sigma-1}\,=\,\sum_{n\geq 1} \frac{1}{(\lambda^*)^n}\phi_{n,\alpha}(a)\,<\,\phi_\alpha(a),\end{equation}

but this contradicts the fact that $a\geq a^2_c(\sigma,\alpha)$ . Suppose that $a< a_c^1(\sigma,\alpha)$ . Let $\lambda^*$ be an accumulation point of $\lambda$ . Suppose that $\lambda^*=1$ . By Fatou’s lemma, we would have

(6.18) \begin{equation}\phi_\alpha(a)\,=\,\sum_{n\geq 1}\phi_{n,\alpha}(a)\,\,\leq\,\liminf_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } } \sum_{n\geq 1} \frac{1}{\lambda^n}M^n_H\big(\lfloor {\alpha\ell}\rfloor, \lfloor {\alpha\ell}\rfloor \big)\,=\,\frac{1}{\sigma-1},\end{equation}

where the last equality comes from (6.16). Yet the inequality (6.18) stands in contradiction to the fact that $a<a^1_c(\sigma,\alpha)$ . Therefore the infimum limit of $\lambda$ has to be strictly larger than 1.

7. Survival of the flattest

This final section is devoted to the proof of Theorem 6.

7.1. Lower bound on $\lambda_{\text{SP/PL}}(\delta,\sigma,\ell,q)$

The mean fitness $\lambda_{\text{SP/PL}}(\delta,\sigma,\ell,q)$ is also the Perron–Frobenius eigenvalue of the matrix $\smash{\big(f_{\text{SP/PL}}(u)M(u,v), {u,v\in E} \big)}$ . As such, it is a non-decreasing function of the entries of this matrix. Therefore

\begin{align*}\lambda_{\text{SP/PL}}(\delta,\sigma,\ell,q)\,\geq\,\max\big(\lambda_{\text{SP/PL}}(\delta,1,\ell,q),\lambda_{\text{SP/PL}}(1,\sigma,\ell,q)\big).\end{align*}

With the help of Theorems 3 and 4, we conclude that

\begin{align*}\liminf_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }\lambda_{\text{SP/PL}}(\delta,\sigma,\ell,q)\,\geq\,\max\big(\delta e^{-a},\lambda(a,\sigma)\big).\end{align*}

Let us prove next the reverse inequality.

7.2. Lumping

Let $(x(u), u\in \lbrace\,0,1\,\rbrace^\ell)$ be the solution to the quasispecies equations associated with $f_{\text{SP/PL}}$ . We are dealing with a fitness function which depends only on the Hamming classes. So we lump together the genotypes according to their Hamming classes and introduce the new variables y(h), $0\leq h\leq\ell$ , given by

\begin{align*}\forall h\in\{\,0,\dots,\ell\,\}{\qquad} y(h)\,=\,\sum_{u:H(u)=h}x(u).\end{align*}

We use the same computations that led us from (6.2) to (6.5) (these computations relied on the fact that the fitness function depended only on the Hamming classes) to conclude that the variables y(h), $0\leq h\leq\ell$ , satisfy

(7.1) \begin{equation}\forall h\in\{\,0,\dots,\ell\,\}{\qquad}y(h)\,=\,\sum_{0\leq k\leq\ell}y(k)\big(f_{\text{SP/PL}}(k)-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n(k,h),\end{equation}

where the matrix $M_H(k,h)$ is defined in (6.4). Again, with a slight abuse of notation, we denote by $f_{\text{SP/PL}}(k)$ the common value of the fitness function $f_{\text{SP/PL}}$ for all the sequences in the Hamming class H(k). In the rest of the proof, we write simply $\lambda(\ell,q)$ or even $\lambda$ instead of $\lambda_{\text{SP/PL}}(\delta,\sigma,\ell,q)$ . The value $\lambda$ and the variables y also satisfy

(7.2) \begin{equation}\lambda\,=\,\delta y(0)+\sigma y(\ell/2)+1-y(0)-y(\ell/2)\end{equation}

and

\begin{align*}\sum_{0\leq k\leq\ell}y(k)\,=\,1.\end{align*}

7.3. The two main equations

Let us rewrite the two equations of the system (7.1) corresponding to $h=0$ and $h=\ell/2$ :

(7.3) \begin{align}y({0}) &\,=\,y(0)\big(\delta-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n (0,0)\,+\,y\Big(\frac{\ell}{2}\Big)\big(\sigma-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n \Big(\frac{\ell}{2},{0}\Big)\,, \cr y\Big(\frac{\ell}{2}\Big) &\,=\,y(0)\big(\delta-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n \Big({0},\frac{\ell}{2}\Big)\,+\,y\Big(\frac{\ell}{2}\Big)\big(\sigma-1\big) \sum_{n\geq 1} \frac{1}{\lambda^n}M_H^n \Big(\frac{\ell}{2},\frac{\ell}{2}\Big).\end{align}

Let us examine the possible limits of $\lambda$ , y(0), and $y(\ell/2)$ along a subsequence in the long-chain regime. Suppose that, along a subsequence, the following convergences take place:

(7.4) \begin{equation}\lambda\to\lambda^*\,,\quad y(0)\to y_0^*\,,\quad y(\ell/2)\to y_1^*.\end{equation}

In the remainder of the argument, all the limits are taken along this subsequence. Passing to the limit in (7.2), we get

(7.5) \begin{equation}\lambda^*\,=\,\delta y_0^*+\sigma y_1^*+1-y_0^*-y_1^*.\end{equation}

Since $\lambda^*\geq\lambda(a,\sigma)>1$ , necessarily $y_0^*+y_1^*>0$ . Moreover, the following convergences take place:

(7.6) \begin{equation}M^n_H(0,0)\to e^{-na}\,,\quad M^n_H\Big(\frac {\ell}{2},\frac {\ell}{2}\Big)\to \phi_n(a)\,,\quad\end{equation}

while

(7.7) \begin{equation}M^n_H\Big(\frac {\ell}{2}, {0}\Big)\to 0\,,\quad M^n_H\Big({0}, \frac {\ell}{2}\Big)\to 0.\end{equation}

The convergences (7.6) have already been proved. The first one is obvious and the second one is the purpose of Lemma 3. To prove the convergences (7.7), we use the expression for $M_H$ computed in (6.4) and get

\begin{align*}M^n_H\Big(\frac {\ell}{2}, {0}\Big)&\,=\,\mathbb{P}\bigg( \text{Bin}\Big(\frac{\ell}{2}, \frac{1+\rho^n}{2}\Big) +\text{Bin}\Big(\frac{\ell}{2}, \frac{1-\rho^n}{2}\Big)\,=\,0\bigg)\cr&\,\leq\,\mathbb{P}\bigg( \text{Bin}\Big(\frac{\ell}{2}, \frac{1+\rho^n}{2}\Big)\,=\,0 \bigg)\,\leq\,\mathbb{P}\bigg( \text{Bin}\Big(\frac{\ell}{2}, \frac{1}{2}\Big) \,=\,0\bigg)\,,\cr M^n_H\Big({0}, \frac {\ell}{2}\Big)&\,=\, \mathbb{P}\bigg(\text{Bin}\Big(\ell, \frac{1-\rho^n}{2}\Big)\,=\,\frac{\ell}{2} \Bigg)\,\leq\, \mathbb{P}\bigg(\text{Bin}\Big(\ell, \frac{1}{2}\Big)\,=\,\frac{\ell}{2} \Bigg) .\end{align*}

It is well known that the right-hand quantities go to 0 as $\ell$ goes to $\infty$ .

7.4. Upper bound on $\boldsymbol\lambda_{\textbf{SP/PL}}(\boldsymbol\delta,\boldsymbol\sigma,\boldsymbol\ell,\boldsymbol{q})$

Proceeding as in the proof of Theorem 4 (see Lemma 3 and thereafter), we take advantage of the fact that $\lambda^*>1$ to interchange the limit and the four sums appearing in (7.3). So, passing to the limit in (7.3), and using the convergences (7.4), (7.6), (7.7), we get

(7.8) \begin{align}y_0^* &\,=\,y_0^*\big(\delta-1\big)\sum_{n\geq 1} \frac{e^{-na}}{(\lambda^*)^n}\,=\,y_0^*\frac{\delta-1}{{\lambda^*}{e^{a}}-1}\,, \cr y_1^*&\,=\,y_1^*\big(\sigma-1\big)\sum_{n\geq 1} \frac{\phi_n(a)}{(\lambda^*)^n}.\end{align}

If $y_0^*>0$ , then the first equation implies that $\lambda^*=\delta e^{-a}$ . If $y_1^*>0$ , then the second equation is equivalent to Equation (1.13). Therefore $\lambda^*$ has to be equal to the unique solution $\lambda(a,\sigma)$ of this equation. In conclusion, there are only two possible limits along a subsequence for $\lambda$ , namely $\delta e^{-a}$ and $\lambda(a,\sigma)$ ; thus

\begin{align*}\limsup_{ \genfrac{}{}{0pt}{1}{\ell\to\infty,\, q\to 0 } {{\ell q} \to a } }\lambda_{\text{SP/PL}}(\delta,\sigma,\ell,q)\,\leq\,\max\big(\delta e^{-a},\lambda(a,\sigma)\big).\end{align*}

The final conclusions of the theorem are obtained as a by-product of the previous argument. Indeed, if $\delta e^{-a}\neq\lambda(a,\sigma)$ , then it follows from (7.8) that we cannot have simultaneously $y_0^*>0$ and $y_1^*>0$ . This remark, in conjunction with the identity (7.5), yields the convergences stated at the end of the theorem.