{'fraction_non_alphanumeric': 0.0461195054945055, 'fraction_numerical': 0.011881868131868131, 'mean_word_length': 3.9091368846931895, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Equilibrium state of molecular breeding

We investigate the equilibrium state of the model of Peng, \textit{et al.} for molecular breeding. In the model, a population of DNA sequences is successively culled by removing the sequences with the lowest binding affinity to a particular target sequence. The remaining sequences are then amplified to restore the original population size, undergoing some degree of point-substitution of nucleotides in the process. Working in the infinite population size limit, we derive an equation for the equilibrium distribution of binding affinity, here modeled by the number of matches to the target sequence. The equation is then solved approximately in the limit of large sequence length, in the three regimes of strong, intermediate and weak selection. The approximate solutions are verified via comparison to exact numerical results.


I. INTRODUCTION
Recent advances in molecular biology have enabled breeding methodologies to be applied at the molecular level to develop novel proteins and DNA sequences with particular desired characteristics [1]. The basic idea is simple: in each generation, a number of molecules are selected from the population according to some criteria of desirability; they are then diversified (via point mutation and/or recombination) and then amplified back to the original population size. This methodology is worthy of study not only for its immediate practical importance but also as a forum for exploring general issues in evolution and population dynamics in a uniquely controllable setting.
Recently, motivated in particular by the experiment of Dubertret, et al. [2] on the in vitro evolution of DNA sequences, Peng, et al. (PGHL) [3] introduced a simple model of molecular evolution. In the experiment, at each round a specified percentage of the population was selected to survive based on its binding affinity to a particular target, the lac-repressor protein. This was accomplished by coating a beaker with the lac-repressor protein into which the DNA was introduced. The beaker was washed out, so that only the most strongly bound sequences remained. The degree of selection could be directly controlled by the strength of the washing procedure. Amplification, accompanied by mutation, was subsequently performed by multiple stages of PCR. The model of PGHL incorporates these basic processes of selection, mutation and reproduction. Each sequence is modeled as a string of L nucleotides, each taken from the set of A = 4 "letters", {A,G,C,T}. The binding affinity of the sequence is taken to be proportional to the number of nucleotide matches between it and the target sequence. The fraction (1 − φ) with the lowest affinities are dropped from the population. The population is then amplified back to its original level, where each daughter sequence in the next round is chosen to be descended from one of the surviving sequences, modified by point mutations at a rate of µ 0 per nucleotide.
PGHL performed simulations on their model and found that, starting from an initial population with low affinity, the average affinity of the population rose, eventually reaching equilibrium. They then constructed a mean-field (i.e. infinite population) treatment of the problem. They showed that the average affinity relaxes to its equilibrium value exponentially in time, with a time constant of Aµ 0 /(A − 1), and verified this via simulation. However, the calculation of the equilibrium value itself was less successful. In their treatment, this value was found to depend on a parameter which was introduced in an ad-hoc manner into their theory, and which needed to be fit from simulation. In this paper, we reexamine the equilibrium problem, locating the source of the difficulty of the PGHL treatment in an inappropriate continuum approximation. In turns out that the problem is fairly subtle, and in fact requires different treatment in the case of very weak, very strong and intermediate selection. The range of the validity of the intermediate selection increases with L, so that for sufficiently large L, the average affinity is given for any φ by the intermediate result.
The plan of the paper is as follows. In Section II, we lay out the general formalism and basic equations that we need to solve. In Sec. III, we solve the strong selection problem, working in the small µ 0 limit. In Sec. IV, we treat the case of intermediate strength selection, also for small µ 0 . In Sec. V, we generalize the intermediate strength solution to the case of arbitrary µ 0 , obtaining the central result of the paper, the average density of mismatches for any φ and µ 0 , given that L is sufficiently large. In Section VI, we tackle the weak selection limit, reverting back to small µ. In Sec. VII, we return again to the case of intermediate selection, calculating the next order corrections. We discuss the connection to the PGHL treatment in Sec. VIII. We conclude in Sec. IX with some observations.

II. BASIC FORMALISM
We are interested in this paper in calculating the equilibrium distribution of binding affinity in the population of sequences. Since the dynamics only depends on the binding affinity, determined by the number of matches between the sequence and the target, we can characterize each sequence by this quantity, or equivalently, the number of mismatches, 0 ≤ k ≤ L, in the length L sequence. We will work in the infinite population limit, so that fluctuations are irrelevant, and focus on the master equation for the normalized mismatch distribution, P k . We need to track the changes induced in this distribution by each of the two processes, selection and mutation; amplification in and of itself having no effect on the distribution.
Selection of the fraction φ of the highest affinity sequences throws out those states with the highest k's, keeping the lower k states. Thus, after selection, some set of states 0 ≤ k ≤ n are retained in their entirety, along with some fraction α of the k = n + 1 state. The renormalized distribution after selection, P s , is given by where the equation implicitly determines n and α.
If the mutation rate is sufficiently small, we need only consider single nucleotide mutation events, so that only transitions k → k ± 1 are allowed. We will take up the case of larger mutation rates in Sec. VII. The distribution after mutation, P µ , is given by where µ ≡ µ 0 L is the genomic mutation rate, ν ≡ 1/(A − 1),nu ≡ 1 − ν, and we have set P s −1 ≡ 0. Mutation repopulates the k = n + 2 state, which was emptied by selection. At equilibrium, the post-mutation distribution P µ k must reproduce the initial distribution P k , giving us a closed set of equations for P k . These equations involve only the finite set of states 0 ≤ k ≤ n + 2. In fact, P n+2 is not involved in any of the equations, since its initial value is zeroed out by selection. Thus, the equilibrium equations close on the n + 2 quantities P k , 0 ≤ k ≤ n + 1, with P n+2 being slaved to the solution for P n+1 . This system is (0 ≤ k ≤ n + 1) with the filling factor f k = 1 for 0 ≤ k ≤ n, f n+1 = α and f −1 = f n+2 = 0. These equations for the P 's, with a given n and α are linear and homogeneous, with φ playing the role of an eigenvalue condition for the existence of a nontrivial solution. As α increases for a given n, φ increases, implying weaker selection, till α reaches its limiting value of 1. At this point, we increase n and start α again at 0. Thus the solution is piecewise analytic in φ. This piecewise behavior is seen clearly in Fig. (1), where the numerical solution for L = 50, µ = 0.1, A = 2 is presented. For simplicity here, we will only treat the case of α = 1 − , where the last selected level is in fact completely selected. Writing N = n + 2 and setting P N = 0 for convenience (its true value being determined afterwards), we can write the equations in their final form with the boundary conditions P −1 = P N = 0, where we have introduced the scaled eigenvalue After the solution of this system, the P 's can be normalized via φ = N −1 k=0 P k , and then P N can be constructed from P N −1 via We see that the system equilibrium is controlled by the single combined parameter χ, rather than the two parameters φ and µ separately. Increasing φ or µ both act to increase χ, and as we shall see, increase the average mismatch. Since the rate of approach to equilibrium is governed solely by µ, one can achieve the same equilibrium in less time by increasing µ and the selection strength (thereby decreasing φ) simultaneously.
We see that there is a special role to the value φ = 1 − µ, where χ = 0. This is the limit of maximal selection, since at this point only the k = 0 is selected. Lower values of φ just result in a partial selection of the k = 0 state, which changes nothing since, after amplification, the result is the same. Thus, for φ ≤ 1 − µ the solution is P 0 = 1 − µ, P 1 = µ, and the mean mismatch k = µ. The other limit of χ is given by the weak selection limit φ → 1 − , where χ also approaches 1.
The system of equations Eq. (5) is not analytically solvable in general, but simplifies in the limit of large L, which we now proceed to study.

III. STRONG SELECTION, N ≪ L
We first consider the case of strong selection, where the number of selected states, N , is much less than the sequence length L. The system Eq. (5) can then be written as This homogeneous set of equations has a solution when the determinant, D N vanishes, providing an eigenvalue condition on χ. Expanding in minors about the last column gives a recursion relation Dropping theν (N −1) L in the first term and the N −2 L in the second term of the right side of Eq. (9) (because N ≪ L) gives the equationD whereD N = (L/ν) . We are thus interested in the zeroes of He N (y). In fact, we have to take the maximum root of He N (y), which we call y * N , because any other solution give rise to negative P n 's. Using Maple, for example, to compute the y * N 's, we obtain our approximate solution for N ≪ L: We now can computek, the average number of mismatches. Going back to Eq. (8), one can see that P 1 ∼ √ LP 0 , P 2 ∼ √ LP 1 etc., so that P N −1 is the single dominant term and is thus 1 to leading order (P N is of order µ, and so is also irrelevant in our small µ limit). This gives usk ≈ N − 1, so that This is presented in Fig. 1, along with the exact numerical results for the case L = 50, for the case A = 2 (ν = 1). As we noted earlier, the numerical results are piecewise continuous, due to the partially filled levels, with the transitions between levels apparent as cusps in the curve. We see that the agreement is quite good for the smaller N , and gets worse as N increases, as expected.
It is interesting to investigate the behavior of y * N for large N . We remember that the N 'th Hermite polynomial, times a Gaussian, is the solution of the Schroedinger equation for the harmonic oscillator, with ω = 1/2 and energy level E N = (N + 1 2 )ω, ( = m = 1): The zeroes of He N are the zeros of the wave function ψ. For large N , the maximal zero is near the classical turning point, where the energy equals the potential, so that To achieve a more accurate approximation, we expand the potential around the turning point, writing x = 2 √ N + a(z − z 0 ), and dropping the quadratic term gives an Airy equation The first zero of the Airy function is at z = −2.338, so that . The next order correction due to the dropped quadratic piece can be seen to be of order N −2/3 . This is an excellent approximation to y * N even for N = 2, as can be seen in Fig. 1, where this approximation is used to achieve a completely analytic expression for φ(k).
The upshot of all this is that in the first region the scaling variables are y = L/ν(φ − (1 − µ))/µ andk. Plotting the results for different L in these variables demonstrates the collapse beautifully, see Fig. 2.
The only question left to answer is how far the first region extends. We see from Fig. 1 that for L = 50 it covers roughly 70% of the range of φ above the threshold for single level selection, φ = 1 − µ, which is surprisingly high. The explanation of this depends on the results of the next section, with the intermediate range asymptotics. We anticipate the results by noting that the leading order correction to N (φ) is, for A = 2, δN ∼ N 2 /L. Thus, the first-region results start failing for δN/N ≈ 0.05, or N /L ≈ 0.05. For large L, χ ∼ 2 N/L, we expect that for large L the first-region results are trustworthy out to χ ≈ 2 √ .05 ≈ 0.45, i.e.; 45% of the region of interest. For smaller L, if fact things are better as we have seen. It can be shown that a more accurate criterion is (N − 2.338N 1/3 ) 2 /(N L) ≈ 0.05. Applying this to L = 50 gives N ≈ 10 which compares quite well with the numerical results in Fig. 1. For ν = 1, the first region is much narrower, though still a finite fraction of the allowed range of χ, independent of L. The first region approximation here fails when χ ≈ 0.1ν/ν, which for ab = 4, (ν = 1/3), reads χ ≈ .05, very much smaller. We can see this in Fig. 3, where again our approximation is over-conservative for smaller L.
When N gets too large, of order L, we can no longer neglect theν (N −1) L and N −2 L terms in the recursion relation for the determinants. It is more convenient now to consider the P n 's directly, using the difference equation 5. Motivated by the result that, in the first region, the P k 's grow rapidly with k, we assume that here too the same rule applies: The relevant P 's are those with k ∼ N . Writing N = αL, we thus assume that k ∼ αL. Plugging this into Eq. (5), we get the difference equation where There are 3 options for λ + and λ − : 1. They can both be real, but this is incompatible with P −1 = P N = 0. 2. They can both be complex, but then we get an oscillating solution, which is also not acceptable, since the P 's must be non-negative. 3. The degenerate solution, where λ + = λ − ≡ λ, and the solution that satisfies P N = 0 is This solution indeed decays rapidly on a scale of order 1 as n decreases from N , justifying our initial assumption. It also means that the change in the solution from one k to the next is also of order 1, and thus the difference equation can not be approximated by a differential equation, as PGHL attempted. The condition for degeneracy of the λ's is The last task is to calculatek. Since the P n 's decay exponentially (with an L-independent width) as n decreases, to leading order in 1/L,k = N .
Since here N = αL, there is a finite probability, q =k/L of a mismatch per nucleotide in the intermediate selection region.
Combining results, we have A comparison of the formula in Eq. (23) and numerical results is given in Fig. (4) for A = 2, and in Fig. (5) for A = 4 (ν = 1/3). We see that our formula overlaps that of the first region for α ≪ 1, so that χ ≈ 2 √ να = 2 νN/L, which is indeed the first region result for large N . The leading correction in the overlap can be calculated, leading to the results discussed above in Sec. III. The second region solution has q increasing from 0 at χ = 0 (φ = 1 − µ) to the fully random value of 1/(1 + ν) = (A − 1)/A at χ = 1 (φ = 1), so that this second region result covers the entire range of q and φ. Indeed, for fixed φ, this result is asymptotically correct for sufficiently large L. At small χ, our result only breaks down when N is not large, which only happens for χ ∼ O(L −1/2 ). Similarly, it breaks down α is too close to 1/(1 + ν), since then λ goes to unity and the relevant k's are no longer concentrated in a region around k = N . We shall see in Sec. VI, that this occurs for 1/(1 + ν) − α of order L −1/2 , so that λ is of order √ L. Thus, except for small regions of φ near the very weak and very strong limits, the size of which vanish in the large L limit, our result Eq. (23) is reliable.

V. INTERMEDIATE SELECTION, GENERAL µ
Up to now we have only considered the case of small mutation rate, so that only single transitions had to be considered. We now take up the case of large µ, dealing with multiple transitions in a single round. Again, we are working in the intermediate selection region where only the states with k ≈ N = αL are relevant. We take the probability of a single site mutation µ 0 to be small, but µ = µ 0 L to not necessarily be small, so that we may use the Poisson distribution for the number of up and down mutations, considered separately. The mean number of up (increasing mismatch) mutations is µ(1 − α) and the mean number of down mutations is µνα. Consider the case, for example, where after the mutation phase you wind up at the same level where you began. This could result, of course, from having no mutations of any kind, which has a probability e −µ(1−α) · e µνα = e −µ+µαν . In addition, however, one can wind up back where one started by having two mutations, one up and one down. The probability for this is (e −µ+µαν ) · µ(1 − α) · µνα. One can also return to the original state after 4, 6, . . . mutations. The total probability of remaining in the original state is then where I 0 is the modified Bessel function. Henceforth, we shall write 2µ αν(1 − α) ≡ t. Proceeding likewise, the probability of mutating down from state n + j (j > 0) to state n is and up from n − j to n is Putting this all together yields the recursion relation (I −j = I j ) The characteristic polynomial (here of infinite order) is This sum can be performed (see, e.g. [4], Eq. 9.6.33), yielding The condition for degeneracy is ∂F ∂λ = 0, yielding λ = 1−α αν , exactly as in the small µ case! Plugging this into Eq. (26) gives In the limit of small µ, φ ≈ 1 + µ 2 αν(1 − α) + α(ν − 1 which is of course what we got before. As with small µ, the mismatch density q is just α to leading order, so that Turning this around gives A comparison between Eq. (30) and numerical results for ν = 1 3 and µ = 1 is shown in Fig. (6). The formula for general µ can be seen to reproduce exactly that for small µ with the replacement (1 − φ) by − ln(φ), with the control parameter in general being − ln(φ)/µ, so that stronger mutation can always be exactly compensated by stronger selection.

A. Determining χ(N )
We return to Equation (5) describing the equilibrium state. We saw from the second region solution that what N approaches L/(1 + ν), the width of the distribution diverges and one can no longer consider k as a constant. However, the fact that the distribution is very broad (of order √ L as we will see) means that it changes slowly with k, which allows us to expand P k+1 and P k−1 in a Taylor series, yielding: where x is a continuous variable replacing the discrete variable k in Eq. (5). We can eliminate L by going over to the scaled variables yielding to leading order ν 1 + ν P ′′ y + ((1 + ν)y) P ′ y + (δχ + 1 + ν) P y = 0.
We again transform into a Schroedinger equation, via then similarity transformation P y = f y g y with Rescaling z = (1 + ν)y/ √ ν, we obtain (yet again) the Schroedinger equation for a harmonic oscillator, in which m = = 1, ω = 1/2, with the energy term reading The boundary conditions, P −1 = P N = 0, are equivalent to infinite walls at z = − L/ν and z = ((1 + ν)N − L)/ √ νL. Because the factor g z decays rapidly as z decreases away from 0, we can replace the left boundary condition by P (z = −∞) = 0. at the cost of an exponentially small error. Thus, for a given χ, or equivalently E, as we come in from y = −∞ we want to find the first z = z * for which f * z = 0, which then determines N through Thus, as we noted in the beginning, N is a distance of order √ L from L/(1 + ν) in the third region. There are of course a set of E's for which we know z * analytically. If E = (n + 1 2 )ω for some whole number n, then the solution of the Schroedinger equation is a Gaussian times the Hermite polynomial, He n , and so z * in this case is just the largest (negative) root of the He n . It is remarkable that, as in the first region, the solution is determined by the zeroes of the Hermite polynomial. However, whereas there the quantum level n increased with φ, here it is just the opposite. Now, the root z * decreases with n, so that as χ decreases from 1, N decreases, as we expect. At n = 1, the wave function has a single root at the origin, so z * vanishes, and so N = L/(1 + ν) at this point. For the case of general E, the solution for f is a parabolic cylinder function, and the root z * has to be determined numerically. In particular, for E < 3 2 ω, z * is positive, so N > L/(1 + ν). As the ground state wave function has no zeros, z * must go off to ∞ as E approaches the ground state energy of E = 1 2 ω. In fact, the case E = 1 2 ω is precisely the limiting case χ = 1 where there is no selection and so N = L. The divergence we obtain is consistent with with the fact that the exact z * is not order 1, but rather of order √ L at this point, and is a result of our approximation of the left boundary condition. Thus, the scaling variables in the third region are E and z * , which are related to the physical variables χ and N via Eqs. (37) and (38). The results for the numerical calculation of z * vs. E is presented in Fig. 7, together with the exact solution for integer n in terms of the zeros of the Hermite polynomial and the asymptotic formula for large E derived in the section on the first region:

B. Findingx(N )
Using the relation between the variables x and z above, we will calculate the averagex viā where z * is as above the value of z for which f z = 0. Using the numerical solution for f for various E's allows us to construct the scaling functionz(E), the scaledk, (i.e. (1+ν)k−L √ νL ) vs. the scaled selection parameter, E = 1/4 + (1 − χ)L/2(1 + ν) as shown in Fig (8).
We can generate an approximation toz for large E. Here as we saw, z * is large and negative and so the Gaussian factor in the integral restricts the relevant range of z to those near z * , which as we saw in the Region 1 calculation, is itself near the classical turning point z 0 ≡ −2 √ 2E. Proceeding as in Sec. III, we introduce z = z 0 + Az, A = (−z 0 /2) − 1 3 , and requiring f (−∞) = 0 yields We do not worry about normalization since we explicitly divide by the normalization integral. This leads tō wherez * = 2.338 is the first zero of f . Changing the integration variablez to s =z * −z yields: Because of the fast decay of the exponential factor, (|z 0 | ≫ 1), the main contribution for the integral in (43) is around s = 0 This allows us to expand Ai(s −z * ) in (43) in a Taylor series: Ai(s −z * ) ≈ sAi ′ (−z * ) + . . .. This gives: On the scale of Fig. (8), the difference between this formula and the exact numerical calculation is not visible.

VII. NEXT-ORDER CORRECTION TO INTERMEDIATE SELECTION
In order to better understand the origin of the degeneracy condition we invoked to solve the second region problem, it is worthwhile to explore the next-order correction. In addition, we will be able to understand the relatively poor accuracy of our result forq(φ). We will for simplicity restrict ourselves in this section to the case A = 2. In the second region, we could not replace the difference equation by a differential equation, since the change in P from step to step was of order 1. We will fix this by rescaling P by the geometrical series behavior we found, by putting Q n = Pn λ n , where λ is as determined above, λ = (1 − α)/α, so that the differential equation: becomes We now assume the Q's vary smoothly and expand Q n−1 and Q n+1 in a Taylor series, (writing n = αL + y), and dropping all terms of order L −1 , though not y/L) yielding where As usual, we transform into a Schroedinger equation by a similarity transformation, Q = f g, with: giving where we have again dropped a term of order L −1 . We see that the coefficient of f in Eq. (49) vanishes at y = 0 if χ is taken equal to χ 0 , the zeroth order solution. Thus the zeroth order solution corresponds to choosing the "energy" −χ so that the classical turning point is at y = 0, or n = N . However, the true boundary condition is that f vanish at n = N , which happens some small distance (relative to the length scale for changes of the potential L) behind the turning point. Thus, we must choose χ slightly less than χ 0 in order to move the turning point to negative y. Assuming y/L and δχ ≡ χ − χ 0 are small, we get an Airy equation: The solution of this equation that decays for negative y is The requirement that the first zero of f lie at y = 0 implies A comparison of the formula in (52) and the exact numerical results is given in Fig. (9). Remembering that χ + δχ = φ−1+µ µ yields: where N = αL. Now, let's remember we got for the second region at zeroth order P n ∝ (N − n)λ n . We can calculaten as follows:n = αL n=0 nP n αL n=0 P n . (54) with N = αL, and we get, up to exponentially small terms so we get a correction forn − N . Combining this correction with the correction for φ in (53), We get a more accurate graph of n as a function of φ.

VIII. COMPARISON TO PGHL
Before concluding, we wish to make some comments on the connection of this analysis to that attempted by PGHL. As we have noted, their approach was doomed by the assumption of a continuum limit for the distribution function. This is true, as we have seen, only in the weak selection region, which is valid only in a very limit range of φ near 1, of order µ/L. Even there, their answer as stated is incorrect, due to their neglect of the dependence of the effective diffusion constant on the local mismatch. If one uses for the diffusion constant that obtained at the average mismatch, one obtains an implicit equation for the average mismatch at equilibrium, which agrees with our results for the overlap region between the second and third regions. Neglecting this dependence forced PGHL to introduce a parameter τ which they fitted numerically to in essence restore the correct diffusion constant. Achieving agreement with the asymptotic behavior for weak selection would have dictated a value of τ = 2 in the infinite L limit. Fitting τ at the finite value L = 170 gave them instead a value of τ = 2.77. Actually, their formula neglecting the mismatch dependence of the diffusion constant and using τ actually results in better agreement with the exact results, this agreement is the fortuitous outcome of a partial cancellation of the errors induced both by the mismatch dependence of the diffusion constant and the continuum approximation. At strong selection, their approach in any case is unreliable. This situation is summarized in Fig. (12), where we show the results of our theory together with the stated predictions of PGHL, along with the (partially) corrected form of their theory taking account of the mismatch dependence of the the diffusion constant. case they considered. We see that the partially corrected form of PGHL agrees near φ = 1 and subsequently systematically deviates from the simulations. The τ = 2.77 original formula numerically does better overall, but is nowhere correct. Nevertheless, the basic mechanism responsible for selecting the equilibrium state, namely the degeneracy condition discussed in Sec. IV, was correctly identified by PGHL. This degeneracy condition is intimately connected with the marginally stable nature of the dynamic problem [5].

IX. CONCLUSIONS
We have seen that the steady-state equation is quite rich in its behavior, with different behavior in each of the three separate regimes of the selection parameter φ. As a summary, we present in Fig. 13 a graph with all three approximations displayed, as well as the exact numerical solution for N = 200, µ = 0.1, A = 2. We see that, as already noted, the strong selection approximation covers roughly the first half of the nontrivial range of φ, while the second-order intermediate selection result works essentially from the smallest φ = 1 − µ up to φ of 0.99. The weak selection theory works in this case already from φ = 0.97.
A few qualitative features stand out from the solution and are worth pointing out. First, the average mismatch is generically a finite fraction of the sequence length, except for very strong selection, extremely near the threshold for selecting only the best state. Second, the distribution of mismatches is very skewed toward high mismatch. The most probable state is very near the bottom of the barrel of the states that survive at all. The high affinity states are exponentially rare in the equilibrium distribution. Signs are this behavior are perhaps discernible in the experimental data of Ref. [2]. Last is the exact tradeoff between selection strength and mutation rate, due to the existence of a single lumped control parameter. This opens the possibility for optimizing the breeding process by increasing the mutation rate to speed the approach to equilibrium, at no cost in the quality of the final distribution of affinities. There is of course a limit to this imposed by the finite population size, which we have not considered at all in this work. If the selection criteria is too rigid, the finite population size will eventually begin to play a role.