1　Models

The Crow-Kimura quasispecies model was proposed to study the Drosphilla’s evolution before the Eigen quasispecies model, which was given some consideration in Chapter 6 in the textbook by Crow and Kimura (1970)^[2]. Just like the Eigen model, the Crow-Kimura model was traditionally formulated in the language of chemical kinetics and can describe the basic processes of mutation and selection for molecules with self-replicating and information decoding functions, for instance, DNA and RNA. In the Crow-Kimura model, different individuals are represented by macromolecular sequences with length N. If one only considers the purine and pyrimidine each site of a sequence is a binary variable, and the i-th sequence can be denoted by Si=(s1i,s2i,...,sNi) , where sn= -1 or 1 represents the purine or the pyrimidine. The total number of all possible sequences is 2 ^N . In the population, the sequence without mutation is called the wild sequence or reference sequence and the other sequences are the mutant sequences. The evolution dynamics of the sequences is given by

Here, xi and ri are the relative concentration and fitness of the i-th sequence, respectively. The latter reflects the replication rate of the sequence. A simple single peak fitness landscape is adopted in the present paper, in which r1=A0 , ri=Ai=A1 1<i≤2N and for all the mutant sequences they have the same fitness. The factor ∑j=12Nrjxj is the mean fitness of the population (or the dilution flux) which keeps the total concentration constant. mij represents the mutation rate from Sj to Si per unit period of time. We adopt the specification for the mutation rate given by Baake et al. ^[8].

where μ is the mutation rate for each site, hi, j represents the Hamming distance between sequences Si and Sj , h(i, j)=(N-∑k=1Nskiskj)/2 , which is the number of the different sites between the two sequences. The above coupled set of equations governs the evolution of the master sequence and mutant sequences, with which one may discuss the dynamics and steady state properties of the population. The coupled set of equations contains 2 ^N equations, which precludes any direct calculation even for moderate values of N. Nevertheless, we may classify all possible sequences according to the Hamming distance from the master sequence. The sequences with the same Hamming distance form a class and one has therefore N+1 classes. The classification simplifies the problem greatly. In addition, the transformation defined by the Eq. (3)^[8] is used in the present paper which turns the coupled set of equations into a set of linear equations. The set of linear equations can be written in the matrix form and its coefficient matrix determines the evolution dynamics of the classes (the population). The relative concentration for each class in its steady state is uniquely determined by the right eigenvector of the largest eigenvalue of the coefficient matrix^[5,28]. The steady state is the state when time goes to infinite, which is an asymptotic state thereof. In the present paper, we only focus on the steady state properties of the relative concentrations of the classes.

In the mutation randomized Crow-Kimura model, we regard the site mutation rate as a Gaussian distributed variable. Its probability density distribution is given by

Here, u is a variable, u¯ and ω2 denote the averaged value and variance of the random variable. In order to compare with the deterministic Crow-Kimura model, u¯ is identical to the mutation rate μ used in the deterministic Crow-Kimura model. The fluctuation strength of the random variable is defined as d=ω/u¯ , representing the relative width of the Gaussian distributed random variable. As a result of the randomization of the site mutation rate, the relative concentration of each class becomes a random variable and is characterized by its ensemble averaged value and the deviation from the mean value. We are only interested in the steady state properties of the randomized concentration.

2　Results and discussions

2.1　Characteristics of the error threshold in the mutation randomized Crow-Kimura model

In the numerical simulations, we set the sequence length N=20, the fitness of the wild class A0=10 and the fitness of all the mutant classes Ai=A1=1 , so a single peak fitness is used here. For the deterministic Crow-Kimura model, the relative concentration of each class at its steady state versus the site mutation rate ( μ ) is calculated (Figure 1) in which number 0 represents the wild class, and numbers 1, 2, 3… denote the mutant classes 1,2,3…. When the site mutation rate( μ ) is small, the quasispecies distribution appears, namely, the mutant classes gather around the wild class which peaks there. As μ increases, the relative concentration of the wild class decreases, and the relative concentrations of the mutant classes increase. If μ goes beyond the error threshold, all classes have the same negligible concentration. The error threshold is located at μ0 =0.49, which is a sharp point similar to a phase transition in physics.

Fig. 1 Relative concentrations of the classes in their steady states based on the deterministic Crow-Kimura model versus the site mutation rate (μ)

NOTE: The error threshold is located at μ0 =0.49. Number 0 in the figure represents the master class,and numbers 1,2,3… denote different mutant classes.

For the mutation randomized Crow-Kimura model, the site mutant rate is replaced by a Gaussian distributed random variable. We take ten thousand (or more) random samples and perform the ensemble average. The ensemble averaged relative concentrations of the classes at their steady states versus the site mutation rate μ are computed (Figure 2) where the fluctuation strength d=0.05. One can see that when the fluctuation strength is small (d=0.05), the change of the error threshold is slight. As a whole, the averaged relative concentrations in this case are basically consistent with those obtained from the deterministic Crow-Kimura model. Nevertheless, with the increment of the fluctuation strength, the change of the error threshold becomes more and more obvious. The error threshold extends downwardly as well as upwardly. The downward extension suppresses the dominance of the wild class in the population and therefore destroys the quasispecies structure. The upward extension pushes the real error threshold to a larger value. This situation is illustrated (Figure 3) in which d=0.25, the error threshold has completely turned into a smooth crossover region. Although the error threshold is modified significantly for a sizable fluctuation strength, the relative concentrations in the region outside the crossover region are almost consistent with those given by the deterministic Crow-Kimura model, which implies that they are relatively stable against the site mutation rate fluctuation.

Fig. 2 The ensemble averaged relative concentrations of the classes in their steady states versus the site mutation rate (μ) with d=0.05 and the number of random samples n=10 000

NOTE: Number 0 represents the master class,and numbers 1,2,3… denote different mutant classes.

Fig. 3 Same as Figure 2 but d=0.25

The error threshold can be changed by the fluctuation of the site mutation rate and could be even sizable when the fluctuation strength is large. This can be seen from an alternative angle. Let us consider an order parameter m=1-2h/N , where h is the averaged Hamming distance of the population from the wild class weighted by the relative concentrations of the classes in their steady states^[29]. Since the averaged Hamming distance carries the information of the relative concentrations the above order parameter may indicate the change due to the randomization of the site mutation rate. For m=1, the population is composed entirely of the wild class. When m=0, all classes of the population are distributed randomly in the sequence space. In the case of m=-1, the population consists completely of the complementary class to the wild class. The order parameter m as a function of the site mutation rate μ is presented (Figure 4) for the cases of d=0, 0.05 and 0.25. It can be seen that the order parameter in the deterministic Crow-Kimura model declines dramatically at the error threshold, which demonstrates that quantity m really acts as an order parameter to indicate the phase transition. When d is small, the deviation of the order parameter from that obtained from the deterministic Crow-Kimura model is small even around the error threshold, as shown in the case of d=0.05 (Figure 4). The deviation becomes more and more significant, especially around the error threshold as d gets larger and larger, and the phase transition at the error threshold is replaced by a smooth crossover around the error threshold. We show the above situation in the case of d=0.25 where the smooth curve of the order parameter passes through the error threshold and extends to a region far away from the error threshold. Therefore, both the relative concentrations and the order parameter indicate the error threshold is no longer a phase transition point but a smooth crossover region in the presence of a sizable fluctuation of the site mutation rate.

Fig. 4 The order parameter in the population steady state versus the site mutation rate μ for different values of the fluctuation strength （d） and the number of random samples n=10 000

NOTE: The solid line is for d=0(for the deterministic Crow-Kimura model). The dashed line is for d=0.05,and the dashed-dotted line is for d=0.25.

2.2　The characteristics and comparison of the crossover regions appearing in the randomized Crow- Kimura and Eigen models

The position of the error threshold in the deterministic Crow-Kimura model is easily determined, since the distribution of the classes has a dramatic change at the error threshold. However, in the randomized Crow-Kimura model, the change is smooth and the phase transition point is replaced by the smooth crossover region. Precisely determining the error threshold position is then impossible and unnecessary. Instead, we use the width of the crossover region to describe the quantitative change of the error threshold due to the randomization, which was firstly introduced by Li et al. ^[26,27]. The width is defined as follows. Its starting point is the phase transition point (the error threshold) in the deterministic Crow-Kimura model, and its endpoint is the location where the relative difference of the relative concentration of two complementary classes is less than 0.01. The relative difference is given by

with i + j = N. For different complementary classes, their relative concentrations behave in a similar manner in the crossover region. With the above definition, the width of the crossover region w can be evaluated. When the fluctuation strength d is taken to be 0.05, 0.10, 0.15, 0.20, 0.25 and 0.30 the corresponding values of the width w are 0.04, 0.115, 0.215, 0.310, 0.440 and 0.598, respectively. It is shown that as the fluctuation strength increases the width gets wider and wider. It can be found that the relationship between the width and the fluctuation strength of the site mutation rate is nonlinear (exponential) by fitting the above values (Figure 5) in which the width of the crossover region versus the fluctuation strength in the cases of the random fitness and random site mutation rate is presented. Feng et al. ^[25] demonstrated in that there is a crossover region around the error threshold as the fitness of the master and mutant classes are replaced by Gaussian distributed random variables in the Crow-Kimura model. In order to fully examine the characteristics of the error threshold in the randomized Crow-Kimura model (Figure 5), our mutation randomization results are compared with those given by Feng et al. By fitting the width values, we found that the relationship between the width and the fluctuation strength is linear in the randomized fitness case.

Fig. 5 The width of the crossover region varies with the fluctuation strength in the two different randomizations of the Crow-Kimura model.

NOTE: The solid line with circle denotes the width in the random site mutation rate case. The line with triangle represents that in the random fitness case which is obtained by fitting the width values^[25].

For the randomized Crow-Kimura model, in the case of the randomized site mutation rate, when the fluctuation strength d is taken to be 0.15 and 0.25, the width of the crossover region wμ is 44% and 90% of the value of the error threshold, respectively. In the case of the randomized fitness, when the fluctuation strength d is equal to 0.15 and 0.25 the corresponding width wf is 32% and 54% of the error threshold. Therefore, the width caused by the randomized fitness is comparable to that caused by the randomized site mutation rate. Nevertheless, for the randomized Eigen model, when d equals 0.15 and 0.25 the width wμ is 47% and 104% in the case of the randomized site mutation rate. In the case of the randomized fitness, when d is taken to be 0.15 and 0.25 the width wf is only 9% and 16%^[26]. So, the width caused by the randomized fitness is much less than that caused by the randomized site mutation rate, namely, the extension of the error threshold is mainly caused by the randomized site mutation rate for the same fluctuation strength. The above results are shown in Table 1. The relationship between the width and the fluctuation strength in the randomized Eigen model is linear for the randomized fitness and nonlinear (actually exponential) for the randomized site mutation rate, which is the same as that in the randomized Crow-Kimura model. Meanwhile, the fitness fluctuation in the randomized Crow-Kimura model causes much more extension of the error threshold than that in the randomized Eigen model and the site mutation rate fluctuation in the randomized Eigen model induces more extension of the error threshold than that in the randomized Crow-Kimura model. Therefore, the two randomized models are the same qualitatively but quite different quantitatively in the sense of the relationship between the extension and the fluctuation strength.

3　Concluding remarks

In the present paper, in order to complement the randomization of the Eigen model and the Crow-Kimura model the site mutation rate in the Crow-Kimura model has been treated as a Gaussian distributed random variable. We paid attention to the characteristics of the error threshold as well as the relationship between the extension of the error threshold and the fluctuation strength of the randomized mutation rate. It has been shown that both the relative concentrations and the order parameter indicate the error threshold is no longer a phase transition point but a smooth crossover region in the presence of a sizable fluctuation of the site mutation rate. The error threshold extends not only downwardly but also upwardly. The downward extension suppresses the dominance of the wild class in the population and destroys the quasispecies structure thereof. The upward extension pushes the real error threshold to a larger value. The relationship between the width and the fluctuation strength in the mutation randomized Crow-Kimura model is nonlinear. The obtained results were compared with those from the randomized Eigen model and it is found that for the two randomized models the relationship between the width and the fluctuation strength is linear for the randomized fitness and nonlinear (exponential) for the randomized site mutation rate. Therefore, they are the same qualitatively in the sense of the relationship between the extension and the fluctuation strength. Nevertheless, they are quite different quantitatively. For the randomized Crow-Kimura model the width caused by the randomized fitness is comparable to that caused by the randomized site mutation rate. For the randomized Eigen model the extension of the error threshold is mainly caused by the randomized site mutation rate. Up to now we have had a full picture about the randomization effects of the fitness and site mutation rate on the error threshold based on the Eigen model and the Crow Kimura model.

The full picture might have the following implications for anti-viral strategies, cancer therapy and breeding of animals and plants. (1) The downward and upward extension of the error threshold suggests that one may use the fluctuation of the fitness and site mutation rate to destroy the quasispecies structure and roughly kill viruses. The extension itself implies a kind of residual life which appears in the crossover region^[27]. Therefore, the upper limit of the crossover region should be reached to completely drive biological populations to go extinct^[30]. (2) Radiation mutagens or chemical mutagens have been used to increase the mutation rate of viruses to eliminate them, which means one may also change and control the width of the crossover region externally, for instance, by electromagnetic field, radiations and thermal fluctuations for many purposes.

References