text:
# Four-state quantum chain as a model of sequence evolution

## 1 Introduction

One prominent phenomenon in the theory of molecular evolution that has also attracted considerable attention in statistical physics is the so-called error threshold. It describes the breakdown of genetic order in mutation-selection models for mutation rates surpassing a certain critical value. The prototype model for the description of the error threshold is Eigen’s quasispecies model in sequence space  (which is effectively equivalent to a coupled mutation-selection model in population genetics, cf ), originally designed for the description of prebiotic RNA evolution. However, the threshold is supposed to be a phenomenon that should occur in a rather general class of mutation-selection models.

In order to set up a mutation-selection model that is tractable by analytical (or at least numerical) methods, severe simplifications of the original biological situation seem to be indispensable. Analytical approaches generally have to restrict to the treatment of infinitly large populations and rather simple fitness functions, such as the sharply peaked landscape of Eigens original model. Another common approximation, also used in previous studies of the quasispecies model, amounts for the simplified representation of genotypes as binary strings. In the context of molecular evolutionary theory, this may be thought of as representing DNA or RNA strands by sequences of purins and pyrimidins, hence with only two states per site, neglecting the fact that genetic information is really given by a four-letter alphabet. In this article, we present a four-state mutation-selection model which is capable to describe the full nucleotide alphabet and incorporates the standard mutation schemes of molecular phylogeny. In particular, the phase diagramms are discussed in detail which are more polymorphic than for the two-state model. This shows that, for a full understanding of the error threshold behaviour in molecular evolution, investigations can not be restricted entirely to the study of two-state models.

One important step towards an understanding of the threshold phenomenon has been its identification with an equilibrium phase transition in physics by the translation of a time-discrete version of the quasispecies model into the transfer matrix of an anisotropic two-dimensional Ising model . This equivalence was further exploited to study various aspects of the error threshold with methods from statistical physics . It turns out, however, that the anisotropy of that model is not so easy to handle and the analysis of the relevant biological quantities (which correspond to certain surface properties of the Ising model) remains an involved problem. Due to the complications of the model, almost all results obtained so far are approximate or numerical. The only exact result for the sharply peaked landscape  has been worked out via a different analogy to a model of directed polymers, using the specific properties of that very special fitness landscape.

An alternative approach to the analysis of mutation-selection models and the error threshold which avoids some of the problems of the anisotropic Ising model has been brought up in . Here, the starting point on the biological side is a slightly changed model which describes the evolution of a population with overlapping generations in continuous time. It turns out that, after a reformulation in tensor products, the two-state version of this model is equivalent to the Hamiltonian of an Ising quantum chain. Thereby, the change to continuous time in the biological description corresponds to the anisotropic limit that connects the two-dimensional Ising model and the quantum chain in physics (cf. ). The quantum chain model is technically easier to handle, and exact results for two non-trivial fitness landscapes, namely Onsager’s landscape and the quadratic fitness function, have been worked out .

Accordingly, we extend this latter approach to a full four-state model in this study. The quantum chain analogy allows to use well-known methods from statistical mechanics for the solution of the model, so that we do not have to dwell on technical details here. For an extended presentation of methods (with regard to the two-state model) using techniques from rigorous mean field theory, we refer to . The main focus is instead on the discussion of the threshold behaviour and in particular the increased complexity of the phase diagram due to the consideration of the four-state nature of biological information and the refined schemes of molecular mutation rates.

In the following section, we start with a presentation of the biological foundations of our model. Only thereafter, we will introduce the quantum chain model in Section 3. In Section 4, analytical and numerical results are presented for a number of specific four-state models with permutation invariant fitness landscapes. Also the properties of finite sequences and the evolution dynamics will be studied. We close with a summary of our results and a discussion of open problems in Section 5.

## 2 Biological foundations

Genetic information is coded in DNA (and RNA) molecules. These are heteropolymers of four units (nucleotides) which differ in a specific base. The essential aspect of a DNA sequence is captured in a string over a four-letter alphabet

$$𝝈VV_1\times V_2\times \mathrm{}\times V_N;V_i=\{A,C,G,T\}$$
(1)
where each letter represents a particular base: $`A`$ and $`G`$ for adenine and guanine (the purins), $`C`$ and $`T`$ for cytosine and thymine (the pyrimidins). In RNA sequences, $`T`$ is replaced by $`U`$ for uracil. We will therefore treat the $`4^N`$ different sequences of a fixed, finite length $`N`$ as our genotypes (which may be thought of as coding for something, such as a virus or an enzyme). Disregarding environmental effects, we may identify a collection of genotypes with a population of haploid ‘individuals’. Evolution then describes the change of the population composition in time.

A standard model for the evolution of an infinite, asexually reproducing population under the basic forces of mutation and selection which works in continuous time is given by the following system of non-linear differential equations

$$\dot{p}_𝝈(t)=\left(r_𝝈\overline{r}(t)\right)p_𝝈(t)+\underset{𝝈^{}}{}m_{𝝈𝝈^{}}p_𝝈^{}(t).$$
(2)
Here, $`p_𝝈(t)`$ denotes the relative frequency of genotype $`𝝈`$ at time $`t`$ with corresponding Malthusian fitness (replication rate minus death rate) $`r_𝝈`$, and

$$\overline{r}(t)=\underset{𝝈}{}r_𝝈p_𝝈(t)$$
(3)
is the mean fitness of the population. It is the origin of the non-linearity in (2). Finally, $`m_{𝝈𝝈^{}}`$ is the (time independent) rate at which $`𝝈^{}`$ mutates to $`𝝈`$. This framework has originally been defined in classical population genetics . In the sequence space context, it has been introduced in  and has been called the para-muse (parallel mutation-selection) model, since it assumes mutation and selection to act independently and in parallel at each instant of time. The model ignores recombination and genetic drift due to finite population size. Both assumptions can be considered as fairly reasonable at least in the context of the evolution of viruses or bacteria where populations can be huge and recombination is absent, or the nucleotides are tightly linked. In the following subsections, the basic processes of mutation and selection shall be described in some detail.

### 2.1 Mutation

We take mutation as a point process acting independently on all sites, ignoring more complicated mechanisms, such as insertions or deletions. Molecular mutation rates shall be chosen according to the following scheme, known as the Kimura 3 ST model in molecular phylogeny :

Within this general setup, a number of simpler models is contained, which treat mutation at different levels of sophistication. In the simplest approach, the mutation rates between all four nucleotides are assumed to be equal $`(\mu _1=\mu _2=\mu _3)`$. This is the so-called Jukes-Cantor mutation scheme. While this simple frame already seems to be sufficient for a number of applications, measurements reveal that there are indeed pronounced differences in the mutation rates that should be accounted for in more realistic models. In particular, the transitions between the two purins (A,G) and the two pyrimidins (C,T) are much more frequent than the purin–pyrimidin mutations which are called transversions. This may range up to relative differences of $`\mu _1\mu _3\mu _2/2`$ in the nucleus and $`\mu _1\mu _3\mu _2/40`$ in mitochondrial DNA . A mutation scheme with $`\mu _2>\mu _1=\mu _3`$ is known as the Kimura 2 parameter model. The full Kimura 3 ST scheme, finally, also accounts for the small difference between $`\mu _1`$ and $`\mu _3`$, such that $`\mu _2>\mu _1>\mu _3`$.

Implementing this mutation model into the evolution equation (2), we obtain the following mutation rates between genotypes ($`i\{1,2,3\}`$)

$$m_{𝝈𝝈^{}}=\{\begin{array}{cc}\hfill \mu _i,& d_i(𝝈,𝝈^{})=d_{𝝈𝝈^{}}=1\hfill \\ \hfill N_i\mu _i,& 𝝈=𝝈^{}\hfill \\ \hfill 0,& d_{𝝈𝝈^{}}>1\hfill \end{array}.$$
(4)
Here,

$`d_1(𝝈,𝝈^{})`$ $`=`$ $`\mathrm{\#}_{AC}(𝝈,𝝈^{})+\mathrm{\#}_{GT}(𝝈,𝝈^{})`$
$`d_2(𝝈,𝝈^{})`$ $`=`$ $`\mathrm{\#}_{AG}(𝝈,𝝈^{})+\mathrm{\#}_{CT}(𝝈,𝝈^{})`$ (5)
$`d_3(𝝈,𝝈^{})`$ $`=`$ $`\mathrm{\#}_{AT}(𝝈,𝝈^{})+\mathrm{\#}_{CG}(𝝈,𝝈^{})`$
are restricted Hamming distances between $`𝝈`$ and $`𝝈^{}`$. In (5), $`\mathrm{\#}_{XY}(𝝈,𝝈^{})`$ counts the positions at which $`X`$ and $`Y`$ are exchanged in $`𝝈`$ and $`𝝈^{}`$. Finally,

$$d_{𝝈𝝈^{}}=d_1(𝝈,𝝈^{})+d_2(𝝈,𝝈^{})+d_3(𝝈,𝝈^{})$$
(6)
is the total Hamming distance. Note that the choice of the diagonal term $`m_{𝝈𝝈}`$ in (4) just accounts for probability conservation ($`_𝝈\dot{p}_𝝈=0`$) in the mutation part of the evolution equation (2).

### 2.2 Selection and fitness landscape

Whereas the mutational part of the dynamics is fairly well understood at least on the microscopic (molecular) level, the relation of genotype and fitness, which defines the respective selective success, is notoriously complex. Following the standard notion in molecular evolution, we define the fitness function (or fitness landscape)

$$f:𝝈r_𝝈$$
(7)
as a mapping from the configuration space $`V=\{A,C,G,T\}^N`$ into the real numbers, assigning a reproduction rate (Malthusian fitness value) $`r_𝝈`$ to each genotype. Implicitly, the fitness function incorporates all the complicated interactions between the sites. These interactions are typically long-ranged (since RNA strands or proteins fold in three dimensions), highly correlated, and give rise to rather rugged landscapes. Especially in the context of RNA evolution, the construction and characterization of fitness landscapes has motivated numerous studies, see e.g.  for a review.

Below we will show how the evolution equation (2), with an arbitrary choice of the fitness function, can be adapted to the methods from statistical physics by a reformulation in a quantum chain framework. As an application, we then present a study (including analytical and numerical results) for specific examples from the class of permutation invariant fitness functions. Here, due to equivalence of all sites, the fitness of a given genotype is solely a function of its restricted Hamming distances from the so called wildtype sequence with optimal fitness which we choose as the reference genotype. This particularly simple class of fitness landscapes is widely used, as a canonical first approximation, especially in multilocus theory. Also in the context of sequence space evolution, fitness functions of this type have been used in a number of studies on the two-state model . To implement the approach in our four-state model, we fix an arbitrary sequence, denoted by $`𝝈_{++}`$, as the wildtype. We will only consider directional selection here towards a unique genotype with optimal fitness. The fitness of any other sequence is then determined by the restricted Hamming distances $`d_i`$ relative to $`𝝈_{++}`$. Permutation invariance with respect to the position in the sequence thus leads to a drastic reduction of dimensions. For the four-state model, the effective configuration space forms a tetrahedron in 3d (see Fig. 2) and is conveniently represented in Cartesian coordinates which we shall call (following ) the surplus components:

$`s_1(𝝈)`$ $`=`$ $`1{\displaystyle \frac{2}{N}}\left(d_1(𝝈,𝝈_{++})+d_3(𝝈,𝝈_{++})\right);`$
$`s_2(𝝈)`$ $`=`$ $`1{\displaystyle \frac{2}{N}}\left(d_2(𝝈,𝝈_{++})+d_3(𝝈,𝝈_{++})\right);`$ (8)
$`s_3(𝝈)`$ $`=`$ $`1{\displaystyle \frac{2}{N}}\left(d_1(𝝈,𝝈_{++})+d_2(𝝈,𝝈_{++})\right).`$
With this choice, any unstructured random sequence has coordinates $`s_i0`$ (with probability 1 in the limit $`N\mathrm{}`$). Any positive value of a surplus component, on the other hand, signals a non-trivial overlap of the sequence with the wildtype $`𝝈_{++}`$. In particular, $`s_1`$ measures the surplus of sites with purins or pyrimidins as given in $`𝝈_{++}`$ over the purin–pyrimidin mutated sites.

Within this frame, a natural class of permutation invariant fitness functions is

$$f:𝝈r_𝝈=N\underset{i=1}{\overset{3}{}}\left[\alpha _is_i(𝝈)+\frac{\gamma _i}{2}s_i^2(𝝈)\right]$$
(9)
which includes the following special cases

* Setting $`\alpha _i>0`$ and $`\gamma _i=0`$, we obtain the purely additive Fujiyama landscape without genetic interactions. Here, every mutation relative to the wildtype has a fixed deleterious effect, independent of any other mutation that may be present in the sequence. The additive landscape is a canonical zeroth-order approximation, ignoring any kind of genetic interactions. In the context of sequence evolution, this fitness function has been discussed e.g. in .
* With the choice $`\alpha _i\gamma _i>0`$, the model corresponds to a concave quadratic fitness function (with directional selection) as it is frequently met in multilocus theory. Due to the gene interactions, existing mutations tend to aggravate further ones, which is called positive epistasis.
* For $`\alpha _i0`$ and $`\gamma _i>0`$, we finally obtain a convex fitness function for directional selection with long-range gene interactions and negative epistasis (existing mutations tend to alleviate further ones). Since we want to have $`𝝈_{++}`$ as unique wildtype sequence and a fitness function which is monotonous in the surplus components, we restrict $`f`$ to the octant $`s_i0`$ and (smoothly) truncate the fitness function by introduction of a step function $`\mathrm{\Theta }(s_i)`$ whenever frequencies of genotypes with $`s_i<0`$ are non-zero:

$$\stackrel{~}{f}:𝝈r_𝝈=N\underset{i=1}{\overset{3}{}}\left[\left(\alpha _is_i(𝝈)+\frac{\gamma _i}{2}s_i^2(𝝈)\right)\mathrm{\Theta }(s_i)\right].$$
(10)

The variables $`\alpha _i`$ and $`\gamma _i`$ may further be used to distinguish between the effects of the different types of mutations (as defined in Fig 1) on the fitness. In this article, we will present explicit results for the two following cases:

1. For the simplest choice, $`\alpha _1=\alpha _2=\alpha _3`$ and $`\gamma _1=\gamma _2=\gamma _3`$, any mutation away from the wildtype has the same effect. Together with the Jukes-Cantor mutation scheme, symmetry here leads to equal values of the surplus components in the mutation–selection equilibrium. The model may thus also be thought of as a two-state model, where any site is only regarded as occupied either with a wildtype or with a mutant nucleotide. In contrast to the simple two-state model of , however, there is an effectively asymmetric mutation rate between wildtype and mutant in the case considered here.
2. In a more refined model, we distinguish between transitions and transversions. In the mutational part, this is done by applying the Kimura 2 parameter mutation scheme. In the fitness function, we take into account that the deleterious effects of the transversions often dominate over those of the transitions: $`\alpha _1>\alpha _{2,3}`$ and/or $`\gamma _1>\gamma _{2,3}`$.

## 3 Quantum chain model

### 3.1 Symmetries

Since mutation is a random process that is independent of the fitness values of the genotypes involved, the molecular mutation scheme consequently makes no reference to fitness concepts like the wildtype. Biological observables measurable from sequence data, such as the surplus components (8), and also the fitness functions as defined in (9) or (10), on the other hand, are defined relative to the wildtype sequence. In order to set up these concepts in a common framework, it is convenient to reformulate also the mutational part of the evolution equation in coordinates relative to the wildtype. This may always be done due to certain symmetries inherent in the mutation scheme of Fig. 1.

The basic symmetry of the mutation scheme, if all three mutation rates $`\mu _1,\mu _2,\mu _3`$ are pairwise different, is $`C_2\times C_2`$ (Klein’s 4-group), generated by two involutions. If we write the operations in standard permutation notation, we can take as generators the transformations

$$\left(\begin{array}{cccc}A& C& G& T\\ C& A& T& G\end{array}\right)\text{and}\left(\begin{array}{cccc}A& C& G& T\\ G& T& A& C\end{array}\right),$$
(11)
both being the product of two transpositions. This symmetry may now be exploited for a redefinition of the mutation scheme in wildtype coordinates. To this end, we fix, for every site of the wildtype sequence, the element of the 4-group (in the above representation) with the letter of the wildtype nucleotide in the first position (e.g. the string $`(T,G,C,A)`$ for wildtype nuceotide $`T`$). An alternative representation of the configuration space in wildtype coordinates as

$$𝝈V^\pm V_1^\pm \times V_2^\pm \times \mathrm{}\times V_N^\pm ;V_i^\pm =\{++,+,+,\}$$
(12)
is now given by the mapping, on each site, of the string of labels $`(++,+,+,)`$ to the symmetry element of 4-group defined above. With this notation, the three types of mutations included in the Kimura 3 ST scheme simply switch the signs of the labels: $`\pm \pm \pm `$ at rate $`\mu _1`$, $`\pm \pm \pm `$ at rate $`\mu _2`$, and $`\pm \pm `$ at rate $`\mu _3`$.

Higher symmetries of the mutation model are obtained if mutation rates are equal. For the Kimura 2 parameter scheme, $`\mu _1=\mu _3\mu _2`$, the operation

$$ACGTA=\left(\begin{array}{cccc}A& C& G& T\\ C& G& T& A\end{array}\right)$$
(13)
is also a symmetry and generates a cyclic group $`C_4`$. Together with the previous $`C_2\times C_2`$, this generates a dihedral group, $`D_4`$, with 8 elements. Finally, if $`\mu _1=\mu _2=\mu _3`$, we additionally get the simple transposition $`AC`$ and have the full permutation group $`S_4`$ as symmetry. Note that $`S_4`$, which corresponds to the full tetrahedral group with 24 elements, is also the symmetry group of the configuration space of permutation invariant configurations visualized in Fig. 2. The global symmetry (with the same transformation acting at each site simultaneously) of our class of mutation-selection models with fitness functions according to (9) is therefore always a subgroup of $`S_4`$. In particular, the symmetric fitness model with $`\alpha _1=\alpha _2=\alpha _3`$, $`\gamma _1=\gamma _2=\gamma _3`$, and Jukes-Cantor mutation scheme possesses $`C_{3v}`$ symmetry, or the full tetrahedral symmetry if the linear part in the fitness function vanishes ($`\alpha _i=0`$). The transition-transversion model finally, with $`\alpha _1>\alpha _2=\alpha _3`$, or $`\gamma _1>\gamma _2=\gamma _3`$, and Kimura 2 parameter mutation has simple $`C_2`$ symmetry, or $`D_4`$ symmetry if $`\alpha _i0`$. In the latter case, the combination of $`\gamma _2=\gamma _3`$ with $`\mu _1=\mu _3`$ is necessary, not a misprint. Other combinations with global $`D_4`$ symmetry are $`(\gamma _1=\gamma _3;\mu _2=\mu _3)`$ and $`(\gamma _1=\gamma _2;\mu _1=\mu _2)`$.

### 3.2 Construction

With the above preparations, we may now follow the lines of  where the two-state model is treated.

In a first step, we represent the $`4^N`$-dimensional vector space in which we describe the genotype frequencies as the $`N`$-fold tensor product space $`W=_{j=1}^NW_j`$. Hereby, the configuration space $`V^\pm `$ is canonically embedded in $`W`$ by the mapping of the elements of $`V_i^\pm `$ onto the basis vectors $`\{e_j^{++},e_j^+,e_j^+,e_j^{}\}`$ of $`W_j^4`$. Since the nonlinear part in the differential equations (2) only amounts to normalization of the frequencies, a transformation to so-called absolute frequencies

$$z_𝝈(t)=p_𝝈(t)\mathrm{exp}\left(\underset{𝝈^{}}{}r_𝝈^{}_0^tp_𝝈^{}(\tau )𝑑\tau \right)$$
(14)
then reduces the system to the linear equation

$$\dot{z}_𝝈(t)=\left(+\right)z_𝝈(t)$$
(15)
where the mutation and reproduction matrices, $`=(m_{𝝈𝝈^{}})`$ and $`=\text{diag}(r_𝝈)`$, may now be conveniently represented in the frequency space $`W`$. Defining

$$\sigma _j^{(\alpha ,\beta )}:=\left(^{j1}\text{ 1}\text{1}_4\right)(\sigma ^\alpha \sigma ^\beta )\left(^{Nj1}\text{ 1}\text{1}_4\right)$$
(16)
where $`\sigma ^\alpha `$, $`\alpha \{0,x,z\}`$, are the real Pauli matrices and $`\sigma ^0\text{ 1}\text{1}_2`$, we find

$$=\underset{j=1}{\overset{N}{}}\left[\mu _1\sigma _j^{(x,0)}+\mu _2\sigma _j^{(0,x)}+\mu _3\sigma _j^{(x,x)}(\mu _1+\mu _2+\mu _3)\text{ 1}\text{1}\right]$$
(17)
for the mutation matrix. The reproduction matrix $``$ is, for a general fitness landscape, an element of the algebra generated by $`\sigma _j^{(z,0)}`$ and $`\sigma _j^{(0,z)}`$, $`1jN`$,

$$=r_0\text{ 1}\text{1}+\underset{k,\mathrm{}=1}{\overset{N}{}}\underset{[j_1\mathrm{}j_k]}{}\underset{[j_1\mathrm{}j_{\mathrm{}}]}{}\epsilon _{[j_1\mathrm{}j_k],[j_1\mathrm{}j_{\mathrm{}}]}\underset{m=1}{\overset{k}{}}\sigma _{j_m}^{(z,0)}\underset{n=1}{\overset{\mathrm{}}{}}\sigma _{j_n}^{(0,z)},$$
(18)
where $`[j_1\mathrm{}j_k]`$ is an ordered $`k`$-tupel in $`\{1,\mathrm{},N\}`$. Now, from a physical point of view, $`=+`$ is (up to a global minus sign) the Hamiltonian of two coupled Ising quantum chains in a tunable transverse magnetic field (the mutation) and general spin-interactions within the chains.

Translated to our quantum chain model, the fitness function of the permutation invariant landscape defined in (9) results in a (longitudinal) magnetic field and a mean field spin-interaction. We find $`=_\alpha +_\gamma `$, where

$$_\alpha =\underset{j=1}{\overset{N}{}}\left[\alpha _1\sigma _j^{(z,0)}+\alpha _2\sigma _j^{(0,z)}+\alpha _3\sigma _j^{(z,z)}\right]$$
(19)
and

$$_\gamma =\frac{1}{2N}\underset{j,k=1}{\overset{N}{}}\left[\gamma _1\sigma _j^{(z,0)}\sigma _k^{(z,0)}+\gamma _2\sigma _j^{(0,z)}\sigma _k^{(0,z)}+\gamma _3\sigma _j^{(z,z)}\sigma _k^{(z,z)}\right]$$
(20)
Let us stress that, in contrast to most physical applications, the mean field model is a much more natural approach in the biological context where interactions are typically long-range. So, it is a legitimate model here, not an inevitable approximation.

### 3.3 Biological and physical observables

In this subsection, we relate the quantities of biological interest, mean and variance of the surplus components and the fitness, to the physical observables. In what follows, we assume the occuring limits to exist.

#### Genotype composition

According to (15), the Hamiltonian of the quantum chain determines the time evolution of our population of genotypes in an environment that does not constrain the population size. For any genotype-independent regulation of the population size, the relative genotype frequencies are found by statistical normalization. We therefore define the vector of the genotype composition $`|𝒑(t)`$ and the equilibrium composition $`|0`$ as

$$|𝒑(t)=\frac{\mathrm{exp}(t)|𝒑_0}{\mathrm{\Omega }|\mathrm{exp}(t)|𝒑_0};|0:=\underset{t\mathrm{}}{lim}|𝒑(t)$$
(21)
where $`|𝒑_0`$ is the initial composition and $`4^N|\mathrm{\Omega }`$ is the equidistribution of genotypes. Note that the equilibrium composition of the genotype population just corresponds to the ground state of the quantum chain on the physical side (with a different ‘biological’ normalization $`\mathrm{\Omega }|0=1`$).

#### Fitness

The density of the mean fitness (or mean fitness per site) of the population is given by the expression

$$w(t):=N^1\overline{r}(t)=N^1\mathrm{\Omega }||𝒑(t).$$
(22)
Since

$$w:=\underset{t\mathrm{}}{lim}w(t)=N^1\mathrm{\Omega }||0=N^1\frac{0||0}{0|0}$$
(23)
the equilibrium mean fitness (per site) is just given by the (unique) largest eigenvalue of $``$, corresponding to $`|0`$. For an unconstrained population, $`w`$ also determines the growth rate in the long-time limit. In the physical picture, $`(w)`$ is obviously just the ground state energy (per spin).

Using $`|\mathrm{\Omega }=0`$, we derive for the time evolution of the mean fitness

$$\dot{w}(t)=V_r(t)+N^1\mathrm{\Omega }|[,]|𝒑(t)$$
(24)
where $`V_r(t)`$ is the variance of fitness (per site),

$$V_r(t)=\frac{1}{N}\left(\mathrm{\Omega }|^2|𝒑(t)\mathrm{\Omega }||𝒑(t)^2\right).$$
(25)
In the absence of mutation, (24) is of course just a special case of Fisher’s “Fundamental Theorem of Natural Selection”  which states that the rate of increase in fitness is equal to the genetic variance in fitness. For the mutation-selection models considered here, the relation has the following intuitive interpretation: The change in mean fitness is driven by two independent forces. The first one stems from the change of genotype frequencies due to selection and is proportional to the variance of fitness values present in the population. Since variances are positive, it always tends to increase fitness. The second term on the right hand side of (24) typically decreases fitness. It measures the population mean of the change in fitness at time $`t`$ due to the action of mutation. In mutation-selection equilibrium, both terms balance, and the entire residual variance is due to mutation.

#### Surplus

Another quantity that characterizes the genetic order of the population, as it may be measured from sequence data, is the mean surplus. We define, following and generalizing ,

$$u_i(t)=\underset{𝝈}{}s_i(𝝈)p_𝝈(t);u_i=\underset{t\mathrm{}}{lim}u_i(t).$$
(26)
In particular,

$$\mathrm{\#}_m(t):=\frac{1}{4}\left(3(u_1(t)+u_2(t)+u_3(t))\right)$$
(27)
measures the mean number of mutations per site relative to the wildtype while

$$\mathrm{\#}_{tr}(t):=\frac{1}{2}\left(1u_1(t)\right)$$
(28)
denotes the mean number of transversions alone. As a biological order parameter, the mean surplus plays a similar rôle as the physical magnetization. However, as already noted in , both quantities are quite distinct and in many cases not even easily related. In the language of the quantum chain, the equilibrium mean surplus may be derived as

$$u_1=\frac{\mathrm{\Omega }|_i\sigma _i^{(z,0)}|0}{N};u_2=\frac{\mathrm{\Omega }|_i\sigma _i^{(0,z)}|0}{N};u_3=\frac{\mathrm{\Omega }|_i\sigma _i^{(z,z)}|0}{N},$$
(29)
whereas the three-component magnetization is defined as the ground state expectation value

$$m_1=\frac{0|_i\sigma _i^{(z,0)}|0}{N0|0};m_2=\frac{0|_i\sigma _i^{(0,z)}|0}{N0|0};m_3=\frac{0|_i\sigma _i^{(z,z)}|0}{N0|0}.$$
(30)
As we will show below, magnetization and surplus can show rather different behaviour especially near phase transitions. The biological and physical phase diagrams, however, coincide if phase transitions (or error thresholds) are defined as nonanalyticity points of the ground state energy (or mean fitness) $`w`$ in the thermodynamic limit (cf. the discussion in Section 5).

## 4 Results

### 4.1 Fujiyama model

As in the two-letter case , the quantum chain model decomposes into non-interacting one-site Hamiltonians for the additive landscape. The mean fitness and its variance are linear functions in the surplus components. In particular, we obtain from (24)

$$V_r(t)=\dot{w}(t)+2\left((\mu _1+\mu _3)\alpha _1u_1(t)+(\mu _2+\mu _3)\alpha _2u_2(t)+(\mu _1+\mu _2)\alpha _3u_3(t)\right).$$
(31)
For Jukes-Cantor mutation, $`\mu _1=\mu _2=\mu _3\mu `$, this reduces to

$$V_r(t)=\left(4\mu +\frac{\text{d}}{\text{d}t}\right)w(t)$$
(32)
and $`V_r`$ is proportional to the mean fitness in the mutation–selection equilibrium. Exact results are easily found from the solution of the four-dimensional eigenvalue problem of the one-site Hamiltonian. We only give the expression for the mean fitness in the symmetric case, $`\alpha _1=\alpha _2=\alpha _3\alpha `$ with Jukes-Cantor mutation scheme ($`\mu _1=\mu _2=\mu _3\mu `$):

$$w(t)=\frac{\mathrm{exp}[2t(\alpha +\mu )]\mathrm{cosh}[2tQ]\left(\alpha 2\mu +2Q\mathrm{tanh}[2tQ]\right)\alpha 4\mu }{1+\mathrm{exp}[2t(\alpha +\mu )]\mathrm{cosh}[2tQ]}$$
(33)
where

$$Q=\sqrt{\mu ^2+\alpha ^2\alpha \mu }$$
(34)
and the equidistribution of genotypes is chosen as starting configuration.

Means and variances of the fitness and the surplus in mutation–selection balance are shown in Fig. 6 below. A plot of the time evolution of fitness is given in Fig. 8. There is clearly no phase transition (resp. no error threshold behaviour) for the additive Fujiyama landscape, as expected in view of the complete absence of interactions (resp. epistasis).

### 4.2 Quadratic fitness model: Equilibrium results

In contrast to the additive case, no simple relation between surplus and fitness is known in the case of the quadratic landscape as long as $`t`$ or $`N`$ are kept finite. However, due to the permutation invariance of the Hamiltonian, the individual fitness–surplus relation (9) is recovered in the thermodynamic limit for the corresponding mean values of the equilibrium population. We obtain in analogy to :

$$w=\underset{t\mathrm{}}{lim}w(t)=\underset{i=1}{\overset{3}{}}\left(\alpha _iu_i+\frac{\gamma _i}{2}u_i^2\right)$$
(35)
and, from (24), for the equilibrium variance of fitness per site

$$\begin{array}{c}V_r=\underset{t\mathrm{}}{lim}V_r(t)=2(\mu _1+\mu _3)\left(\alpha _1u_1+\gamma _1u_1^2\right)+\hfill \\ \hfill 2(\mu _2+\mu _3)\left(\alpha _2u_2+\gamma _2u_2^2\right)+2(\mu _1+\mu _2)\left(\alpha _3u_3+\gamma _3u_3^2\right).\end{array}$$
(36)
The key to the solution in the thermodynamic limit is now the minimum principle of the physical free energy which translates to a maximum principle for the equilibrium mean fitness. Maximizing

$$𝒙|+|𝒙w\left(𝒙|𝒙1\right)$$
(37)
with respect to $`w`$ and $`𝒙`$, we obtain, taking permutation symmetry of $`𝒙`$ into account, the following variational expression for $`w`$:

$`w(𝜶,𝝁`$ $`,𝜸)=sup_{m_1,m_2,m_3}[\alpha _1m_1+\alpha _2m_2+\alpha _3m_3+{\displaystyle \frac{\gamma _1}{2}}m_1^2+{\displaystyle \frac{\gamma _2}{2}}m_2^2+{\displaystyle \frac{\gamma _3}{2}}m_3^2+`$
$`{\displaystyle \frac{\mu _1}{2}}\left(\sqrt{(1+m_2)^2(m_1+m_3)^2}+\sqrt{(1m_2)^2(m_1m_3)^2}2\right)+`$
$`{\displaystyle \frac{\mu _2}{2}}\left(\sqrt{(1+m_1)^2(m_2+m_3)^2}+\sqrt{(1m_1)^2(m_2m_3)^2}2\right)+`$
$`{\displaystyle \frac{\mu _3}{2}}(\sqrt{(1+m_3)^2(m_1+m_2)^2}+\sqrt{(1m_3)^2(m_1m_2)^2}2)]`$
where $`m_i[1,1]`$ are the components of the physical magnetization. Let us stress that, from the biological point of view, the translation to the physical framework seems a necessary technical step since we do not know of any variational principle for the biological model which works directly in $`L^1`$. We now take a closer look at two special cases.

#### Symmetric fitness model

For the symmetric wildtype–mutant model with $`\alpha _i\alpha `$, $`\gamma _i\gamma `$ and Jukes-Cantor mutation rate $`\mu `$, all components of the order parameters are equal, $`m_im`$ and $`u_iu`$, respectively. Here, the variational expression (4.2) for $`w`$ leads to the following self-consistency condition for $`m`$:

$$m=\frac{1}{3}\left[1+\frac{2(\alpha +\gamma m)\mu }{\sqrt{(\alpha +\gamma m)^2\mu (\alpha +\gamma m)+\mu ^2}}\right].$$
(38)
This is a quartic equation in $`m`$ and can be solved using the standard formulas. However, since the explicit solution is rather lengthly, we do not include it here, but give a qualitative discussion instead.

Obviously, the relation has a unique real solution for any $`\alpha `$ and $`\mu `$ whenever $`\gamma `$ is negative. Like in the case of the two-state model, we thus obtain no phase transition for positive epistasis. In the following, we therefore concentrate our discussion on positive $`\gamma `$ (or negative epistasis). Note that, for calculations in the thermodynamic limit, always the fitness function $`f`$ (9), and hence the reproduction matrix $`_\gamma `$ (20), can be used instead of the truncated form $`\stackrel{~}{f}`$ (10), since the frequencies of genotypes with negative surplus vanish. For $`\alpha _i0`$, this is due to spontaneous breaking of the extra $`C_2\times C_2`$ symmetry of $`=+_\gamma `$.

In contrast to the two-state model, where a phase transition in the thermodynamic limit is only found for zero external field, it turns out that the present model has phase transitions for a whole range of the linear fitness parameter $`\alpha `$ when epistasis is negative: For $`\stackrel{~}{\alpha }:=\alpha /\gamma `$ in the interval

$$0\stackrel{~}{\alpha }<\frac{1}{3}\left(\sqrt{\frac{4}{3}}1\right)0.0515668$$
(39)
we find a first order phase transition of the system at

$$\stackrel{~}{\mu }:=\frac{\mu }{\gamma }=\stackrel{~}{\mu }_c=\frac{2}{3}+2\stackrel{~}{\alpha }$$
(40)
with a finite jump in the magnetization from $`m_+`$ to $`m_{}`$ where

$$m_\pm =\frac{1}{3}\left(1\pm \sqrt{127\stackrel{~}{\alpha }^218\stackrel{~}{\alpha }}\right).$$
(41)
From $`m`$ we derive the mean fitness $`w`$ using (4.2), from $`w`$ we obtain the surplus $`u`$ via (35) and, finally, the variance of the fitness $`V_r=12\mu (\alpha u+\gamma u^2)`$. Looking at the surplus $`u`$, we also find a phase transition at $`\stackrel{~}{\mu }=\stackrel{~}{\mu }_c`$. As $`m`$, it vanishes in the disordered phase for $`\alpha =0`$. Note however that, since $`w`$ is continuous, due to the relation (35), also the surplus is continuous at a phase transition. In  it has been shown that these differences of the biological and physical order parameters arise with the change from classical to quantum mechanical probabilities (resp. the change from $`L^1`$ to $`L^2`$) in translating the biological model into the physical one. We remark that a different, discontinuous behaviour of the biological order parameter at a (physical) first order transition has been observed for the sharply peaked landscape in Eigen’s quasispecies model . Mean fitness and its variance, magnetization, and surplus for different values of $`\alpha `$ are shown below in Fig. 3.

#### Transition–transversion model

In our second example, we wish to distinguish mutations between like and unlike nucleotides. In a first step, we retain the symmetric fitness landscape $`\gamma _1=\gamma _2=\gamma _3\gamma `$ (for simplicity with vanishing linear part $`\alpha =0`$), but let the relative frequencies of transitions and transversions differ by assuming the Kimura 2 parameter mutation scheme, $`\mu _1=\mu _3\mu \mu _2`$.

In the extended parameter space of the reduced mutation rates $`\stackrel{~}{\mu }=\mu /\gamma `$; $`\stackrel{~}{\mu }_2=\mu _2/\gamma `$, we now obtain a phase diagram with three distinct phases (see Fig. 4).

* For $`\stackrel{~}{\mu }`$ and $`\stackrel{~}{\mu }_2`$ sufficiently small, all three surplus components are positive, indicating genetic order with respect to the entire 4-letter alphabet of the nucleotides: ACGT phase.
* If we increase the mutation rate $`\stackrel{~}{\mu }_2`$ for low $`\stackrel{~}{\mu }`$, the system crosses over to a phase which does no longer distinguish between the different kinds of purins (A,G) and pyrimidins (C,T), but is still ordered with respect to transversions. This is the limiting case described by the two-state model. We call this the PP phase.
* For higher mutation rates $`\stackrel{~}{\mu },\stackrel{~}{\mu }_2`$, we finally enter a completely disordered phase with vanishing fitness and surplus.

In a second step, we now also let the mutation effects of transitions and transversions differ and assume a fitness landscape with $`\gamma _2=\gamma _3\gamma `$, but $`\gamma _1\gamma `$ in general. The changes in the phase diagram for increasing $`\stackrel{~}{\gamma }_1=\gamma _1/\gamma `$ are shown in Fig. 5. The phase transitions between the three phases may be first or second order. In general, we obtain the following phase space structure:

* Phase transitions between the disordered and PP phase are second order and located on the line $`\stackrel{~}{\mu }=\stackrel{~}{\gamma }_1/2`$. This phase transition corresponds to the one also seen in the two-state model .
* The phase transition line between the ACGT and PP phases in general changes from first to second order with increasing mutation rate $`\stackrel{~}{\mu }_2`$ (see Figs. 4, 5). For the second order transitions we derive, on expanding (4.2) to lowest order in $`m_2=m_3`$,

$$\mu =\frac{\gamma _1}{\gamma _1+2\gamma }\sqrt{(\gamma _1+\mu _2)(2\gamma \mu _2)}.$$
(42)
Numerically, we find that the first order transitions are located on a straight line up to $`\stackrel{~}{\mu }=\stackrel{~}{\gamma }_1/2`$ where the PP phase changes into the disordered phase. The $`\stackrel{~}{\mu }_2`$-interval of first-order transitions decreases for increasing $`\stackrel{~}{\gamma }_1`$. For $`\stackrel{~}{\gamma }_18.45`$, all phase transitions between the ACGT and PP phases are second order.
* Finally, for $`\stackrel{~}{\gamma }_14`$, there are direct first order phase transitions between the ACGT phase and the disordered phase (for $`\stackrel{~}{\mu }_2`$ sufficiently small). For higher values of $`\stackrel{~}{\gamma }_1`$, these two phases are separated by the PP phase.

As for the symmetric fitness function discussed above, there are no compact analytic expressions for the fitness or the surplus in the ACGT phase. In the PP phase, however, the following values for the mean fitness and the non-zero components of the mean surplus and the magnetization are found:

$$w=\frac{\gamma _1}{2}(1\frac{2\mu }{\gamma _1})^2;u_1=1\frac{2\mu }{\gamma _1};m_1=\sqrt{1\left(\frac{2\mu }{\gamma _1}\right)^2}.$$
(43)
The variance in fitness per site, finally, is proportional to the mean fitness in the PP phase: $`V_r=8\mu w`$. Note that all these expressions are independent of the transition rate $`\mu _2`$ and directly comparable to the results of the two-state model  by idebtifying $`\{++,+\}`$ with ‘$`+`$’ and $`\{+,\}`$ with ‘$``$’.

### 4.3 Quadratic fitness model: Finite sequence length

For the Fujiyama model with independent sites, all the quantities calculated here, means and variances per site in infinite populations, are independent of the assumed length $`N`$ of the sequences. This is no longer the case for models including epistasis. In this subsection, we therefore present a quick numerical investigation of the symmetric fitness model for finite system sizes and compare the results with those in the thermodynamic limit. Since the frequencies of genotypes with negative values of the surplus no longer vanish for finite sequences, we use the truncated fitness function (10), with $`\gamma _i\gamma >0`$ and $`\alpha _i=0`$ for our calculations.

All results are obtained by direct numerical solution of the eigenvalue problem in the $`[(N+1)(N+2)(N+3)/6]`$-dimensional vector space of permutation invariant population vectors. Numerically precise calculations have been performed up to $`N=60`$ (39711-dim.), the results are shown in Fig. 6. It is seen that the mean surplus and the mean and the variance of the fitness rapidly approach the limiting curves and behave qualitatively different from the Fujiyama model even for very small system sizes. We also show the finite-size behaviour of the variance of the surplus $`V_s`$. Since this quantity vanishes as $`1/N`$, it is not obtainable from the leading order terms in the thermodynamic limit. In our finite size calculations, we rescale $`V_s`$ with the sequence length to obtain comparable results. Whereas $`V_s`$ is monotonously increasing for the additive model (where $`NV_s=1u^2`$), it runs through a maximum for quadratic fitness. Note that this maximum, in contrast to the variance of fitness, is located directly at the error threshold. The behaviour is qualitatively similar to the two-state model .

Since there has been some discussion recently on the correct scaling of fitness values and mutation rates with the length of the sequence (cf ), let us finally remark that the finite size results in this and the next section show that our choice, keeping fitness and mutation rate per site fixed, is adequate for all quantities considered here.

### 4.4 Quadratic fitness model: Time evolution

Originally, the error threshold has been defined as an equilibrium phenomenon (cf ): For special forms of the fitness landscape, there is a finite critical value $`\mu _c`$ of the mutation rate beyond which genetic order is no longer maintained by selection. For the four-state model with quadratic fitness, this situation has been discussed above. However, for a suitable fitness function, the threshold is not necessarily connected with high mutation rates. In this subsection, we consider the relaxation of a non-equilibrium population to mutation-selection balance. It turns out that, depending on the starting configuration, an even stronger threshold effect may be observed in the time evolution of the fitness and the surplus for all mutation rates below the critical equilibrium value.

#### Zero-mutation limit of the transition-transversion model

The essence of the threshold phenomenon in the time evolution is already contained in the selection dynamics alone. In a first step, we therefore disregard mutation altogether by working in the zero-mutation limit. Obviously, we then deal with a classical mean-field model on the physical side. As our starting configuration, we choose the completely unstructured population with an equidistribution of genotypes $`|𝒑_0=4^N|\mathrm{\Omega }`$. In this particular situation, some progress is possible also analytically. Noting that

$$\widehat{C}(t)=\frac{\mathrm{\Omega }|\widehat{C}\mathrm{exp}(t)|\mathrm{\Omega }}{\mathrm{\Omega }|\mathrm{exp}(t)|\mathrm{\Omega }}=\frac{\text{tr}(\widehat{C}\mathrm{exp}(t))}{\text{tr}(\mathrm{exp}(t))}$$
(44)
for any element $`\widehat{C}`$ of the algebra generated by $`\{\sigma _i^{(z,0)},\sigma _i^{(0,z)}\}`$, the biological and physical pictures coincide in this case. Using the fitness function of the transition-transversion model with $`\gamma _2=\gamma _3\gamma >0`$, we obtain the following implicit equations for the time evolution of the surplus components:

$`u`$ $`=`$ $`{\displaystyle \frac{\mathrm{sinh}(2\gamma tu)}{\mathrm{cosh}(2\gamma tu)+\mathrm{exp}[2\gamma _1t(2u\mathrm{coth}(2\gamma tu)1)]}}`$ (45)
$`u_1`$ $`=`$ $`{\displaystyle \frac{\mathrm{cosh}[\gamma tQ(u_1)]\mathrm{exp}(2\gamma _1tu_1)}{\mathrm{cosh}[\gamma tQ(u_1)]+\mathrm{exp}(2\gamma _1tu_1)}}`$ (46)
where

$$Q(u_1)=\sqrt{(1+u_1)^2\mathrm{exp}(4\gamma _1tu_1)(1u_1)^2}.$$
(47)
The resulting dynamical phase diagram is shown in Fig. 7. As in the equilibrium situation, there are three phases. Depending on the ratio $`\stackrel{~}{\gamma }_1=\gamma _1/\gamma `$, the system directly crosses to an ordered phase after a sharply defined waiting time $`t_c`$, or performs two consecutive transitions, entering the PP phase in the first one.

As in the equilibrium phase diagram, the dynamical transitions may be of first or second order.

* Second order transitions are located at $`\stackrel{~}{t}=\gamma t=1`$ for $`\stackrel{~}{\gamma }1/4`$ and at $`\stackrel{~}{t}=1/\stackrel{~}{\gamma }_1`$ for the transition from the disordered phase to the PP phase. The transition from the PP phase to the ACGT phase is second order above $`\stackrel{~}{\gamma }_11.9009`$ and implicitly given through $`2\stackrel{~}{t}_c=1+\mathrm{exp}[2\stackrel{~}{\gamma }_1(\stackrel{~}{t}_c1)]`$. A similar second order transition (with a one-component order parameter) has also been observed in the two-state model .
* In an interval around the symmetry point $`\gamma _1=\gamma `$, the system possesses a first order transition (in the sense that there is a finite jump in the magnetization). Note that, in contrast to the equilibrium case, also the surplus and even the mean fitness are discontinous on this line, giving rise to a rather pronounced threshold effect in the evolution dynamics (cf. the solid line in Fig. 8 for $`\stackrel{~}{\gamma }=1`$).

As for the equilibrium values, we also consider the effect of finite sequence lengths on the time evolution. Again, calculations are performed by direct diagonalization of the symmetric fitness model ($`\stackrel{~}{\gamma }=1`$). Fig. 8 shows how the jump discontinouity in the mean fitness (internal energy) and the delta function singularity in the variance of the fitness (specific heat) are approached by the finite systems. A threshold phenomenon is absent in the time evolution of the Fujiyama model which is also shown in Fig. 8.

#### Finite mutation rates and different starting configurations

In a last step, we now discuss the influence of the mutation rate and the starting configuration on the evolution dynamics. Consider first the time evolution of the equilibrium distribution of genotypes $`4^N|\mathrm{\Omega }`$. Although no analytical results are available here, we may give the following intuitive argument that there is a phase transition at finite $`t=t_c`$ for any mutation rate below the critical equilibrium mutation rate $`\mu _c`$: Since mutation alone tries to keep the population in the equilibrium distribution, the evolution dynamics will be slowed down by mutation for small $`t`$. In particular, mean fitness and surplus will remain zero on a finite interval at least up to the threshold value of the corresponding model with vanishing mutation. On the other hand, the limiting values of $`w`$ and $`u`$ are finite for $`\mu <\mu _c`$, giving rise to a non-analytical point of $`w(t)`$ and $`u(t)`$ at some finite $`t=t_c`$. As shown in the upper graph of Fig. 9, this behaviour is clearly visible in numerical results for finite sequence sizes.

In order to contrast the time evolution of the unstructured population with an equidistribution of genotypes as starting configuration, we have also performed calculations for the opposite case of a population with initially homogeneous phenotypes. Here, at $`t=0`$, any ”individual” in the population has the same value $`s_i=0`$ for the three surplus components. The result (for finite sequence length $`N=60`$) is shown in the lower viewgraph of Fig. 9. As for the equidistribution of genotypes, there is a clear threshold effect in the time evolution for any finite value $`0<\mu <\mu _c`$ of the mutation rate. The transition appears to be particularly sharp for small mutation rates. In contrast to the unstructured case, the critical waiting time $`t_c`$ for the transition is no longer monotonously increasing with the mutation rate $`\mu `$, but is separated in two regimes: For mutation rates near the equilibrium threshold value $`\mu _c`$, the situation is similar to the unstructured case: Here, single mutants with higher fitness appear in the population after a short while. Due to the continuing mutation pressure, however, a certain time is needed for these fitter individuals to grow to a finite proportion and to dominate the mean values in the infinite population. For small $`\mu `$, on the other hand, the critical waiting time $`t_c`$ is dominated by the time needed for mutation to explore the configuration space and to generate individuals with higher fitness at a sufficient rate.

## 5 Discussion

When in  a class of models for sequence space evolution was introduced, using the framework of Ising quantum chains, the calculations started with four major simplifications of the biological situation. These are the consideration of a two-state model, the assumption of an infinite sequence length, the use of simplistic fitness landscapes, and the restriction on infinite population sizes. In this paper, we have looked at the first two of these simplifying assumptions. Finally, an extended discussion of the evolution dynamics of these models has also been presented. In the following paragraphs, we give a summary of our findings and an outlook on the remaining open problems.

#### Two-state versus four-state models.

The main concern of this contribution is the generalization of the modelling framework, introduced in , to four states (corresponding to the four nucleotides) on each site. The generalization presented makes use of the $`C_2\times C_2`$ symmetry inherent in the Kimura 3 ST mutation scheme. On the ‘physical side’ this leads to a model of two coupled Ising quantum chains (rather than to a four-state Potts model). Compared with the two-state model, the extension can be thought of as consisting of two steps. In a first step, we represent the four states on each site by the spin values of two spins in decoupled chains. Note that already in this simplified model three phases occur in the phase diagram since the transition lines of the two decoupled chains will not in general coincide. The second step consists of the introduction of a more realistic mutation scheme which also changes the configuration space topology and the corresponding use of a refined fitness landscape. Both these extensions lead to a coupling of the chains, and an even richer phase space structure is found, including first-order transitions. As may be seen from the introduction of a small linear field term into the fitness function in subsection 4.2, this change of the transition to first order leads to an increased robustness of the threshold phenomena with respect to symmetry-breaking perturbations.

#### Finite sequence length.

Typical sequence lengths of enzymes or viruses are of the order $`10^3`$$`10^4`$. While these numbers are certainly far off the typical sizes of macroscopic systems in physics, they are, in principle, large enough to successfully supress $`1/N`$-corrections. However, especially models with simple fitness landscapes describe – at best – the evolution dynamics in a very restricted configuration space of particularly ‘important’ sites, disregarding neutral or altogether lethal mutations. In view of this fact, consideration of finite sequence lengths is indispensible and calculations in the thermodynamic limit even seem to be questionable at first sight. In order to clarify the usefulness of infinite-size methods in this context, we performed a number of numerical calculations for finite sequence lengths. The results are quite encouraging. As shown in subsection 4.3, the characteristic properties of the thermodynamic limit are well visible even for tiny sequence sizes, such as $`N=10`$, and the approximation is already quantitatively reasonable for sequences of length $`60`$.

#### The fitness landscape.

The construction of a tractable fitness landscape which nevertheless comprises the relevant biology is certainly the major task for all these models. In this contribution, in order to obtain at least some analytical results, we have chosen a fitness function from the smooth end of the landscape zoo. Due to its permutation invariance, the quadratic fitness function effectively disregards any local variance in the interaction between sites, but only considers the average epistatic effect. As such, it is in many respects certainly no more than a toy-model for evolution. However, the assumption of permutation invariance of the sites is quite common in evolutionary biology and comprises a large number of standard models for evolution, such as the quadratic optimum model or Eigen’s original sharply peaked landscape. The results show that the essential structure responsible for characteristic effects such as the error threshold is already contained in this simplified framework and may also serve as a reference for future work on fitness functions with increased ruggedness, such as the NK-landscape hierarchy . Here, we expect the results for the quadratic fitness model to be qualitatively stable at least under certain forms of mild ruggedness, such as the introduction of site-randomness in the fields and interactions . Pronounced changes, on the other hand, should be expected when spin-glass effects come into play.

#### Finite population size.

In going from the deterministic limit to the evolution of finite populations, the ordinary differential equation (2) has to be replaced by the master equation of a stochastic process which is no longer covered by the theoretical framework presented in this article. Due to the complexity of the stochastic equations, analytical results seem to be out of reach at present for all but the simplest selection schemes. Monte-Carlo simulations, however, should be possible and could considerably add to theoretical insight here.

Although the general picture of the deterministic case should persist at least for sufficiently large populations, the study of finite population effects is certainly of importance. For related models, such as the quasispecies model with the single peaked landscape, it is has been found  that the deterministic results can be interpreted as the time averages of the stochastic process for mutation rates outside a certain interval around an error transition. Directly at the threshold, however, large fluctuations and a jump in the long-time averages appear in the stochastic system at a critical mutation rate which seems to be lower by an amount roughly proportional to $`1/\sqrt{N}`$ in comparison with the deterministic case. Mainly because of these expected finite population effects we have restricted discussions in this article entirely to the phase space structure of the models and the order of the phase transitions. Any further details of the transitions, even critical exponents, will presumably never be visible in real biological systems and thus seem to be of limited relevance in this context.

Let us finally remark that, although biological populations are certainly finite, the consideration of the infinite population limit is not (only) a technical necessity, but also of direct importance for the study of the error threshold. That is so because this effect, in distinction to the phenomenon of Muller’s ratchet, is by definition not due to genetic drift, but solely due to the form of the fitness function. It has thus always to be shown that the threshold effect persists even for infinitly large population sizes.

#### Error threshold behaviour.

Since there are more than one and sometimes conflicting definitions of the error threshold in literature (cf. the discussion in ), let us start this paragraph with a few clarifying remarks. In this article, following , we use the notion of the error threshold as equivalent to phase transitions. As such, a clear-cut mathematical definition (as non-analytical points in the mean fitness) is possible only in the infinite sites (or thermodynamic) limit. However, since the thermodynamic limit can be considered as an excellent approximation already for rather small systems, the infinite system property gives a valid explanation for prominent features which are observable for finite sequences as well. In our study, we have always considered sequences of a fixed length and have treated the mutation rate per site as the variable driving the transition. In comparing systems of different length, we have scaled the variables such that a well-defined limit is approached as $`N\mathrm{}`$. In particular, the ‘critical’ mutation rate per site in a finite system quickly converges to the limiting value $`\stackrel{~}{\mu }_c`$. Originally, the threshold has been viewed as a limitating factor on the sequence length . This, however, should not be confusing: We switch to this latter picture simply by letting the reduced mutation rate depend linearly on the sequence length, $`\stackrel{~}{\mu }N`$, and obtain a critical length $`N_c\stackrel{~}{\mu }_c`$ (for sufficiently large sequences).

Our results on the error threshold phenomenon fit previous ones for the two-state case and related models in that negative epistasis is needed to observe a transition (cf. ). Contrary to the two-state case, the threshold corresponds to a first-order transition for certain parameter ranges and persists for a sufficiently small linear part in the fitness function. Both, the equilibrium and the dynamical phase diagram of the transition-transversion model (with $`\alpha _i=0`$), possess two ordered phases characterized by non-zero values of one or all three components of the surplus order-parameter and the disordered phase with zero surplus where selection ceases to operate. The threshold effect appears to be especially sharp in the evolution dynamics, where a jump in the mean surplus and fitness and a delta singularity in the variance of fitness occurs.

Besides the threshold effect, however, other properties of mutation-selection models may be studied within the framework presented. After all, exclusive concentration on phase transitions is perhaps too much a physicist’s point of view on these systems. The relations between surplus, mutation rate and the variance of fitness (24), (36), for example, are valid for the entire time evolution and arbitrary mutation rates. Depending on the fitness function applied, they may give rise to characteristic features also far off the transition point. This is particularly explicit for the equilibrium variance of fitness which runs through a pronounced maximum for fitness functions with negative epistasis at a mutation rate much smaller than the threshold value.

## Acknowledgments

It is our pleasure to thank Ellen Baake and Oliver Redner for numerous discussions and comments on the manuscript. Financial support from the German Science Foundation (DFG) is gratefully acknowledged.