Dynamic Fitness Landscapes in Molecular Evolution
We study self-replicating molecules under externally varying conditions. Changing conditions such as temperature variations and/or alterations in the environment’s resource composition lead to both non-constant replication and decay rates of the molecules. In general, therefore, molecular evolution takes place in a dynamic rather than a static fitness landscape. We incorporate dynamic replication and decay rates into the standard quasispecies theory of molecular evolution, and show that for periodic time-dependencies, a system of evolving molecules enters a limit cycle for . For fast periodic changes, we show that molecules adapt to the time-averaged fitness landscape, whereas for slow changes they track the variations in the landscape arbitrarily closely. We derive a general approximation method that allows us to calculate the attractor of time-periodic landscapes, and demonstrate using several examples that the results of the approximation and the limiting cases of very slow and very fast changes are in perfect agreement. We also discuss landscapes with arbitrary time dependencies, and show that very fast changes again lead to a system that adapts to the time-averaged landscape. Finally, we analyze the dynamics of a finite population of molecules in a dynamic landscape, and discuss its relation to the infinite population limit.
The main problem one encounters when dealing with dynamic landscapes is the difficulty to find a good generalization of the quasispecies concept. In the original work of Eigen, the quasispecies is the equilibrium distribution of the different molecular species. It is reached if the system is left undisturbed for a sufficiently long time. Since in a dynamic landscape the system is being disturbed by the landscape itself, the concept of a quasispecies is meaningless in the general case. However, there are special cases in which a meaningful quasispecies can be defined. If, for example, the landscape changes on a much slower time scale than what the system needs to reach the equilibrium, then the system is virtually in equilibrium all the time, and the concentrations at time are determined from the landscape present at that time. Generally, it is certain symmetries in the dynamics of the landscape that allow for the definition of a quasispecies. One example we treat in this paper in detail is the case of time-periodic landscapes, which offer a natural quasispecies definition.
An early investigation of dynamic landscapes has been carried out by Jones [4, 5], who has only considered cases in which all replication rates change by a common factor. Therefore, this approach excludes (among other cases) in particular all situations in which the order of the molecules’ replication rates changes over time, i.e., in which for example one of the faster replicating molecules becomes one of the slower replicating molecules and vice versa. Recent work on dynamic fitness landscapes allow for such changes. Wilke et al. [26, 17] have developed a framework that allows to define and to calculate numerically a quasispecies in time-periodic landscapes. Independently, Nilsson and Snoad  have studied the particular example of a stochastically jumping peak in an otherwise flat landscape. This work has been generalized by Ronnewinkel et al. , who could also define a meaningful quasispecies for a deterministic version of the jumping peak landscape and related landscapes. Very recent work of Nilsson and Snoad deals with self-adaptation of the mutation rate in a jumping peak landscape , and with the low-pass filter effect of dynamic fitness landscapes , a notion put forward by Hirst [27, 28]. Finally, theoretical studies of dynamic fitness landscapes can also be found in the related field of genetic algorithms. Schmitt et al. [29, 30] derive results for finite populations in a relatively broad class of dynamic landscapes. However, only landscapes in which the fitnesses get scaled can be treated, so that the same restriction applies here that applied to Jones’s work. The order of the fitnesses can never change. Also, genetic algorithms with time-periodic landscapes have been studied by Rowe [31, 32]. However, Rowe’s approach has the disadvantage that it is tightly connected to the discrete time used in genetic algorithms, and that the dimension of the transition matrices grows in proportion to the period length of the oscillation. This makes it hard to derive analytical results, and in addition to that, it renders landscapes with large inaccessible to numerical calculations.
In this report, we do not cover the large field of molecular evolution in the particular landscapes induced by RNA sequence-to-structure mappings [33, 34, 35], and, in connection to that, we do not discuss neutral networks [36, 37, 38, 39]. We do so mainly because these topics have so far not been studied in the light of temporal variations in the fitness landscapes, and hence, a discussion of them does not fit very well into the general tenor of this work. In general, it can be argued that neutral networks are of less importance in a dynamic environment, because in that situation a population is constantly on the move to the next local optimum .
The remainder of this report is structured as follows. We begin our discussion in Section 2.2 with a brief summary of the general aspects of dynamic fitness landscapes in the quasispecies model. In Section 3, we will develop the main subject of this work: a general theory of time-periodic fitness landscapes. The theoretical part thereof is presented in Section 3.1, in which we demonstrate how a time-dependent quasispecies can be defined by means of the monodromy matrix, and how this monodromy matrix can be expanded in terms of the oscillation period . In Section 3.2, we present an alternative approximation formula for the monodromy matrix that is more suitable for numerical calculations, and in Section 3.3, we compare, for several example landscapes, the results obtained from that formula with the general theory developed in Section 3.1. The restriction of a time-periodic fitness landscape is weakened in Section 4, where we discuss the implications of our findings for other, non-periodic fitness landscapes. Since our work is based on Eigen’s deterministic approach with differential equations, all results presented up to the end of Section 4 are only valid for infinite population sizes. In order to address this shortcoming, in Section 5 we give a brief introduction into the problems involved when dealing with finite populations. In Section 5.1, some simulation results are shown, demonstrating the relationship between the results from the infinite population limit and the actual finite population dynamics. Finally, an approximative analytical description of a finite population evolving on a simple periodic landscape is developed in Section 5.2. We close this paper with some conclusions in Section 6.
2 The quasispecies model
2.1 Static landscapes
With his model of self-replicating molecules, Eigen showed for the first time that Darwin’s idea of mutation and selection can work in a simple, seemingly “lifeless” system of chemical reactants. His observation that evolution is governed by the laws of physics spawned a fruitful field of work, and hundreds of papers based on his initial ideas have been written since then. Most of that work has been concerned with static fitness landscapes. There exist comprehensive reviews of that work (see [11, 13] for a very detailed coverage of the literature till 1989, and  for a more recent review). Here, we are going to briefly introduce the main concepts developed for static landscapes. In doing so, we restrict ourselves to those concepts that we will refer to later on in our discussion of dynamic fitness landscapes. For a more in-depth discussion of static landscapes, the reader is referred to the above mentioned reviews.
The quasispecies model was originally aimed at describing self-replicating DNA or RNA molecules. Therefore, the molecules were conceived of as sequences consisting of a limited number of elementary building blocks. With that picture in mind, we may represent the molecules as sequences of letters. For RNA molecules, e.g., the conventional alphabet consists of the 4 letters G, A, C, U, representing the 4 bases guanine, adenine, cytosine, and uracil, respectively. Today, most researchers are interested in the information-theoretic aspects of the model. Consequently, the most common alphabet in the molecular evolution literature has become the binary alphabet, consisting of the letters 0 and 1. Throughout this work, we will also adopt this choice. With regard to the RNA example, the binary alphabet can be considered as distinguishing only between purines (which are guanine and adenine) and pyrimidines (which are cytosine and uracil).
The molecules the model describes have the ability to self-replicate. Self-replication is a complicated process, which consumes energy and substrates from the environment. These external resources are supposed to be present, and are not modeled explicitly. The degree to which a macromolecule finds the resources necessary to self-replicate, and is able to exploit them, is expressed in the replication coefficients . A molecule that finds good conditions for self-replication has a high , another molecule which is a bad self-replicator has a much smaller .
The molecules replicate by copying themselves. The copy procedure is generally not error-free. Among the different types of errors one can imagine for the copy of a sequence of letters (substitutions, insertions, deletions), we consider only substitutions. Substitutions can in principle occur with a different probability at every single position in the sequence. However, if we assume the copy procedure to be performed step by step by some sort of molecular machinery, the probability of copying a symbol incorrectly can be expected to be the same for all positions in the sequence. Hence, a letter will change from 0 to 1 or vice versa during the copy procedure with a fixed probability, denoted by . This probability is called the error rate, or, alternatively, the mutation rate.
Molecules are also subject to decay with a particular rate for species . A molecule that decays is assumed to be completely absorbed by the environment, i.e., it does not break into parts that are themselves being described by the model.
The constant production of new molecules due to the ongoing self-replication will drastically increase the concentration of molecules, and will decrease the amount of resources available for further self-replication. We are interested in the description of a steady state, and therefore we have to introduce a constant dilution which lets new resources enter the system and washes away the excess production of those molecules. The total concentration of molecules can thus be kept constant by proper adjustment of the dilution flux.
Finally, we assume that the self-replicating molecules are placed in a well-stirred reactor. As a consequence of this assumption, we can neglect any spatial correlations in the model, and concentrate solely on the molecules’ abundances.
In summary, the model is based on the following postulates:
The molecules are represented by binary sequences of length . They form and decompose steadily. The number of copies of sequence present at time is denoted by .
Sequences enter the system solely as the result of a correct or erroneous copy of another sequence already present.
The substrates necessary for the ongoing replication are always present in sufficient quantity. Excess molecules are washed away by a flux .
The decay of sequences is a Poisson process.
These four assumptions form the basis of Eigen’s theory of molecular evolution. A quantitative analysis of these assumptions can be done in terms of differential equations in the molecules’ occupation numbers . In the following paragraphs, we recall the quantitative analysis of the model that has been developed by Eigen and coworkers.
We begin by writing down an expression for the change in the number of copies of sequence . The abundance of sequence increases because a proportion of the molecules of type replicates faithfully, while some of the other molecules of other types produce offspring of type as the outcome of fruitless attempts to self-replicate. Let the matrix express the probability that molecule copies into molecule . The associated increase in the abundance then amounts to .
The decrease of a molecule’s abundance can have two reasons: its decay, expressed by , and its removal due to the flux term, expressed by . Here, is the total number of molecules at time , i.e., . Putting all the different terms together, we end up with the net change in the number of copies of sequence ,
The quantity gives the probability that a molecule replicates faithfully. is sometimes referred to as the copy fidelity.
For the further development of the theory, it is useful to introduce an average excess production , defined by
The conservation law
expressing the fact that every (possibly erroneous) copy of a molecule represents another molecule in the chemistry, allows us to rewrite the average excess production in terms of the total production rate and the dilution flux, viz.
It proves to be useful to define the matrix as
Now we can rewrite Eq. (1) as
Let us introduce concentration variables , defined by
The quantity measures the relative concentration of molecule in the population. The change in is given by
The total number of molecules grows with
Typically, one assumes that the flux is adjusted such that vanishes, as expressed by Assumption 3 on page 3. However, this assumption is not really necessary for the further development of the theory. Since from this point onwards, we will only be concerned with the relative concentrations , we will disregard the flux altogether, and work with Eq. (9) exclusively.
At this point, it is useful to introduce vector notation, by lumping the concentrations together into a single vector . Equation (9) then becomes
where is the identity matrix. The matrix W can be decomposed into
if we write the replication and the decay coefficients in matrix form as well. Both A and D are diagonal matrices of the form
The matrix Q is the mutation matrix introduced above. Note that the average excess production can also be transformed into vector notation. It takes on the form
where is a vector of 1s, i.e., . The scalar product between and a concentration vector [say ] gives exactly the sum over all components of that vector.
Equation (11) is nonlinear, since the term is quadratic in . Nevertheless, a straightforward solution of Eq (11) is possible, because a transformation exists which removes the nonlinearity. The strength of this transformation lies in the easy reconstruction of the concentration variables from the transformed system. Following Thompson and McBride , or Jones et al. , we introduce
The new variables satisfy the linear equation
Note that if all decay constants are equal, i.e., with a single scalar value , then the transformation
leads to the even simpler equation
The concentration vector can again be obtained from Eq. (17).
The mutation matrix Q has so far remained undefined. The choice of Q specifies the relationship between the different molecules in the chemistry. As was stated above, we conceive the molecules as bitstrings of length , and we assume that copy errors are equally likely on all positions on the string. In that case, the mutation matrix can be written down straightforwardly. Suppose that two bitstrings differ in positions. A mutation from one to the other occurs only if exactly these positions are copied erroneously, while all others are copied faithfully. Such a copy error appears with probability . Hence, we can write the mutation matrix Q as
where is the Hamming distance between two sequences. The Hamming distance is the number of bits in which two sequences differ.
The matrix Q is typically very large, since its number of rows and columns grows as with increasing sequence length . This makes Eq. (16) almost intractable, numerically as well as analytically, for all but the shortest sequences. An alternative matrix is often used, in which certain sequences are grouped together, such that the number of concentration variables can be reduced to . The main idea of this grouping, developed by Swetina and Schuster , is to define a single sequence (which should have the highest replication coefficient of all sequences) as master sequence, and to group all other sequences into error classes, according to their Hamming distance from the chosen master sequence. All the sequences with the same Hamming distance from the master form a single error class. This procedure has the advantage of greatly reducing the number of variables, but it also restricts considerably the number of fitness landscapes that can be studied. Since all sequences in an error class have to share the same replication coefficient, a fitness landscape, for example, in which two peaks have a small or moderately large Hamming distance cannot be defined.
The error class matrix has been given by Nowak and Schuster  in a relatively simple form:
Note that in , the indices and are interchanged inadvertently. The error class matrix can be derived from Eq. (20) by deletion of all but one column of every group of columns whose index yields the same , and by the subsequent summation of all rows whose index yields the same .
From that expression, we can read off that the (unnormalized) variables will grow exponentially over time. Mathematically, this growth can be accompanied by either exponentially damped or exponentially amplified oscillations. For all cases of biological interest, however, there can be at most exponentially damped oscillations. First of all, we will almost always deal with a symmetric matrix Q. A symmetric mutation matrix implies that the probability with which a sequence mutates into a sequence is the same as the probability for the inverse process. Rumschitzki has noted that in this case, the spectrum of W is real . Namely, we can transform the non-symmetric matrix W into a symmetric one by means of a similarity transformation,
The spectrum of the transformed matrix is real because of the matrix’s symmetry, and hence the spectrum of the untransformed matrix must be real as well. In that case, oscillations will be absent in Eq. (22). For non-symmetric Q, we can apply the Frobenius-Perron theorem  if the decay rates satisfy
The Frobenius-Perron theorem is applicable to a matrix with all positive entries, and it guarantees a real largest eigenvalue for that matrix. Consequently, we have at most exponentially damped oscillations as long as (24) is obeyed. In addition to that, the Frobenius-Perron theorem states that the eigenvector corresponding to this largest eigenvalue has only strictly positive entries, and hence, that this eigenvector can be interpreted as a vector of chemical concentrations if normalized appropriately.
Let us write down a more explicit solution to (16), under the assumption that the spectrum of W is known (exact spectra of evolution matrices have been derived in [9, 43, 44]). Let be the eigenvalues of W, and let be the associated eigenvectors. Without loss of generality, we order the eigenvalues such that . The solution to Eq. (16) can then be expressed as
where the coefficients have to be chosen such that the initial condition is satisfied. In other words, the give the expansion of the initial condition in terms of the eigenvectors ,
The principal eigenvector of W has been called the quasispecies by Eigen. The reason for this name will become clear shortly.
In terms of the concentration variables , the solution Eq. (25) becomes
where is again a vector with all entries equal to 1.
As we already mentioned, the Perron-Frobenius theorem assures that the largest eigenvalue is non-degenerate, and that all components of the corresponding eigenvector are positive. We factor out the exponential of the largest eigenvalue in the numerator and denominator of Eq. (27) and obtain
As long as , all contributions apart from the one corresponding to the largest eigenvalue disappear in the limit . Hence, the system evolves towards a steady state, characterized by the dominant eigenvector of the matrix W.
The case is somewhat artificial, because it implies that the system has been started in an exact superposition of eigenstates excluding the dominant quasispecies. In that case, the system cannot evolve towards . Instead, it converges towards the eigenvector corresponding to the next-largest eigenvalue. In a real chemical system with more or less random initial conditions and under the presence of noise, the case is of no relevance.
The appearance of a steady state distribution of concentration variables has important implications for the understanding of Darwinian evolution. In general, the eigenvectors are a mixture of several molecular species . As a consequence, a number of species is present simultaneously in the asymptotic distribution. This means that selection combined with random mutation does not remove all but one molecular species, even if there exists one (the master sequence) whose replication coefficient exceeds all others. Instead, the interplay of selection and mutation creates a cloud of mutant species around the master sequence. It is this cloud that is termed the quasispecies. The mutant distribution forms a unit competing with other similar units, represented by the other eigenvectors. Selection acts on these mixtures of sequences, and ultimately singles out the dominant one, i.e., the one corresponding to the largest eigenvalue.
Interestingly, the master sequence does not necessarily occupy a large fraction of the quasispecies. In fact, as the mutation rate increases, it is common that the master dwindles while the one-mutant and the two-mutant sequences form the prevailing part of the quasispecies. Figure 1 shows the asymptotic sequence distribution versus the error rate for a sharply peaked landscape in which all but one of the replication coefficients are identical, while the remaining one exceeds them significantly. The mutants are grouped into error classes. If copy errors are not present, i.e., for , the master sequence is the only species present, and all others have a vanishing concentration. As soon as takes on values slightly above 0, the concentration of the master sequence starts to decrease. At the same time, the concentrations of the mutants begin to grow, first that of the one-mutant sequence, than that of the two-mutant sequence, and so on. At we observe a sharp transition, beyond which the concentrations take on constant values. The mutation rate at which this transition occurs is called the error threshold. Beyond the error threshold, the outcome of a replication event can be considered a random sequence, and therefore, all molecular species take on the same concentration in this regime. The fact that all concentrations are equal beyond the error threshold is somewhat blurred in Fig. 1 because of the use of error classes. In fact, the heights of the different curves merely reflect the sizes of the corresponding error classes. For , the sequences replicate stochastically in any landscape, since faithfully copied and erroneously copied symbols are then equally likely. The outcome of a single copy process is therefore a random sequence in that case. When comes close to the value 1, almost every symbol is copied incorrectly. This implies that the copy is, apart from mutations, the inverse of the original sequence (note that this is a peculiarity of binary strings). As a consequence, order can again be seen for large . In Fig. 1, the inverse error threshold occurs at . Beyond that point, the inverse master sequence dominates the quasispecies.
If an evolving system displays an error threshold, the amount of information that can be acquired by the system is severely limited. In the case of the sharply peaked landscape, the critical mutation rate at which the error threshold occurs decreases as with increasing string length . The quantity in that expression stands for the relative advantage of the master sequence over the other sequences. For large , the critical mutation rate is therefore very close to 0. This implies that for a given mutation rate, the molecules can grow to a certain length, while beyond that length, all information is lost. Incidentally, it follows that the spontaneous formation of complex self-replicating molecules seems to be impossible. Eigen and Schuster attempted to solve this problem with the concept of the hypercycle . However, it must be said that not all landscapes do present an error threshold. In particular, multiplicative landscapes [46, 47, 48] show a smooth crossover from the completely ordered situation at , where the master sequence is the only sequence present, to , where all molecules take on the same concentration. Moreover, the decay of the master sequence’s concentration depends only on the mutation rate and the selective advantage, but not on the length of the strings. Hence, in a multiplicative model, sequences can grow in length without bound, and self-replicating molecules might form spontaneously.
The error threshold can be viewed as the critical point of a phase transition. This fact is well established through a rich body of work, in which the quasispecies model and related models have been mapped onto spin lattices [49, 10, 50], spin chains [51, 48], and more recently also onto the statistical mechanics of directed polymers in a random medium . The order parameter that is generally being used to describe that phase transition is the average overlap of the sequences in the population with the master sequence, given by the expression , where is the Hamming distance between a sequence of type and the master sequence. The overlap takes on the value if we compare the master sequence with itself, the value if we compare it with its inverse, and intermediate values for all other sequences. The order parameter reads
2.2 Time-dependent replication rates
Having developed a good understanding of the concepts of molecular evolution theory in static environments, let us move on to the non-static case. The basic quasispecies equation (11) changes only in so far as the matrix W now becomes time dependent,
In the most general case, the time dependency may stem from the replication rates, from the decay rates, or even from the mutation matrix:
The average excess production is then
As in the static case, we can transform away the nonlinearity in Eq. (30) with the introduction of unnormalized variables
The resulting equation is again a linear differential equation, however, this time with non-constant coefficients:
In the previous subsection, we were able to immediately write down a solution for the linearized quasispecies equation. In the case of a dynamic fitness landscape, however, no such a simple closed form solution exists. Instead, we have to be satisfied with partial solutions, approximations, and limiting cases. In particular, we cannot generally define a quasispecies as the steady state the system approaches for . Therefore, a central question in relation to dynamic landscapes is the appropriate definition of a quasispecies concept.
For now, we start collecting information about Eq. (34). First of all, we note that we can map the quasispecies model onto a linear system with a symmetric matrix if is symmetric for all . This can be seen by introducing
Unfortunately, we cannot write down a solution for Eq. (36) from the knowledge of the eigensystem of if has an arbitrary time dependency. However, the above mapping shows that for symmetric , we may safely assume that the matrix has only real eigenvalues.
There are two limiting cases that we can discuss without specifying a landscape, namely very fast changes in on the one hand, and very slow changes in on the other hand. We begin with the case of very slow changes. For the rest of this work, we will assume that has a real spectrum for all . To be on the safe side, we also assume that (24) is satisfied for all . In that way, the Perron eigenvector of can always be interpreted as a vector of chemical concentrations.
For every time , we can define a relaxation time
where and are the largest and the second largest eigenvalue of , respectively. The time gives an estimate of how long a linear system with matrix needs to settle into equilibrium. Therefore, if the changes in happen on a timescale much longer than , the system is virtually in equilibrium at any given point in time. Hence, for large enough , the quasispecies will be given by the Perron eigenvector of . Strictly speaking, this is only true if there is always some overlap between the largest eigenvector of and the one of , but in all but some very pathological cases we can assume this to be the case.
The situation of fast changes in is somewhat more difficult, because, as we are going to see later on, we have to define a suitable average over in order to make a general statement. Therefore, we postpone this case for a moment. A detailed discussion of fast changes will be given for the particular case of periodic fitness landscapes in the next section, and later on, we will discuss fast changing landscapes in general.
3 Periodic fitness landscapes
3.1 Differential equation formalism
In this section, we study periodic time dependencies in , for which we can prove several general statements.
If the changes in are periodic, i.e., if there exists a such that
then Eq. (34) turns into a system of linear differential equations with periodic coefficients. Several theorems are known for such systems . Most notably, if is the fundamental matrix, such that every solution to Eq. (34) can be written in the form
then we can define the monodromy matrix ,
which simplifies Eq. (40) to
for the decomposition with the phase . In particular, we have
so that for every phase , we have a well defined asymptotic solution, given by the eigenvector associated with the largest eigenvalue of . Hence, periodic fitness landscapes allow the definition of a quasispecies, much in the same way as static fitness landscapes do. However, here the quasispecies is time-dependent. In other words, a system that evolves in a periodic fitness landscape runs for into a limit cycle with period length .
3.1.1 Neumann series for X
With the knowledge of the monodromy matrix X, the calculation of the system’s limit cycle becomes a simple eigenvalue problem. The monodromy matrix, however, can in general not be given in a closed from. Consequently, we have to rely on expansions in various parameters. In this section, we derive a formal expansion in for the monodromy matrix. This formal expansion is similar in spirit to the Neumann series which gives a formal solution to an integral equation, and it is based on the Picard-Lindelöf iteration for differential equations. As the first step, we have to rewrite Eq. (34) in the form of an integral equation, i.e.,
Our goal is to solve this equation for by iteration. Our initial solution is
which we insert into Eq. (44). As a result, we obtain the 1st order approximation
Further iteration yields
and so on. Now we define
and, in general
and obtain the formal solution
For suitably small , the infinite sum on the right-hand side is guaranteed to converge. When we compare this equation for to the definition of the monodromy matrix Eq. (41), we find that [introducing ]
In particular, since is identical to the time-average over , regardless of , we have the high-frequency expansion
Equation (53) reveals that for very high frequency oscillations, the system behaves as if it was subject to a static landscape. That static landscape is given by the dynamic landscape’s average over one oscillation period.
The radius of convergence of the expansion Eq. (52) can be estimated as follows. Since all entries of are positive, the tensor
can be bound by
from which it follows that
The matrix norm induced by the sum norm
is the column-sum norm
With that norm, and using Eq. (57), we can estimate
Hence, the expansion Eq. (52) converges certainly for those that satisfy
Since all entries in W are positive, we have further
where the bar in and indicates that these quantities are averaged over one oscillation period. The second equality holds because of (24) and because of . Without loss of generality, we assume that the maximum is given by . Then, Eq. (61) is satisfied for
It is interesting to compare this expression to the relaxation time of the time-averaged fitness landscape, . To th order, the principal eigenvalue of is given by . The second largest eigenvalue is to the same order given by the second largest diagonal element of , which we assume to be without loss of generality. Hence, the relaxation time is approximately given by
which is in general larger than the radius of convergence of Eq. (52). In particular, if the largest and the second largest eigenvalue of lie close together, the relaxation time may be much larger than the largest oscillation period for which the expansion is feasible. This restricts the usefulness of Eq. (52) to substantially high frequency oscillations in the landscape. The interesting regime in which the changes in the landscape happen on a time scale comparable to the relaxation time of the system can unfortunately not be studied from Eq. (52).
3.1.2 Exact solutions for and
The two extreme cases (no replication errors) and (random offspring sequences) allow for an exact analytic treatment. The second case is identical to the case of static landscapes, and therefore we will mention it only briefly. At the point of stochastic replication , the population dynamics becomes independent of the details of the landscape. As a consequence, temporal changes in the landscape must become less important as approaches . However, this is not very surprising, since in most cases, an error rate close to 0.5 implies that the population has already passed the error threshold, which in turn implies that it does not feel the changes in the landscape any more.
The case of , on the other hand, is more complex than the corresponding case in a static landscape. Since the matrix Q becomes the identity matrix for , Eq. (34) reduces to
The matrices and are diagonal by definition, and hence, a solution to Eq. (65) is given by
and, in particular,
The integral in the second expression is taken over a complete oscillation period, and hence, it is independent of . Thus, we find for arbitrary
With a vanishing error rate, the monodromy matrix becomes the exponential of the time-average over . Since the exponential function only affects the eigenvalues of a matrix, the quasispecies is given by the principal eigenvector of , irrespective of the length of the oscillation period . In other words, in the absence of mutations will the sequence with the highest average value of take over the whole population after a suitable amount of time, provided it existed already in the population at the beginning of the process. By continuity, this property must extend to very small but positive error rates . So, similar to the case of , the temporal changes in the landscape loose their importance when approaches 0.
There is, however, a caveat to the above argument. In case the largest eigenvalue of is degenerate, temporal changes in the landscape may continue to be of importance for . A degeneracy of the largest eigenvalue of is possible, because the Frobenius-Perron theorem applies only to positive error rates. For degenerate quasispecies, the initial condition determines the composition of the asymptotic population. In this context, let us consider the general solution for periodic fitness landscapes, Eq. (3.1). We have
So even if X becomes independent of for , this need not be the case for , because of Eq. (71). If the largest eigenvalue of is degenerate, variations in will remain visible for arbitrarily large times . Hence, we will observe oscillations among the different quasispecies which correspond to the largest eigenvalue. Clearly, this effect is the more pronounced the larger the oscillation period .
3.1.3 Schematic phase diagrams
The results of the previous two subsections allow us to identify the general properties of the quasispecies model with a periodic fitness landscape at the borders of the parameter space. As parameters, we shall only consider error rate and oscillation period , since all other parameters (replication rates, decay rates, details of the matrix Q) do not influence the above results. In Fig. 2, we have summarized our findings. Along the abscissa runs the oscillation period. For very fast oscillations, the evolving population sees only the time-averaged landscape. For very slow oscillations, on the other hand, the population is able to settle into an equilibrium much faster than the changes in the landscape occur. Hence, the population sees a quasistatic landscape. Along the ordinate, we have displayed the error rate. We discarded the region above , in which anti-correlations between parent and offspring sequences are present, as it is redundant. For , all sequences have random offspring, and hence, all sequences replicate equally well. Therefore, for this error rate, the landscape becomes effectively flat. On the other side, for , we have again the time-averaged landscape. However, for large , the fact that we see the average landscape does not mean that the concentration variables are asymptotically constant. Degeneracies in the largest eigenvalue may cause a remaining time dependency due to oscillations between superimposed quasispecies. The exact form of these oscillations is dependent on the initial condition . For small , the oscillations disappear, because the ratio of newly created sequences during one oscillation period and remaining sequences from the previous oscillation period decays with [Eq. (53)].
From the above observations, we can derive generic phase diagrams for periodic fitness landscapes. There are two main possibilities. The fitness landscape may average to a landscape that has a distinct quasispecies, or it may average to a flat landscape. These two cases are illustrated in Fig. 3. Note that the diagrams are meant to illustrate the qualitative form and position of the different phases. In their exact appearance, they may differ substantially from the exact phase diagram of a particular landscape.
If an averaged landscape sports a distinct quasispecies, then for every oscillation period and every phase of the oscillation , we have a unique error threshold . For small , the error threshold converges towards that of the average fitness landscape, , irrespective of the phase . For larger , the error threshold oscillates between and . In the limit of an infinitely large oscillation period, converges towards , which is the largest error threshold of all the (static) landscapes . Similarly, converges towards in that limit, where is accordingly defined as the smallest error threshold of all landscapes . For a fixed oscillation period and a fixed error rate with , we have necessarily for some phases , and for the rest of the oscillation period. As a result, a quasispecies will form whenever , but it will disappear again as soon as . This phenomenon has for the first time been observed in , where the region of the parameter space in which it can be found has been called the temporarily ordered phase. In this phase, whether we observe order or disorder depends on the particular moment in time at which we study the system. Correspondingly, we will call a phase “ordered” only if order can be seen for the whole oscillation period, and we will call a phase “disordered” if during the whole oscillation period no order can be seen. The relationship between the ordered phase, the disordered phase, and the temporarily ordered phase for the first type of landscapes is displayed in Fig. 3a. Compare also the phase diagram of the oscillating Swetina-Schuster landscape in Fig. 4.
In a landscape that averages to a flat one, on the other hand, the disordered phase must extend over the whole range of for sufficiently small , and order can be observed only above a certain . The behavior of the system for small above cannot in general be predicted solely from the knowledge of Fig. 2. For a landscape with a flat average, the eigenvalues of the monodromy matrix are degenerate for . Therefore, in the limit , the Perron eigenvector can, in principle, converge to any superposition of the eigenvectors of . This situation is visualized in Fig. 3b. If the limit corresponds to a non-homogeneous sequence distribution, the ordered phase extends to arbitrarily small [indicated by the solid line in Fig. 3b]. If, on the other hand, the limit would correspond to a homogeneous sequence distribution, we might find a lower error threshold below which the system would again be in disorder [this is indicated by the dashed line in Fig. 3b]. Since for longer oscillation periods, the oscillations in the degenerate quasispecies at become important, order would be observed for much smaller with increasing . Hence, the lower disordered phase would fade out for . A study investigating under which situations the ordered phase extends to arbitrarily small will be presented elsewhere. As a preliminary result, we can state that the limit does in general not lead to a homogeneous sequence distribution.
3.2 Discrete approximation
The differential equation formalism we have used so far allows for an elegant discussion of the system’s general properties. However, if we want to obtain numerical solutions, this formalism is not very helpful, because we do not have a general expression for the fundamental matrix from Eq. (40), nor for the monodromy matrix from Eq. (41). Therefore, for a numerical treatment we need to move over to the discretized quasispecies equation,
In the case of constant W, the quasispecies obtained from this equation is identical to that of Eq. (34), and it is also identical to that of equation
Equation (73) has been studied by Demetrius et al. , and has been employed by Leuthäusser [49, 10] for her mapping of the quasispecies model onto the Ising model. In the general time-dependent case, however, the additional factor and the identity matrix of Eq. (72) are important, and cannot be left out. The analogue of the fundamental matrix for Eq. (72) reads
where indicates that the matrix product has to be evaluated with the proper time ordering . Similarly, the analogue of the monodromy matrix becomes
where we have assumed that is an integral multiple of , and where . The influence of the size of on the quality of the approximation has been investigated in . A more in-depths discussion of the relationship between the continuous and the discrete quasispecies model can also be found in .
3.3 Example landscapes
For the rest of this section, we are going to study several example landscapes, in order to illustrate the implications of our general theory. In all cases considered, we represent the molecules as bitstrings of fixed length . Moreover, we assume that a single bit is copied erroneously with rate , irrespective of the bit’s type and of its position in the string.
3.3.1 One oscillating peak
In previous work on the quasispecies model with periodic fitness landscapes [26, 17], most emphasis has been laid on landscapes with a single oscillating sharp peak. As a generalization of the work of Swetina and Schuster , the master sequence has been given a replication rate , where is the replication rate of all other sequences. The replication rate has been expressed as
with a -periodic function . The generic example for that function is , leading to
The parameter allows a smooth crossover from a static landscape to one with considerable dynamics, and the exponential assures that is always positive.
In [26, 17] it has been found that the behavior at the border regions of the parameter space is indeed as it is depicted in Fig. 2, and that a phase diagram of the form of Fig. 3a correctly describes the relationship of order and disorder in an oscillating Swetina-Schuster landscape. Here, we put less emphasis on the numerical simulations that lead to these conclusions, but instead show that the phase borders in such a phase diagram can, for an oscillating Swetina-Schuster landscape, be calculated approximately.
For static landscapes with a single peak, the assumption of a vanishing mutational backflow into the master sequence allows to derive an approximate expression for the error threshold [45, 11, 13]. A similar formula can be developed to calculate the error threshold as a function of time in a landscape with a single oscillating peak. But before we turn towards the dynamic landscape, we shall rederive the expression for the master sequence’s concentration in a static landscape, based on neglecting mutational backflow. The expression we shall find is slightly more general than the one that was previously given, and it will be of use for the periodic fitness landscape as well.
The 0th component of the quasispecies equation (30) becomes, after neglecting the mutational backflow,
The average excess production can be expressed in terms of and W as
With that expression, the solution of Eq. (78) requires the knowledge of the stationary mutant concentrations , which are usually unknown. To circumvent this problem, we make the somewhat extreme assumption that all mutant concentrations are equal. Although this assumption, which is equivalent to the assumption of equal excess productions in the usual calculation without mutational backflow, will generally not be true, it works fine for Swetina-Schuster type landscapes. With this additional assumption, Eq. (79) becomes
where is the number of different sequences in the system. When we insert this into Eq. (78) and solve for the steady state, we find
The expressions involving sums over matrix elements in Eq. (81) can be identified with the excess production of the master,
and with the average excess production without the master,
if W has the standard form . Therefore, Eq. (81) corresponds to the often quoted result
However, Eq. (81) is more general in that it can be used even if W is not given as .
The idea here is to insert the monodromy matrix into Eq. (81) in order to obtain an approximation for in the case of periodic landscapes. But under what circumstances can we expect this to work? After all, Eq. (81) has been derived from an equation with continuous time, Eq. (78), whereas the monodromy matrix advances the system in discrete time steps, as can bee seen in Eq. (43). The important point is here that we are only interested in the asymptotic state, which is given by the normalized Perron vector of the monodromy matrix, whether we use discrete or continuous time. Therefore, we are free to calculate the asymptotic state in a periodic landscape for a given phase from
even if this equation does not have a direct physical meaning for finite times. The asymptotic molecular concentrations are then given by the limit of
When we neglect the backflow onto the master sequence, the 0th component of that equation becomes identical to Eqs. (78) and (79), but with the matrix instead of W. This shows that we may indeed use Eq. (81) as an approximation for the asymptotic concentration of . Of course, since we have neglected mutational backflow, this approximation works only for landscapes in which a single sequence has a significant advantage over all others. But this restriction does equally apply to the static case. Numerically, we have found that Eq. (81) works well for a single oscillating peak, and that it breaks down in other cases as expected.