{'fraction_non_alphanumeric': 0.03818976886235857, 'fraction_numerical': 0.010411605469563406, 'mean_word_length': 4.007368274711525, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Finite Population Model of Molecular Evolution: Theory and Computation

This paper is concerned with the evolution of haploid organisms that reproduce asexually. In a seminal piece of work, Eigen and coauthors proposed the quasispecies model in an attempt to understand such an evolutionary process. Their work has impacted antiviral treatment and vaccine design strategies. Yet, predictions of the quasispecies model are at best viewed as a guideline, primarily because it assumes an infinite population size, whereas realistic population sizes can be quite small. In this paper we consider a population genetics-based model aimed at understanding the evolution of such organisms with finite population sizes and present a rigorous study of the convergence and computational issues that arise therein. Our first result is structural and shows that, at any time during the evolution, as the population size tends to infinity, the distribution of genomes predicted by our model converges to that predicted by the quasispecies model. This justifies the continued use of the quasispecies model to derive guidelines for intervention. While the stationary state in the quasispecies model is readily obtained, due to the explosion of the state space in our model, exact computations are prohibitive. Our second set of results are computational in nature and address this issue. We derive conditions on the parameters of evolution under which our stochastic model mixes rapidly. Further, for a class of widely used fitness landscapes we give a fast deterministic algorithm which computes the stationary distribution of our model. These computational tools are expected to serve as a framework for the modeling of strategies for the deployment of mutagenic drugs.


Introduction
The rapid genomic evolution of viruses such as HIV has made the design of drugs and vaccines with lasting activity one of the most difficult challenges of our time. A novel intervention strategy which potentially outplays viruses in this evolutionary game was suggested by the pioneering work of Eigen and coworkers [Eig71,EMS89]. Eigen and coworkers considered the asexual evolution of a haploid organism and found that when the mutation (or evolutionary) rate was small, the organism survived as a collection of closely related yet distinct genomes together termed the quasispecies. Viral populations in infected individuals are known to exist as such quasispecies [LA10]. Remarkably, this quasispecies model predicted that when the mutation rate increased beyond a critical value, called the error threshold, the collection of genomes in the quasispecies ceased to be closely related; in fact, a completely random collection of genomic sequences was predicted to emerge. This transition with increasing mutation rate thus induced a severe loss of genetic information in the quasispecies and has been referred to as an error catastrophe. The generic antiviral drug ribavirin has been shown to act as a mutagen-an agent that induces mutations-against poliovirus and trigger a severe loss of viral infectivity in culture [CCA01]. This strategy of enhancing the viral mutation rate thus appears promising and particularly advantageous because it is unlikely to be susceptible to failure through viral evolution-driven development of drug resistance. Mutagenic drugs that attempt to induce an error catastrophe are thus being explored as a potential antiviral strategy [CCA01, GPLL + 05, ADL04], and one such drug for HIV is currently under clinical trials [MHH + 11].
The success of mutagenic strategies relies on accurate estimates of the error threshold of the pathogens under consideration. Notwithstanding the tremendous insights into viral evolution the quasispecies model provides, important gaps remain between the quasispecies model and the realistic evolution of viruses and other haploid asexual organisms. First, whereas the model assumes an infinite population size and, hence, adopts a deterministic approach, real populations are often small enough to lend themselves to substantial stochastic effects. For instance, the effective population size of HIV-1 in infected individuals is about 10 3 − 10 6 [KAB06,BSSD11], which is thought to underlie the strongly stochastic nature of HIV-1 evolution. Second, the model assumes a single-peak fitness landscape, where one genomic sequence is assumed to be the fittest and all other genomes are equally less fit. Realistic fitness landscapes are far more complex [HMC + 11]. There have been significant efforts in the last 30 years to close these gaps [Wil05]. While more general landscapes have been successfully considered in the quasispecies case [SH06], a rigorous treatment of the finite population case has remained elusive (see Section 4). Importantly, it still remains to be established whether the insights gained from the quasispecies model, such as the occurrence of an error catastrophe, translate to the more realistic, finite population case.
Here, we consider a finite population model of the asexual evolution of a haploid organism. Following standard population genetics-based descriptions [HC06], the model considers evolution in discrete, nonoverlapping generations. Within each generation, genomes undergo reproduction (R), selection (S), and mutation (M), yielding progeny genomes for the next generation. We analyze this RSM model formally and establish the following key results. We show that in the infinite population limit, the expected structure of the quasispecies predicted by the RSM model converges to that of the quasispecies model. Thus, insights from the quasispecies model may be translated to the finite population scenario. Indeed, we show further that the error threshold predicted by the RSM model also converges in the infinite population limit to that of the quasispecies model. Finite population models, such as the RSM model, appropriately tuned to mimic specific details, such as the fitness landscape, of the pathogens under consideration may thus be employed to obtain realistic estimates of the error threshold.
Unlike the quasispecies model, where the quasispecies is identified readily using black-box eigenvector determination algorithms, identifying the expected quasispecies of the RSM model is computationally prohibitive even for the smallest realistic genome and population sizes. Monte Carlo sampling techniques are therefore often resorted to [Wil05, AB05, BKP + 11, GD10,TBVD12]. Here, going beyond the ideas of the quasispecies model, we examine the mixing properties of the RSM model. We establish constraints on the model parameters under which the RSM model exhibits rapid mixing and therefore allows fast estimation of the expected quasispecies structure. Finally, we suggest an algorithm that uses the Markov Chain Monte Carlo paradigm to estimate the error threshold in the RSM model. Our study thus serves as a framework for elucidating quantitative guidelines for the modeling of intervention strategies that employ mutagenic drugs. The paper is organized as follows. In Section 2 we briefly describe the quasispecies model and the notion of the error threshold. In Section 3 we setup the finite population RSM model, present our main results and outline the techniques employed. In Section 4 we discuss our results in the context of previous studies and highlight open problems arising from work. Formal statements of our results are presented in Section 5. Detailed proofs are contained in the Appendix.

Preliminaries
We consider the evolution of a population of haploid organisms that reproduce asexually. In this evolutionary process the genome of each individual is modeled as a string of nucleotides. During reproduction, the genome is copied with possible mutations, which can be insertions, deletions, or substitutions. In applications, such as the modeling of viral evolution, it is often convenient to neglect insertions, deletions, and substitutions other than transitions (A to G or C to T , and vice versa). Under this assumption, a genome may, without loss of generality, be represented as a binary string. We thus represent a genome as an L-bit string σ = (σ 1 , σ 2 , . . . , σ L ), where σ i ∈ {0, 1}.
The fitness of a genome is then modeled in terms of its propensity to produce copies of itself. Specifically, the fitness of the genome σ is defined by the number of copies a σ of itself that it produces in one round of replication (also called one generation). However, during replication, each bit of each of the a σ offsprings is copied incorrectly with probability µ (called the error or mutation rate), thus potentially giving rise to an L-bit string different from σ . The fittest genome, also termed the master sequence, is without loss of generality assumed to be 0 = (0, . . . , 0) so that a 0 > a σ for all σ = 0.
The primary cause of the complexity and diversity in the evolution of such organisms is the variety of possible fitness landscapes, which a priori can be arbitrary functions from {0, 1} L to the set of nonnegative integers. Several special classes of fitness landscapes have been employed in the literature and we list the important ones below. We will assume that a σ ≥ 1. The case a σ = 0 for some σ 's has been used [WK93,TH07,GD10], and will be discussed in Section 4.

(General)
Here the only condition is that a σ ≥ 1.

(Class
Other landscapes such as the simpler additive or linear landscapes and more complex correlated landscapes have also been employed in the literature [BS93,Wie97,vNCM99].

The Quasispecies Model
Eigen and coworkers [Eig71,EMS89] gave the following differential equations for the time-evolution of the fractional population of the genome σ at time t, denoted by x σ (t): Here, is the probability that σ mutates to τ and d H (σ , τ) is the Hamming distance between σ and τ.Ā(t) is the average fitness ∑ σ a σ x σ (t) at time t. Defining A σ τ def = a σ when σ = τ and 0 otherwise, they showed that the vector of stationary frequencies, v σ µ def = lim t→∞ x σ (t), is the dominant right eigenvector of the value matrix QA at mutation rate µ such that v µ 1 def = ∑ σ v σ µ = 1. 1 The collection of genomes determined by this dominant eigenvector, which marks the culmination of the evolutionary process, is called the quasispecies. It is important to note that the vector v µ is independent of the starting population distribution.
We will mostly be concerned with the discrete time version of the quasispecies model. In the discrete time case, t = {0, 1, . . .}, denoting the fraction of genomes of type σ at time t by m σ t , Eigen's equations can be written as: (1) In vector notation, given the fractional population m t at time t, the fractional population m t+1 at time t + 1 is given by m t+1 = r(m t ), where the σ co-ordinate r σ of the vector valued function r is defined as Again, m t can be shown to converge to v µ irrespective of the starting population distribution as t goes to infinity. However, at any finite t, m t depends on the initial state m 0 .

The Error Threshold
With the single peak landscape, Eigen et al. observed empirically that there is a critical value µ c ≤ 0.5 for the mutation rate µ such that for µ µ c , the quasispecies is dominated by the master sequence, i.e. v 0 µ ≥ v σ µ ∀σ , whereas when 0.5 ≥ µ > µ c the quasispecies is dispersed approximately uniformly. The critical mutation rate µ c is called the error threshold because the uniform dispersal for µ > µ c implies a severe loss of representation in the quasispecies of the genetic information encoded by the master sequence. Evidently, this dispersal also decreases the mean fitness, A def = ∑ σ a σ v σ µ . We note here that despite the notion of the error threshold being widely recognized, no consensus exists on its definition; see Wilke [Wil05]. Since, in most cases, v µ will never become exactly the uniform distribution on {0, 1} L , it is clear that the goal is to find the smallest µ such that v µ is close to the uniform distribution on all genomes, i.e., the vector with every coordinate equal to 2 −L , which we denote by U. To define the error threshold we also need a function that measures closeness: e.g. v µ − U 1 , v µ − U ∞ or the difference in Shannon entropies of v µ and U, namely |H(v µ ) − L|. Hence, for a given distance function d which measures closeness of v µ and U, and a bound on closeness ε > 0, we define First note that at µ = 1 /2, the steady state vector v µ is exactly U. Hence, µ d c (ε) ≤ 1 /2 for all d, ε > 0. Second, note that changing the distance function d will change the error threshold quantitatively. Eigen and coworkers presented a heuristic argument that the error threshold should scale as 1 /L for the single-peak model without any rigorous proofs of its existence and without mentioning any closeness function.

A Finite Population Model and Our Main Results
In this section, we describe at an informal level the salient features of our work, which comprises a finite population model to capture molecular evolution, and our theoretical and computational results associated with it. We give a high-level technical overview of the methods used to prove our results in Section 3.3, while precise definitions and formal statements of our results appear in Section 5. Proofs have been moved to the Appendix due to considerations of space.

A Finite Population (RSM) Model
We consider the following stochastic discrete time finite population model of evolution which we call the RSM model. The parameters are the same as in the quasispecies model: the genome length L, the sequence space {0, 1} L , a per bit mutation rate µ and the fitness landscape {a σ } σ ∈{0,1} L with all a σ ≥ 1 and integers. At any time t, let N σ t be the number of genomes (a random variable) of type σ , and fix the total population ∑ σ N σ t to be N. This fixing of the population size to N at each step is the key distinction from the quasispecies model and is a new parameter. In each time step, the ensuing evolution is then described in terms of the following three steps.
1. (Reproduction) First, in the reproduction step, each genome σ produces a σ copies of itself, giving rise to an intermediate population 2. (Selection) Second, in the selection step, N genomes are chosen at random without replacement from this intermediate population of size I t , resulting in the selection of S σ t genomes of type σ where 3. (Mutation) Third, in the mutation step, each selected genome is mutated with probability µ per bit, giving rise to the next generation of N σ t+1 genomes of type σ , such that ∑ σ ∈{0,1} L N σ t+1 = N.
The starting state is denoted by N 0 which is typically ggiven by N 0 0 = N and N σ 0 = 0 for all σ = 0, but we will often not use this assumption. This RSM model is best viewed as a Markov chain where the state space is the set of functions f : {0, 1} 2 L → {0, 1, . . . , N} such that ∑ σ f (σ ) = N. Thus, the number of states of this Markov chain is N+2 L −1 N which is roughly N 2 L . It can be shown (see Fact 5.1) that for any 0 < µ < 1 the transition matrix of this Markov chain has a unique stationary vector, denoted by π. π is indexed by all f satisfying the property above and ∑ f π( f ) = 1, i.e., π is a probability distribution over the state space of the RSM Markov chain.
,1} L denote the random vector which captures the fractional population at time t. It can also be shown that lim t→∞ E [D t |D 0 ] exists and is independent of D 0 . We denote this limit as E [D ∞ ] and it can be computed from the stationary vector π as follows: Finally, we will subscript D t with parameters such as µ, N when we want to highlight dependence on them, e.g. D t,µ,N . The key questions of interest, especially given the fact that computing π would be prohibitive even for small values of L and N, are: We present theoretical results that address all of the above questions. Note that if the answer to the first question is that E [D t |D 0 ] converges to the prediction of the quasispecies model m t with the same starting states (m 0 = D 0 ), then it is important as we can leverage the significant understanding obtained from the study of the quasispecies model while incorporating stochastic effects with finite populations. For the second question, we need to first define a notion of the error threshold in the RSM model. We do so formally, given a distance function d.
where U is the uniform distribution over all genomes of length L.
What we want to understand is lim N→∞ µ d c (ε, N). Again, if we can show that the answer here is µ c (ε), then we can translate the insights from the quasispecies model to the RSM model.
For the third question, first note that if we want to estimate the error threshold, we need to be able to compute E [D ∞ ] . Secondly, we consider the standard notion of efficiency: the algorithm to estimate E [D ∞ ] should be polynomial in the input size. As we noted, the state space of the RSM Markov chain is prohibitively large and computing the stationary state is prohibitive. We employ the Markov Chain Monte Carlo method and run the RSM process for some time τ such that it is guaranteed that distribution from which D τ is drawn comes statistically close to the distribution from which D ∞ is drawn, irrespective of D 0 . Simulating each step of the random walk can be done efficiently. Hence, we are led to the question of bounding the mixing time of the RSM Markov chain: the smallest time the finite time distribution needs to come close to the steady state distribution for all starting configurations.
The issue of how the input is presented is also important and we briefly discuss it here. In one model, one can be given all a σ which would require bit length about ∑ σ log a σ and can, in principle, be as large or even larger than 2 L . Often, this is not the case and either the values a σ are given by a simple equation, or only some fixed number, say k 2 L of the values a σ are strictly bigger than 1. In the latter case the input has bit length roughly O(k log max σ a σ + log L + log 1 /µ). Another case for the input is when a σ are class invariant. In this case the input is of length roughly O(L log max σ a σ + log 1 /µ). We now proceed to summarize our results.

Our Results
We now give informal statements of our main results, before describing the mathematical techniques employed in the proofs of these results. The formal statements of the theorems described here appear in Section 5, after a discussion in Section 4, while the formal proofs are deferred to the Appendix.

Convergence of the Quasispecies and the RSM Model
Theorem 3.2 (Convergence of the RSM and Quasispecies Models). Fix a fitness landscape A with positive entries and a mutation transition matrix Q. Consider the RSM process started with the initial state D 0 and consider the evolution of the quasispecies model started with the initial state m 0 = D 0 . Then for any fixed time t, lim where m t is the fractional population vector at time t starting from m 0 predicted by the quasispecies model.
The theorem shows that in the infinite population limit, the stochastic fluctuations of the RSM process disappear, and the model converges to the quasispecies model. Informally, the main technical difficulty in proving the above theorem is to establish a result of the form lim N→∞ E[D t |D t−1 ] = lim N→∞ E [D t |D 0 ] with probability 1, which would establish convergence to the quasispecies model. The full proof is deferred to the Appendix. This convergence result allows us to show that for any distance function d, a finitary version of the error threshold, µ d c (ε, N) as defined above, converges to the error threshold µ d c (ε) of the quasispecies model, as the population size goes to infinity. These two results provide validation for the finite population RSM model by establishing that in the infinite population limit, the predictions from the RSM model converge to those of the quasispecies model. We now move on to problems concerning the mixing time and other computational issues of the RSM model.

Computational Results in the RSM Model
As noted before, a primary computational question in both the quasispecies model and the RSM model is the determination of the quasispecies, or the expected population profile at stationarity in the RSM model, which can then be used to estimate error thresholds (see Section 3.3 for an overview and Sections 5.3.2 and 5.3.3 for details). For the quasispecies model, a satisfactory solution to this problem is obtained via the observation that the quasispecies is the leading right eigenvector of the QA matrix. The QA matrix is of dimension 2 L ×2 L , and the above observation can thus be used to obtain efficient algorithms using black-box eigenvector finding algorithms for moderate values of L. In the case of class-invariant fitness landscapes, it is known [SS82] that one only needs to find the leading eigenvector of an (L + 1) × (L + 1) matrix.
However, similar approaches are not as effective for the RSM model. In this case, the stationary distribution is the leading eigenvector of the transition matrix M of the RSM process which is of dimension roughly N 2 L . Using ideas similar to those referred to above, one can reduce the running time for computing the stationary distribution to N O(L 2 ) . However, as mentioned before, in many applications, as in the case of HIV, for instance, where N ∼ 10 3 − 10 6 and L ∼ 10 4 , the problem is still computationally prohibitive. In these cases, one typically resorts to Monte Carlo simulations of the RSM process for estimating the population profile at stationarity [TBVD12], and thus we are led to considering the mixing time of the RSM process. The following theorem derives conditions on the parameters of evolution under which the RSM model mixes rapidly.
Theorem 3.4 (Mixing Time of the RSM Process). Given a fitness landscape A, mutation rate µ, the RSM process exhibits fast mixing if (1 − 2µ) max τ a τ min τ a τ L + 1 N < 1.
Having stated our results, we now highlight the techniques employed in the proofs.

Overview of Our Technical Contributions
As before, we will denote by D t the random variable of fractional populations after t steps of the RSM process, and by S t the random variable of the populations of genomes after the replication and selection steps of the (t + 1)-the step of the RSM process.
Our convergence result (Theorem 3.2). The starting point of the proof of our convergence result is to observe that E [D t+1 |D t ] has the same functional form r (as a function of D t ) as the evolution equation of the discrete time quasispecies model, with r as defined in Equation (2): Using the Lipschitz continuity of the evolution function r, we can then chain these concentration results inductively to show that the evolution of D t is tightly concentrated around the evolution of the discrete time quasispecies model, which allows us to show that E [D t ] converges to the quasispecies as N → ∞. To illustrate the ideas involved, we consider the case L = 1. Here the two genomes are {0, 1}. After the replication phase in the t-th step, there are a 0 D 0 t N copies of 0. For the i-th copy, let R i denote the indicator variable for this copy being selected in the selection phase, so that Since the R i 's are not independent, we cannot directly apply a Chernoff bound. However, since they are negatively correlated, one expects concentration to hold, and this can indeed be shown using the so-called method of bounded differences. The same reasoning works for S 1 t , and thus we get that given D t , the intermediate population S t after the replication and selection steps is concentrated around its expectation with high probability. We now look at the mutation step. Let M i be the indicator variable for the ith genome being 0 after the mutation step. We then have ND 0 t+1 = ∑ N i=1 M i . Since the M i 's are independent random variables, it can be shown using a Chernoff bound that given S t , The two steps can then can be combined to show that given The same reasoning works for D 1 t+1 . With some more work, this argument can be generalized to work for arbitrary L. The concentration guarantee we obtain is of the following form: there are quantities ε t and p t which are both o N (1) such that In the next step, we chain these stepwise bounds inductively in order to remove the conditioning and show that for all t ≤ t 0 , D t is concentrated around m t . An important component of the induction is the observation that r is Lipschitz continuous, which allows us to control the propagation of the errors ε t in each step. By the induction hypothesis, we have that |D t − m t | ≤ ε t with probability at least 1 − p t , where ε t and p t are both o N (1). Assuming the Lipschitz constant of r to be K, with probability at least 1 − p t . Applying the convergence result from the first step, we then have that with probability at least 1 − p t+1 = 1 − p t − p t , |D t+1 − m t+1 | ≤ ε t+1 = Kε t + ε t of m t+1 . The quantities p t , ε t for t ≤ t 0 can be chosen to be o N (1), which is sufficient to show the required convergence. The details appear in Appendix B.2. We now give an overview of the proofs of our computational results.
Computing the stationary distribution in the class invariant case (Theorem 3.3). Recall that the state space of the RSM Markov chain is roughly N 2 L . However, if the fitness function is class invariant, we can show that the number of distinct coordinates in the state space is about N L . To see this, first we define an equivalence relation on the states of the RSM Markov chain. We say that f , g, which are functions from have the same statistics for every Hamming class, i.e., for every 0 ≤ i ≤ N, Thus, the state space of the RSM Markov chain gets partitioned into about (N + 1) L+1 different classes. Then, due to the fact that the fitness function is class invariant, it can be shown that the transition probability of f to any other equivalence class is the same as that of g to the same class. Hence, one only needs to compute the transition probability from one equivalence class to another. This probability is a large binomial sum and one has to be careful in its computation and keep track of the number of bits required to represent each entry of this Markov chain over the equivalence classes. Once we have the transition matrix of this Markov chain, one can compute its largest eigenvector which corresponds to the stationary state. We show that, if one does this carefully, one can compute the eigenvector in time roughly N O(L 2 ) . The details appear in Section B.5.
Algorithm to compute the error threshold. Once we have the ability to either compute the stationary state of the RSM process or derive independent samples from its stationary state (which allows us to estimate the relative frequencies of the genomes at stationarity with a good precision by taking an average of the sampled states), the algorithm to estimate the error threshold is simple. The idea is to start with a small value ( 1 /L) of µ, and to estimate/compute the stationary distribution of the RSM process for the current value of µ. The algorithm then checks if the estimate of the stationary distribution is close to the uniform distribution on the genomes in the measure of closeness of one's choice. If so, it stops and outputs the current value of µ as an estimate of the error threshold. Else, it increases µ by a very small amount and repeats the above steps. In case direct computation of the stationary state of the RSM process is computationally prohibitive, independent samples from the stationary distribution of the RSM process are derived by simulating the RSM process up to its mixing time. The number of samples required can be estimated from a simple application of Chernoff bound on the random variable corresponding to the stationary state distribution of the RSM process. Hence, to establish bonds on the running time of the error threshold estimating algorithm, it is important to be able to bound the time it takes for the RSM process so that D t comes close to the stationary state, D ∞ . Our next result is towards this.
Mixing time result (Theorem 3.4). Since the stationary distribution of the RSM chain is not very well understood, it is not clear how to apply conductance-based geometric tools or the canonical paths method (see, for example, [JS88]) in order to prove the mixing time result. We are thus led to more combinatorial coupling based methods. Here one starts with an integer valued metric d on the state space of the Markov chain, and then one runs two copies X t and Y t of the chain. To show fast mixing, it is then sufficient to prove that X t and Y t can be coupled so that E [d(X t+1 ,Y t+1 )|(X t ,Y t )] ≤ α < 1. In general, defining a coupling can be tricky because one needs to carefully argue that the marginals of the coupling agree with the original Markov chain.
We define the coupling in two phases: the first phase includes the replication and selection steps and the second phase includes only the mutation step. We begin with the easier mutation step. Let I t and J t denote the state of the two RSM chains after the replication and selection steps. For most natural choices of the distance metric d, it is possible to couple the mutation step using the standard coupling for the random walk on the hypercube so that while coupling the replication and selection steps, since because of the global nature of the replication and selection steps, can become quite large. We control this increase by a careful choice of the metric d, and by appealing to the path coupling methods of Bubley and Dyer [BD97]. The path coupling theorem says that for integer valued d, it is sufficient to only for states X t and Y t satisfying d(X t ,Y t ) = 1 in order to establish fast mixing. Our coupling is then defined as follows. Fix a permutation of the N genomes in the chain X t : d(X t ,Y t ) is then the minimum over all possible permutations of the N genomes in Y t of the sum of the Hamming distances between the genomes at the same positions in the two permutations. The main technical step is to show that for this d, the replication and selection steps can be coupled in such a way that starting It is in this step that we use the form of the distance metric d crucially; the details of the coupling are somewhat technical and involve arguing carefully that the coupling is valid, and are given in Appendix B.4. We then combine this with the coupling for the mutation step described above to show contraction in the expected distance under the condition

Previous Work
The notion of the quasispecies and the existence of an error threshold were recognized first by Eigen and his coworkers in the 1970s and 1980s [Eig71,EMS89]. Translation of these ideas into intervention strategies requires overcoming two key limitations of the quasispecies model. First, the model assumes an infinite population size, whereas realistic population sizes can be quite small. With HIV, for instance, the effective population size is estimated to be ∼ 10 3 − 10 6 [KAB06,BSSD11]. Second, the theory assumes a singlepeak fitness landscape, whereas realistic landscapes can be far more complex [BCP + 04, HMC + 11]. Efforts over the last several decades have attempted to overcome these limitations of the quasispecies model [NS89, BS93, WH96, Wie97, AF98, SH06, TH07, PMnD10] (see Wilke [Wil05] for a recent review). The finite population case, however, has remained difficult to solve in full generality. Most studies resort to simulations or use approximate or heuristic approaches to describe the finite population case, and we discuss some of these here.
Nowak and Schuster [NS89] used a birth-death process to model the underlying evolution in finite populations and using simulations predicted that the error threshold scales as 1/ √ N. Their model, however, does not converge to the quasispecies model as N goes to infinity. Alves and Fontanari [AF98] present a model which employs a two-stage sampling with replacement in the selection process: first sampling uniformly from the population, and then sampling from the obtained sample with biases proportional to the fitness. They note, however, that sampling with replacement destroys the negative correlation between the selection of two individuals of the same species induced by the finite population constraint when selection is implemented using sampling without replacement. They find that the error threshold scales as 1/N. Further, they only analyze a heuristic deterministic approximation to their model, and do not consider rigorously the question of convergence of their original model to the quasispecies model. The closest to our convergence result is that by van Nimwegen et al. [vNCM99] who show convergence to a deterministic Eigen-like dynamics but again employ sampling with replacement and use special cases of additive fitness landscapes. With finite populations, there can be a significant statistical difference between sampling with and without replacement, the latter (which we employ) being more realistic. It is well known that as N → ∞ the difference between sampling with and without replacement shrinks, but then as we prove, so does the difference between our population genetics model and the quasispecies model. Their convergence proof has a similar structure as ours but we are able to use Chernoff bound type inequalities which are much stronger than the second moment inequalities used by them. Consequently, our convergence results are quantitatively stronger. We additionally prove convergence of the error threshold and fast mixing, questions not considered by [vNCM99]. More recently, Musso [Mus11] presented the transition matrix for sampling with replacement in the case L = 1 and also claimed convergence to the quasispecies model in the deterministic limit. No attempt, however, is made in [Mus11] to make this latter claim rigorous. Another class of studies relies on approximations and heuristics inspired from physics, and in particular statistical mechanics, to render the finite population case mathematically tractable (e.g., [BK98,SRA08,PMnD10]).
While previous studies have focused extensively on the fractional distribution of genomes at stationarity, little is known of the time to reach the stationary state. Campos and Fontanari [CF99] show that in the limit of infinitely large genome lengths (L → ∞) and population sizes (N → ∞) and with the single peak fitness landscape, the timescale associated with the decline of the master sequence is 1 /ln(qa) where q is the probability that a genome is replicated without error, and a is the relative fitness of the master sequence. Further, they show that with finite populations, this timescale is proportional to √ N. The mixing time when L and N are both finite and when the fitness landscape is more general than the single peak remains unknown. The latter mixing time has practical significance in the modeling of the action of mutagenic drugs, as it respresents the duration of therapy required to ensure completion of the transition to the error catastrophe. Our study presents conditions when the mixing is rapid and hence the transition to error catastrophe occurs quickly. Further, for computational studies that attempt to realize this transition in silico, our study presents an algorithm that allows efficient Monte Carlo sampling-based estimation of the error threshold.

Applications of the RSM Model
The motivation behind the RSM model and the algorithms discussed here is to get a basic framework for understanding the evolution of viruses of current interest such as HIV. Making concrete predictions relevant to the clinical setting requires super-imposing the specifics of the viruses of concern on the present framework. This often involves subtle modifications of the RSM process along with validation against data. For example, in related recent work, two of the authors and their co-workers applied the RSM model to mimic the within-host genomic evolution of HIV-1 [TBVD12]. It has been shown before [BSSD11] that these simulations quantitatively capture data of the evolution of viral genomic diversity in patients over extended durations (∼ 10 years) following infection and the approach is extended in [TBVD12] to estimate the error threshold of HIV-1. We envision that similar adaptation of our model will prove useful in elucidating the evolution and treatment guidelines for other asexual haploid organisms of interest.

Critique of the RSM Model
We note that our structural and computational results are independent of the nature of the fitness landscape so long as there are no lethal mutations (a σ = 0 for any σ ). Our model, however, does not consider lethal mutations. While letting some a σ be 0 does not affect the quasispecies model due to the constant rescaling involved, it introduces an absorbing state in the RSM Markov chain, thus making it non-ergodic, and causing the population to eventually decrease to zero. While Wilke and others [Wil05,WK93,TH07] have commented on the role of lethal mutations in extreme cases, establishing their full implications lies beyond the scope of the present paper. Although lethal mutations do occur, it turns out that in many important scenarios such as the evolution of HIV-1, a Hamming class invariant landscape without lethal mutations appears to capture key features of the underlying fitness interactions [BCP + 04], rendering our RSM framework applicable.
Finally, we note that our assumption of a fixed population size, N, is consistent with the widely accepted population genetics-based models of evolution, where a constant effective population size is employed to quantify the strength of stochastic effects [HC06]. Note that allowing N to vary with time (generations), does not increase the complexity in our model. The distinction between an infinite population model and a finite population model arises from the culling of the population in the latter model in order to maintain a finite population size. A fixed N or varying N will only result in different extents of culling in different generations, but will not change the overall structure of the model. The advantage in keeping N constant for our present study is that it allows easier examination of the convergence to the quasispecies model.

Open Problems
Our study of the quasispecies and RSM models has revealed several interesting and important problems. We list the main ones here.
Structure of the Quasispecies. Perhaps the most attractive feature of finite population models as opposed to the quasispecies model is that they can be used to study the effect of random genetic drift on inter-patient variations. Inter-patient variations in disease progression and response to treatments are known to be significant with HIV infection [NBS + 98, GKB + 05]. The collection of viral particles in an infected individual may be thought of as one realization of the random viral evolutionary process, and lim t→∞ Var [D σ t ] then provides an estimate of inter-patient variations in viral evolution due to the effect of the finite population size. Thus, in addition to the structure of the quasispecies in the finite population case, defined by the expected frequencies lim t→∞ E [D t ] when N < ∞, the variance of the frequencies lim t→∞ Var [D σ t ] as a function of the population size N is also an important quantity to be studied.
Error Threshold. In the quasispecies model with the single-peak fitness landscape, µ c has been found, without a rigorous proof, to be O(1/L), so that an error catastrophe occurs for µ 0.5 (e.g., see [EMS89]). Further, the transition is sharp, so that a small increase in µ from below to above µ c induces a dramatic change in the quasispecies structure. With other fitness landscapes, such as the multiplicative landscape, however, the quasispecies approaches the uniform distribution gradually as µ approaches 0.5 [WH96]. Further, lethal mutations, where a σ = 0 for some σ 's, appear to show the existence of an error threshold only if multiple mutations in a single replication are allowed [WK93,TH07,Wil05]. Thus, the conditions under which a sharp transition leading to an error catastrophe at µ c 0.5 would occur remain to be established. Second, the dependence of µ c on N remains to be identified. While some simulations suggest a 1/ √ N dependence [NS89,BS93], others find the dependence to go as 1/N [AF98]. As pointed out before, knowledge of µ c for finite N is important in the modeling of antiviral strategies based on mutagenic drugs.
Mixing Time. The main outstanding question here is to get a tight bound on the mixing time of the RSM Markov chain for a full range of evolutionary parameters. We notice that our result shows a good mixing time bound only under certain conditions on the parameters. Though we conjecture that the chain is rapidly mixing for other values of the parameters too, we believe that novel methods would be needed to extend our results in this direction. Apart from being useful in determining the time required for simulations to produce samples from the stationary distribution, the mixing time bounds also have biological significance. For example, when modeling the effect of a mutagenic drug under the RSM model, the convergence rate would models the minimum required duration of treatment before the error catastrophe occurs.

Preliminaries and Definitions
In this section we present rigorous statements of our results. Several definitions may be found repeated here in the interest of the readability of this section. We recall that genomes of length L are denoted by L-bit 0-1 strings. We will denote the Hamming distance between genomes σ and τ by d H (σ , τ), and the Hamming weight of a genome σ by w H (σ ). A population is defined as a multiset of genomes of the same length. While discussing the RSM model, we will fix the size of the population to be N.
The Markov Chain for the RSM Model. We will denote the evolution of the RSM process using a timeindexed sequence of vector valued random variables (N t ) ∞ t=0 . The entries of N t are indexed by genomes σ , and the entry N σ t denotes the number of genomes of type σ at time t. At every time t, ∑ σ ∈{0,1} L N σ t = N. The random variables D σ t def = N σ t /N denote the fractional population of the genome σ at time t.

Reproduction
Step and the Fitness Landscape In the reproduction step, each genome σ produces a σ copies of itself, so that the number of genomes of type σ after this step is I σ t def = a σ N σ t , and the total number of genomes is I t def = ∑ σ a σ N σ t . The matrix A defined by A σ σ def = a σ and A σ τ def = 0 for σ = τ is called the fitness landscape. The fitness landscape is said to be class-invariant if a σ depends only on the Hamming weight of σ . By a slight abuse of notation, we will denote by a i the fitness of all genomes with Hamming weight i in the class invariant case.

Selection and Mutation Steps and the Mutation Rate
In the selection step, N genomes are sampled without replacement from the genomes obtained after the reproduction step. In the mutation step, each bit of each of the N genomes obtained after the selection step is flipped with a probability µ, called the mutation rate. The mutation transition matrix Q defined by Q σ τ τ) gives the probability that a genome of type σ mutates to one of type τ in the mutation step. The RSM process as described above is a Markov chain on the state space of functions f : {0, 1} L −→ N, satisfying ∑ σ ∈{0,1} L f (σ ) = N. The transition matrix M of this chain is described in Section A.
Fact 5.1. When the mutation rate µ ∈ (0, 1) and a σ > 0 for all σ , the Markov chain M corresponding to the RSM process is ergodic, and hence has a unique stationary distribution. This is a simple consequence of the fact that µ and A are positive. See Section B.1 for a proof. Thus, we can consider the projected Markov chain M w with state space Notice that |Ω w | = N+L L . Also, if π w is the stationary distribution of M w , then by the projection property This property implies that the expected populations for every Hamming class of genomes and, hence, the expected average fitness at stationarity, are the same for M w and M .
Mixing Time. We will denote by π the stationary distribution of the RSM process, and let N ∞ be a random variable distributed according to π. We know that the distributions of the random variables N t converge in total variation distance (and hence in distribution) to π, due to the ergodicity of the RSM process. We fix our notation for mixing times in this section.
Definition 5.3. The total variation distance between two probability distributions D 1 and D 2 on the sample space Ω is defined by D 1 − D 2 TV = max A⊆Ω |D 1 (A) − D 2 (A)|.
Definition 5.4. Fix a Markov chain N on a state space S. We define For 0 ≤ ε ≤ 1/2, the mixing time of N is defined by Notice that by the projection property, the mixing time of the projected RSM chain M w is at most the mixing time of the original RSM chain M .
Error Thresholds. In the following definition, we specifically emphasize the dependence of the random variables N t , D t on µ, N by denoting them as N t,µ,N and D t,µ,N . We denote by D ∞,µ,N a version of D t,µ,N distributed according to the stationary distribution of the RSM process, and by U the uniform distribution over genomes. Given a distance function d, one can define the error threshold with respect to d as follows.
where U is the uniform distribution over all genomes of length L.

Convergence to the Quasispecies Model
Our first main result is that the RSM model converges to the quasispecies model.
where m t 0 is the state of evolution of the quasispecies model at time t 0 starting from m 0 .
The proof of the above theorem is relegated to Section B.2. As a corollary to the theorem above, one can show that there is convergence of a finitary version of the error threshold µ h c (ε, N) to the error threshold µ h c (ε) for the quasispecies model, as the population size goes to infinity. Formally, we have the following: Corollary 5.7. Fix a mutation rate µ ≤ 1/2 and an error parameter ε. For every δ > 0, there exists a time t 0 > 0 such that for t > t 0 , one can find an N δ ,t such that for N > N δ ,t , , and d h E D t,µ,N , U ≤ ε + δ when µ = µ h c (ε), Here we use the subscripts µ and N to emphasize the dependence of the distribution of D t on µ and N.
Although we will prove our results for the error threshold in terms of the average Hamming distance, it is easy to translate our results to other common dispersal measures as described in Section B.1.1. The proof of the above Corollary follows easily from Theorem 5.6 and is given in Section B.3. We note here that extending the above corollary to get convergence of finite population error thresholds depends upon proving a strengthened version of our convergence result (Theorem 5.6), which we leave as an open problem. In fact, on the basis of simulation results, we conjecture that for fixed ε, µ h c (ε, N) monotonically increases to µ h c (ε).

Mixing Time Bounds on the RSM Process
We give a coupling argument in Section B.4 which allows us to prove the following theorem.

Computing the Stationary Distribution in the Class Invariant Case
Theorem 5.9. For every A which is class invariant, and µ which can be represented using b bits, there is an algorithm running in time T given by which computes π w for the Markov chain M w described above. For fixed L, T = O(N O(L 2 ) ).
The proof of the above theorem appears in Section B.5, and is based on the projected RSM process discussed in Section 5.1. The above theorem immediately gives a grid-search based algorithm that given a grid resolution δ and ε > 0, outputs a approximation µ 0 to the error threshold in time T · 1 /2δ such that µ 0 ≥ µ h c (ε, N) and d h (D ∞,µ 0 −δ , U) > ε. We now consider Markov Chain Monte Carlo based grid-search methods.

Markov Chain Monte Carlo Methods
The general strategy for Monte Carlo based grid search methods for determining error thresholds is described in the algorithm ERRORTHRESHOLD in Figure 1 in the Appendix. We will denote the mixing time τ mix (ε) for parameters L, N, A and µ as τ(L, N, A, µ, ε). We consider the projected chain M w described above which contains enough information to compute the average Hamming weight, and whose state can be maintained as a tuple in {0, 1, . . . , N} L+1 .
The proof of the above theorem appears in Section B.6, where we also point out some technical subtleties about the definition of error thresholds.

A Starting State and Transition Matrix of the RSM Markov Chain
As stated before, the RSM Markov chain starts with the "fittest" possible population with all the weight concentrated on the master sequence, so that N M 0 = N and N σ 0 = 0 for all σ = M. We now proceed to set up some notation for writing out the transition matrix M .
Definition A.1. (Multivariate Geometric distribution). Let g(σ ) denote the number of genomes of type σ in an urn. Consider the process of choosing, without replacement, N genomes from this urn. Then P hyp (g → f ; N) denotes the probability of obtaining f (σ ) genomes of type σ for each σ . We have, Definition A.2. (Multivariate Binomial Distribution). Let f (σ ) denote the number of genomes of type σ . Consider a stochastic process in which each genome of type σ independently mutates into a genome τ (possibly equal to σ ) with probability Q(σ , τ). We denote by P bin ( f → D; Q) the probability that D(σ , τ) genomes of type σ mutate to type τ under this process. We have We can now write the entries of M . Proof Sketch of Fact 5.1. When µ ∈ (0, 1) and a σ > 0 for all σ , it can be verified easily that this chain is irreducible and also has a non-zero self-loop probability at every point in the state space. Thus, the chain is ergodic and hence by the Fundamental theorem of Markov chains, has a unique stationary distribution to which it converges as t → ∞.
We now give a proof of Lemma 5.2.
Proof of Lemma 5.2. We will show that under class invariance, we can project the RSM Markov chain so that its state space consists of equivalence classes indexed by functions h : We will find it convenient to consider the reproduction and selection phases separately from the mutation phase, show that the projection described above can be done for both of them, and then combine the two results using the following general fact about projected Markov chains, the proof of which we include for completeness.
Fact B.1. Let P and R be the transition kernels of two Markov chains on the same state space Ω, and let S denote the composition PR of the two chains. Suppose that there is a partition of Ω into equivalence classes Ω , such that for any f ≡ f , and any equivalence class Then, we also have S( f , [g]) = S( f ([g]), for all f , f and g as described above.
Proof. The proof is by direct computation. We have, Just as in the derivation of equation (6)  In order to use the last fact, we now decompose the matrix of the RSM process into the following two Markov chains on Ω: 1. The Reproduce-Select Chain. We denote the transition matrix of this chain as P, such that P( f , g) is the probability of obtaining the state g starting from state f after the reproduction and selection phases. Notice that P( f , g) = P hyp (A f → g; N) .

Assume that A is class invariant and let
Since R( f ,

B.1.1 Relationships between Error Thresholds
We first define error thresholds according to various dispersal measures.
Definition B.2 (Error Thresholds). Let ε ≥ 0, and U be the uniform distribution over the set of genomes.
Here H denotes the Shannon entropy, using the base e. Definition B.3 (Error Threshold for the Quasispecies Model). Let µ ∈ (0, 1) and let v µ denote the the stationary expected fraction vector with 1 norm 1. The error thresholds are defined as follows.
where H denotes the Shannon entropy, using the base e.
Our results are mostly stated in terms of the error threshold µ h c . However, we now describe how the different definitions above are related to each other. The following inequalities relate the different distance measures that we have considered. The first of these follows from the definition of the 1 norm, while the second is the well known Pinsker's inequality relating the 1 norm to the entropy.
This gives us the following relationship between the error-thresholds: Using the fact that the distributions involved are defined over a state space of size 2 L , we can show the following weak converse to inequality (10): This gives us a further relationship between the error thresholds: However, we notice that one cannot in general close the loop in inequalities (9) and (10) (and hence in inequalities (11)) above. To see this, consider for example the following two distributions P and Q for L > 1.
1. P: puts total weight 1 − ε on weight 1 strings and weight ε on the string 0, so that the average Hamming weight is 1 − ε.
The average Hamming weight in both cases is 1 − ε, so that in that metric, the distance between P and Q is zero. However, the total variation distance between P and Q is at least 1 − ε.

B.2 Proof of Theorem 5.6
In the rest of this section, we will use the following concentration inequalities about the multivariate hypergeometric distribution: Fact B.4. Consider the hypergeometric distribution P hyp (g → f ; N) defined in equation (5) above. Let D σ be the random variable denoting the fraction of genomes of type σ which are drawn in the process starting with g(τ) genomes of each type τ. We then have: 1. E [D σ ] = g(σ ) g,1 .
2. The following concentration inequality holds for ε ≥ 0: Similarly for the multivariate binomial distribution, we have the following: Fact B.5. Consider N genomes with f (σ ) genomes of each type σ . Let D σ be the random variable denoting the fraction of genomes of type σ after a mutation step. Then 2. The following concentration inequality holds for ε ≥ 0: 4 is a consequence of Azuma's inequality, and a proof can be found in the book by Dubhashi and Panconesi [DP09]. Fact B.5 is essentially a restatement of the Chernoff-Hoeffding bound. Combining the above bounds, we can deduce the following concentration inequality for each step of the RSM process: Lemma B.6. Consider a state N t of the RSM process. We then have with r σ as defined in equation (2).
Proof. For ease of notation, let g σ = a σ N σ t . Let I σ be the random variable denoting the fraction of genomes of type σ left after the selection step. Thus, we have Using a union bound over all genome types with the concentration inequality in Fact B.4, we get that with probability at least 1 − 2 L+1 exp −ε 2 1 N conditioned on N t , we have We denote the above event by E . Now, we consider the concentration of D t conditioned on I. Using a union bound over all genome types along with the concentration inequality in Fact B.5, we get with probability at least 1 − 2 L+1 exp −2ε 2 2 N conditioned on I, we have We denote the above event by F . With probability at least 1 − 2 2L+1 (exp −ε 2 1 N + exp −2ε 2 2 N ), conditioned on N t , both E and F occur, and then we have, for all σ , which is what we sought to prove.
Before proceeding, we need the following lemma: Lemma B.7. Fix a fitness landscape A with positive entries and a mutation transition matrix Q with µ < 1/2. The functions r τ defined in equation (2) are Lipschitz with Lipschitz constant on the set of probability distributions over genomes.
Proof. For a probability distribution x over genomes, we have where the last line follows by noticing that fact that for all x and all σ , min σ ,τ Q σ τ ≤ r σ (x) ≤ max σ ,τ Q σ τ , and min σ τ Q σ τ = µ L , while max σ τ Q σ τ = (1 − µ) L . Thus, by the mean value theorem, for any probability distributions x and y over genomes, we get |r(x) − r(y)| ≤ K ||x − y|| 1 .
Proof (of Theorem 5.6). Fix a time t 0 . In the rest of the proof, we drop the conditioning on the initial state being concentrated on the master sequence for ease of notation. We will prove the following claim by induction for 0 ≤ t ≤ t 0 : We now proceed to prove Claim B.8.
Proof of Claim B.8. At t = 0, we can set l σ t = u σ t = m σ t , and p t = 0. By the definition of the starting state, this satisfies the conditions claimed in the claim. Now suppose that we have shown that with probability 1 − p t , we have D σ t ∈ [l σ t , u σ t ] for all σ . We call the latter event E t . Recall that Notice that m σ t+1 ∈ l σ t+1 , u σ t+1 . Also, because of the Lipschitz condition on the function r σ shown in Lemma B.7, we have u σ t+1 − l σ t+1 ≤ 2 L Ks t . Now, we condition on the event E t defined above, and in this case, we have , for all σ . Choose ε(N) = N −1/3 = o N (1), and set l σ t+1 = l σ t+1 − ε/2, u σ t+1 = u σ t+1 + ε/2. Using the concentration result quoted in Lemma B.6, we get that conditioned on E t , with probability at least 1 − p(N) where p(N) = exp(−Ω(N 1/3 )) = o N (1), D σ t+1 ∈ l σ t+1 , u σ t+1 , for all σ . Now, we saw above that E t occurs with probability at least 1 − p t . Hence, by a union bound, we get that with probability at least 1 − p t+1 , where p t+1 = p t + p(N), D σ t+1 ∈ l σ t+1 , u σ t+1 , for all σ . This proves the induction hypothesis, except that we need to make sure that s t , p t are o N (1). We first consider s t . From above, we have the following recurrence for s t :

B.3 Proof of Corollary 5.7
We begin by noticing that for 0 < µ < 1/2, we can choose a time t 0 such that for t > t 0 , the state m t of the quasispecies model satisfies where v is the unique stationary vector of the quasispecies model. Now fix t > t 0 . Since the distance function d h is continuous, Theorem 5.6 allows us to choose an N δ such that for N > N δ , Combining equations (14) and (15), we get , and, and when µ = µ h c (ε), d h E D t,µ,N , U ≤ ε + δ when µ = µ h c (ε).

B.4 Proof of Theorem 5.8
Before proving Theorem 5.8, we first set up some notation for the coupling argument.
Definition B.9. A coupling of two probability distributions D 1 and D 2 is a pair of random variables (X,Y ) defined on a single probability space such that the marginal distribution of X is D 1 and the marginal distribution of Y is D 2 .
Definition B.10. A coupling of Markov chains with transition matrix M is defined to be a process (X t ,Y t ) ∞ t=0 with the property that both (X t ) and (Y t ) are Markov chains with transition matrix M , although the two chains may possibly have different starting distributions.
Any coupling of Markov chains with transition matrix M can be modified so that the two chains stay together at all times after their first simultaneous visit to a single state: more precisely, such that if X s = Y s , then X t = Y t for t ≥ s. In the following, we only consider such couplings. The following well known facts are the basis of coupling based methods for proving mixing time bounds.
Theorem B.11. Let {(X t ,Y t )} be a coupling satisfying the definition above for which X 0 = α and Y 0 = β .
Let τ couple be the first time the chains meet: Lemma B.12 (Coupling Lemma). Let X,Y be random variables defined on a finite sample space Ω and let C be any coupling of C and Y. Then Definition B.13. Let d : Ω × Ω −→ R ≥0 be a distance metric on the state space Ω of the two Markov chains {X t } t and {Y t } t . Suppose C is a coupling such that for every t ≥ 0, for every starting distributions X 0 ,Y 0 , then we call C a (θ , d) coupling. Note that this implies that where D = max σ ,τ∈Ω d(σ , τ).
Fix a integer valued distance function d. Let let {X t } t be a realization of the Markov chain starting from X 0 and Y t be another realization starting from the stationary distribution π of the Markov chain. If C is a (θ , d) coupling, then it follows from the Coupling Lemma that This implies that the mixing time τ mix (ε) = O log (D/ε) log (1/θ ) when θ < 1.

B.4.1 A Coupling for the RSM Process
The coupling C we will construct will have two independent parts C = (C S , C M ), C S for the Reproduce-Select phase and C M for the mutation part. We first describe the somewhat simpler mutation coupling.
The mutation coupling follows the following algorithm, with M set to be the permutation which achieves the minimum in the above definition. Lemma B.14. For µ ≤ 1/2, C M is a ((1 − 2µ), d match ) coupling.
Proof. To prove that C M is a valid coupling one just needs to note that if r is distributed uniformly at random in [0, 1] then so is 1 − r. Hence, for i = 1, . . . , N and j = 1, . . . , L, each bit σ i ( j) (respectively, ν i ( j))) flips with probability exactly µ. Further, these flips are independent by construction.
To prove that C M is a ((1 − 2µ), d match ) coupling, let X t ,Y t be the states of the two Markov chains with By definition of d match , there is some matching M which achieves d. Without loss of generality assume that M is identity, i.e., M (i) = i, for all 1 ≤ i ≤ N.
Proof. By construction d match (X,Y ) ≥ 0 with equality if and only if X = Y . Now consider states X = in the state space Ω. Let α and β be permutations of [N] such that We will use the following general path coupling result of Bubley and Dyer [BD97] to define the coupling C S . Theorem B.16 (Path Coupling [BD97]). Consider a Markov chain M on state space Ω and a distance function d on Ω such that d (x, x) = 0 for all x ∈ Ω. Consider a connected undirected graph G on Ω such that the length of each edge {x, y}, if present in G, is d(x, y), and let d be the shortest path metric on G. Suppose there exists a coupling C for M such that for some α < 1, and all X t ,Y t ∈ Ω which are adjacent in G, If every edge of G is a shortest path under the metric d described above, then the coupling C can be extended to a coupling C such that for all X t ,Y t ∈ Ω.
We will first show now that the path metric resulting from an application of the above theorem to d match is d match itself, since this is crucial for composing the C S coupling with the coupling C M described above. Proof. For brevity we will denote d match (·, ·) by d. Notice that since each edge is of length 1, it is also a shortest path by construction. Since d is a metric, we also have d (X,Y ) ≥ d(X,Y ) for all X,Y ∈ Ω. We now proceed by induction to show that d(X,Y ) ≥ d (X,Y ) for all X,Y ∈ Ω. Notice that when d(X,Y ) = 1, this is true by definition of d . Now, suppose that d(X,Y ) ≤ k − 1 implies d (X,Y ) ≤ d(X,Y ), and consider the case d(X,Y ) = k > 1. We claim that there exists a Z such that d(X, Z) ≤ k − 1 and d(Y, Z) ≤ 1. The existence of such a Z implies using the induction hypothesis that It only remains to construct such a Z.
. Without loss of generality, we may assume that d( Since d(X,Y ) = k > 1, there exists a j such that d H (σ j , ν j ) ≥ 1. Let s be a string obtained by flipping a single bit of ν j such that Claims about the Coupling. Suppose we find a coupling C for the selection phase such that when d match (X t ,Y t ) = 1, then the intermediate states I(X t ) and I(Y t ) satisfy then using Theorem B.16 and the coupling for the mutation phase described above, we get This will give us fast mixing as long as α(1 − 2µ) < 1.
We now describe such a coupling for the selection process, for general X t and Y t , which we will analyze only in the simple but sufficient case when d match (X t ,Y t ) = 1. Suppose that X t and Y t contain, respectively, n x σ and n y σ genomes of type σ . After reproduction, the number of genomes of type σ in the two chains is a σ n x σ and a σ n y σ respectively. Let the total number of genomes be M x = ∑ σ a σ n x σ and M y = ∑ σ a σ n y σ respectively. Without loss of generality let us assume that M x ≥ M y , and set M = M x . We now construct a bag of M balls as follows. For each σ , the bag has exactly a σ min n x σ , n y σ balls with label (σ , (x, y)). If n x σ ≥ n y σ , then the bag has exactly a σ (n x σ − n y σ ) balls with label (σ , x), otherwise it has exactly a σ (n y σ − n x σ ) balls with label (σ , y). Thus the total number of balls in the bag is M. We now take a random permutation of the M balls, and take the intermediate state I X (respectively, I Y ) to be the multiset of genomes given by the first N balls carrying the label (x, y) or x (respectively, (x, y) or y). Notice that a ball carrying a label (x, y) can contribute a genome to both I X and I Y .
Claim B.18. The above coupling is a valid Markovian coupling for the selection phase of the RSM chain.
Proof. Notice that sampling without replacement a objects from a set of b objects is equivalent to taking the first a elements from a uniform random permutation of the b objects. Also note that given a subset S of a set of b objects, and a uniform random permutation α over the b objects, the restriction of α to the elements of S is a uniformly random permutation of the elements of S. Now consider the set of M labeled balls constructed above, and define S X (respectively, S Y ) to be the set of balls carrying a (x, y) or x (respectively, x(x, y) or y) label. By the observations above, see that the set I X (respectively, I Y ) has the same distribution as if it was sampled without replacement from S X (respectively, S Y ). This proves the claim.
Lemma B.19. Suppose d match (X t ,Y t ) = 1. Then under the above coupling, we have Proof. For brevity, we will denote X t and Y t as and d match (·, ·) by d. Since d(X,Y ) = 1, we can assume without loss of generality that σ i = ν i for i > 1, and σ 1 and ν 1 differ in exactly one bit.
We now consider the coupling described above, and let S be the set of balls with label (x, y). Notice that |S| = ∑ 1<i≤N a σ i . Notice that |S| ≥ (N − 1) min σ a σ . We assume without loss of generality that a σ 1 > a ν 1 , so that the number of balls M in the bag is |S| + a σ 1 . Consider a random permutation α of the M balls, and let I be the (random) multiset of balls with label (x, y) occurring in the first N positions of α. Notice that |I X ∩ I Y | ≥ |I|. We observe that if the intersection of I X and I Y (seen as multisets) is of size at least |I|, then d match (I X , I Y ) ≤ L(N − |I|), hence we have Now, we have Plugging this into equation (16), we get which establishes our claim.
Using the above discussion, we get fast mixing under the condition that (1 − 2µ) 1 1 − 1 N max σ a σ min σ a σ L < 1.
We note here that the main difficulty in designing coupling arguments for the RSM process is the distance expansion property of the reproduction and selection phases. Given two states of the RSM process, the reproduction phase amplifies the distance between them, and the nature of the selection phase tends to keep this distance intact. In this setting, the chain can be described as a noisy random walk on a boolean hypercube, and our bound reflects the intuition that when the noise is small, the fast mixing property of the hypercube should able to enforce fast mixing of the RSM chain. We reemphasize that we consider that achieving a better understanding of the mixing properties of the RSM walk, in terms of both upper and lower bounds appears to be an interesting and challenging open problem.

B.5 Proof of Theorem 5.9
We now proceed to prove Theorem 5.9. We will use notation similar to that used in the proof of Lemma 5.2 in Appendix B.1. By a slight abuse of notation, we will use P and R for the transition matrices of the projected versions of the chains described in the proof of Lemma 5.2 above. We note that for any numbers M and We start with an estimation of the time required to compute the entries of P. Using the above observation, and the form of entries of P, we have the following observation for P.
Observation B.21. Each entry of P is of size at mostÕ(N max a σ ), and can be computed in N max a σ operations over integers of sizeÕ(N max a σ ).
We now proceed to estimate the complexity of computing entries of the matrix R. We first get a bound on the time required to pre-compute the (L + 1) × (L + 1) matrix Q defined in the description of the mutation chain in Appendix B.1.
Observation B.22. Let b be the number of bits required to represent µ. We have Hence, each entry of Q is of size at mostÕ(bL) and can be pre-computed in time O(L 2 ) operations over numbers of sizeÕ(bL).
Proof of Theorem 5.9. Notice that the number of terms in the sum in equation (8)  .
Notice that R is of dimension at most S × S. Computing the products of all the Q i j 's in each of these terms takes N multiplications on numbers of size at mostÕ(bL), and hence produces a number of sizeÕ(NbL) in time at mostÕ(NbL). Computing the products of all the multinomial coefficients in each of the terms takes O(N) multiplications on integers of size log N, thus producing an integer of size O(N log N) in timẽ O (N log N). The total size of each entry of R is thus at most s 1 =Õ (log G + N log N + NbL), and the entry can be computed in time t 1 =Õ(Gs 1 ). All the entries of R can thus be computed in time S 2 t 1 .
Notice that M w = PR, and given the above estimates on the sizes of the entries of P and R and the times required to compute their entries, each entry of M w is of size at most s 2 = O(log S + s 1 + N max a σ ), and can be computed in timeÕ(Ss 2 ), given the entries of P and R. Thus, the time taken for the computation of M w , including the computation of R and P and the pre-computation of Q is T 1 =Õ bL 3 + S 3 s 2 + S 2 Gs 1 and each entry of M w is of size s 2 =Õ(NbL + N max a σ + L 2 ).
The time required to compute an exact solution of M w π w = π w with the restriction ||π w || 1 = 1 using Gaussian elimination is thus of the orderÕ(S 3 s 2 ). Comparing this with T 1 , we get that the total running time isÕ(T 1 ), which is what we sought to prove.
In particular, if we run s independent realizations of the chain up to the time τ(L, N, A, µ, ε/2), and let D i denote the sample obtained by the ith realization, then Proof of Theorem 5.10. Consider the random variables Z s in any of the at most c iterations in ERRORTHRESH-OLD. By the choice of ε and s, we get from the above lemma and a union bound that all the c different variables Z s that we get across all the iterations of the loop satisfy with probability at least 1 − δ 2 . In the rest of the proof, we condition on the above event occurring.
Since the average Hamming distance can be at most L, and we are running the projected chains till τ(L, N, A, µ, δ 1 /(2L 2 )) for each µ, we get that for every µ considered by the algorithm d h (π w,µ , U) − d h (Z s , U) ≤ δ 1 /2.
We therefore deduce that when the algorithm outputs a µ, d h (π w,µ , U) ≤ ε + δ 1 /2, using the conditioning on the event in equation (17). Similarly, by noticing that the algorithm had d h (Z s , U) ≥ ε for µ − δ , we get that d h (π w,µ−δ , U) ≥ ε − δ 1 /2. The estimate of the running time follows by noticing that assuming bits with bias µ can be sampled in O(log(1/µ)) time, it takes time O(NL max a σ log(1/µ)) to simulate each step of the M w chain. The bound on the error probability follows from the conditioning used on the validity of equation (17).
We comment briefly on a technical point about the definition of the error threshold that has been used in the literature (and that we use too). With this definition, there might exist a µ satisfying 1 /2 > µ > µ h c (ε) such that d h (v, U) > ε. An analogous condition might hold in the finite population case too. If we could preclude the occurrence of such anomalous behavior of error-thresholds, we would be able to improve the guarantee on the output µ 0 of the algorithm ERRORTHRESHOLD to be of the form µ h c (ε +δ 1 , N) ≤ µ 0 ≤ µ h c (ε −δ 1 )+δ . We observe, however, that somewhat surprisingly, even for the simpler case of the quasispecies model, to the best of our knowledge, no attempts have been made to rigorously prove that such anomalous behavior cannot occur. We leave the resolution of this point for the finite population case as an open problem.