Increase in Complexity and Information through Molecular Evolution

Abstract

Biological evolution progresses by essentially three different mechanisms: (I) optimization of properties through natural selection in a population of competitors; (II) development of new capabilities through cooperation of competitors caused by catalyzed reproduction; and (III) variation of genetic information through mutation or recombination. Simplified evolutionary processes combine two out of the three mechanisms: Darwinian evolution combines competition (I) and variation (III) and is represented by the quasispecies model, major transitions involve cooperation (II) of competitors (I), and the third combination, cooperation (II) and variation (III) provides new insights in the role of mutations in evolution. A minimal kinetic model based on simple molecular mechanisms for reproduction, catalyzed reproduction and mutation is introduced, cast into ordinary differential equations (ODEs), and analyzed mathematically in form of its implementation in a flow reactor. Stochastic aspects are investigated through computer simulation of trajectories of the corresponding chemical master equations. The competition-cooperation model, mechanisms (I) and (II), gives rise to selection at low levels of resources and leads to symbiontic cooperation in case the material required is abundant. Accordingly, it provides a kind of minimal system that can undergo a (major) transition. Stochastic effects leading to extinction of the population through self-enhancing oscillations destabilize symbioses of four or more partners. Mutations (III) are not only the basis of change in phenotypic properties but can also prevent extinction provided the mutation rates are sufficiently large. Threshold phenomena are observed for all three combinations: The quasispecies model leads to an error threshold, the competition-cooperation model allows for an identification of a resource-triggered bifurcation with the transition, and for the cooperation-mutation model a kind of stochastic threshold for survival through sufficiently high mutation rates is observed. The evolutionary processes in the model are accompanied by gains in information on the environment of the evolving populations. In order to provide a useful basis for comparison, two forms of information, syntactic or Shannon information and semantic information are introduced here. Both forms of information are defined for simple evolving systems at the molecular level. Selection leads primarily to an increase in semantic information in the sense that higher fitness allows for more efficient exploitation of the environment and provides the basis for more progeny whereas understanding transitions involves characteristic contributions from both Shannon information and semantic information.

1. Introduction

Biological evolution progressed and progresses by at least two different mechanisms: (i) Optimization of properties that are relevant for reproduction efficiency, which is measured in terms of fitness; and (ii) Major transitions that lead to the next higher hierarchical level of organismic complexity. Optimization gives rise to the enormous diversity of species and subspecies that mirrors the heterogeneity of environments. Evolutionary optimization follows the Darwinian mechanism of variation and selection [1]. This mechanisms is neither confined to living organisms nor to cells. In vitro evolution of molecules that are capable of replication and mutation, commonly RNA or DNA, follows essentially the same principles [2,3] and became an interesting tool of biotechnology in the design of biomolecules with predefined properties [4,5]. The hierarchical structure of biology and its objects is the most prominent result of major evolutionary transitions. In other words, the fact that we can recognize molecules, supermolecular aggregates, cell organelles, cells, organs, organisms, and societies in nature are the strongest hints that evolutionary transitions have led from one level of complexity to the next. Different major transitions share certain features but in general follow different mechanisms [6]. We illustrate by means of an example: Both transitions, (i) the transition from independent genes to a genome; and (ii) the transition from independent organisms to a symbiontic unit have in common that previously independent elements are integrated into the higher hierarchical unit and give up their autonomy—in total or partially—but the mechanism by which this is achieved is entirely different. In the first case, a sufficiently elaborate ligase joins the genes whereas the second case is built upon specific catalysis as will be discussed later on.

What can be said about the expected changes in information during the two evolutionary processes? First we need to define what features of information we are interested in here and it is appropriate to specify precisely the facets of the diverse meanings of information that we shall require. As outlined in Section 8 it is sufficient for our purpose to consider the simplest possible concepts: Shannon information, which increases with the coding capacity of the information carrier increases, and semantic information based on one or more evaluation criteria. In the simplest case of evolutionary relevance the obvious criterion is fitness being measurable in the form of the number of offspring or the rate parameter of growth. Darwinian optimization is accompanied by non-decreasing fitness corresponding to non-decreasing semantic information. Changes in Shannon information may occur but they are uncorrelated to changes in fitness. We shall encounter cases where increase in fitness is a result of decreasing chain length of RNA molecules [2] (see Section 8). The change in information during major transitions in more subtle. Although the relation between chain lengths and complexity of organisms is obscured by many side effects and anything but obvious even for phylogenetically close by lying species, there is no doubt that genome lengths have increased grosso modo in evolution. The coding capacity of organisms increases in major transitions where two or more carriers of genetic information come together to form a cooperative unit. The bacteria with the largest genomes after all are only partly overlapping in DNA length with the smallest protozoa and non-overlapping with multicellular eukaryotes [7]. At the same time the new units that were formed through integration of two or more simpler elements must compete successfully with their subunits as individuals and this implies that semantic information has increased in major transitions too (see Section 8).

A kind of minimal model (Figure 1 and [8,9]), which is derived from chemical kinetics of replication and which contains all essential features that are necessary for describing both optimization of fitness and transition to a higher hierarchical level, is presented in Section 2. At the same time the model is sufficiently simple to be accessible to mathematical analysis and stochastic simulation. The molecular players are understood as replicators being RNA or DNA molecules that are capable of replication, correct or error-prone in the sense of the sketch shown in Figure 2. By assuming suitable changes in the environmental conditions consisting of tunable accessibility of resources the model allows for a description of evolution on two hierarchical levels as well as the transition from the lower to the higher level. In particular, the occurrence of a transition is bound to the availability of sufficiently rich resources and it is possible to define a threshold below which cooperation between competitors does not happen and then no transition will take place. In accordance with the toy model we shall also consider the most simple concepts measuring changes in information: For Shannon information the contents of information is given by the chain length ν of the replicator. In other words we do not distinguish between coding regions, regulatory regions and other regions on the sequence, and we do not account for redundancies. Semantics assigns values to individual replicators, and consistent with the model the evaluation process uses as measure of semantic information the fitness of the phenotype as expressed by the replication rate parameter.

At present evolutionary optimization is fairly well understood. In fortunate cases the occurrence of optimization can be inferred from the fossil record but the inference is always indirect and the conditions under which it occurred are highly complex. Essentially the same is true for field observations in nature. The principle of Darwinian evolution, however, can be studied in evolution experiments with viruses [10] or bacteria [11,12] as well as in various cell-free in vitro assays. As mentioned already, in the first series of test tube experiments reported on evolution, which have been carried out with wild-type RNA of the Escherichia coli bacteriophage Qβ by Sol Spiegelman and his group [2,13], optimization of fitness was achieved through reduction of the virus genome to the absolutely necessary elements for reproduction: (i) a binding site for the enzyme Qβ-replicase; and (ii) a sufficiently rich secondary structure of both strands—plus and minus—in order to avoid duplex formation [14]. This reduction of the genome in the optimization process is made possible by the reaction medium, which supplies everything that is required for reproduction—except a template. Despite being the result of a complex many step template induced polymerization mechanism [15], RNA replication in case of excess resources and excess replicase follows a simple exponential growth characteristic and hence can be modeled by a simple autocatalytic reaction of the class $\mathrm{A}+\mathrm{X}\to 2\mathrm{X}$. The interpretation of the kinetic data was supported by extensive computer calculations. One basic feature of biology is that all units are obligatory replicators. If an element fails to reproduce, it does not show up in future generations. The reproducing units contain genetic information in the form of DNA, only in case of viruses it might also be RNA serving the same purpose. Asexual reproduction comes closest to simple autocatalysis that will be used for reproduction kinetics in the minimal model. Important experimental achievements in the design of an RNA world as a plausible intermediate scenario for the origin of life were: (i) the discovery of catalysis by RNA molecules called ribozymes [16,17,18]; (ii) the design of autocatalytic networks of RNA molecules [19,20]; and (iii) the creation of an RNA based RNA replicase [21]. Experimental evolution found soon extensive applications in biotechnology: The SELEX method—selection of ligands by exponential enrichment—developed in the nineteen-eighties in the labs of Larry Gold and Jack Szostak [22,23] became an indispensable tool of molecular biology laboratories. In evolutionary biotechnology a whole repertoire of different techniques has been conceived and implemented for the design of biomolecules, for RNA-molecules [4] and for proteins [5].

Major transitions span a very broad spectrum from events described in origin-of-life models to the beginnings of human societies. One particular major transition, the origin of the eukaryotic cell is chosen for the purpose of illustration, because it is better understood than most other transitions and the toy model applied here comes relatively close to it. According to the interpretation of the fossil record exclusively prokaryotic life has covered the Earth for about 1.8 billion years, from 3.5 billion years to 1.7 billion years ago. Eukaryotic cells are not only much larger than prokaryotic cells, they have a much more complex structure with so-called cell organelles that according to the generally accepted hypothesis of endosymbiosis [24,25] were independently living organisms before the transition event. In the eukaryotic cell we have two or three reproducing units: (i) the nucleus; (ii) the mitochondria; and (iii) the chloroplasts in case of plant cells. The nucleus contains the major amount of the DNA of the cell but both mitochondria and chloroplasts have their own specific DNA, which encodes some of the organelle-specific proteins. Symbiosis of nucleus and organelle is, of course not the only form of symbiosis [26]. There are numerous examples of two or three-way symbiosis. Four member and higher symbiontic communities, however, are apparently very rare [27] (p. 71). We shall come back to this fact in Section 4, Section 7, and Section 9. Symbiosis characterizes the mutualistic interaction between species that is mutually beneficial. In the mutualistic system competition between elements is suppressed or, in other words, the individual components cooperate. For the minimal model we introduce cooperation in the form of catalytic action during replication, which in its simplest version has been formulated as catalytic hypercycle [28].

The third process that is basic to biological evolution is variation of phenotypic properties through changes in genotypes. It comes essentially in two forms: recombination and mutation. Consistent with the other simplifications of the minimal model presented here we consider exclusively point mutations, which are mutations in the simplest form: A point mutation changes the sequence of a DNA or an RNA molecules at a single position and already this minimal modification introduces a wide spectrum of possible variations into the phenotype that may range from no change to entirely new phenotypic properties. Introduction of mutation leads to formation of mutant distributions that become stationary after sufficiently long time. These stationary solutions are called quasispecies and have been extensively analyzed in the past [29,30,31,32]. They consist of a fittest variant, called the master sequence surrounded by its most frequent mutants weighted by fitness and by the distance to the master. Quasispecies evolve through the formation of fitter variants through mutation and restructuring of the mutant distribution around a new master sequence. The process of quasispecies evolution is visualized best as optimization or adaptive walk of a population on a fitness landscape. All three basic processes of biological evolution are are sketched in Figure 1.

2. A Minimal Model for Competition, Cooperation, and Mutation

The competition-cooperation-mutation (CCM) model [8,9] considers independent and catalyzed, correct and error-prone reproduction in a population of replicators or subspecies ${\mathrm{X}}_i$ with $i=1,2,\cdots ,n$. Stationary populations resulting from correct replication and mutation will be denoted as quasispecies and accordingly we shall use the term subspecies for individual replicators, which will be here always RNA or DNA molecules. A flow reactor (CSTR) is chosen as open system. The flow reactor is fed from a large reservoir filled with stock solution containing the material A at concentration $[\mathrm{A}]=a_0$ where A stands for the whole set of materials that are required for reproduction. The mechanisms of reproduction and catalyzed reproduction in the flow reactor is encapsulated in $(3n+2)$ chemical reaction equations:

$$*\stackrel{a_{0}r}{\to }\mathrm{A},$$

(1a)

$$\mathrm{A}+{\mathrm{X}}_i\stackrel{k_i{Q}_{ji}}{\to }{\mathrm{X}}_i+{\mathrm{X}}_j;i,j=1,...,n,$$

(1b)

$$\mathrm{A}+{\mathrm{X}}_i+{\mathrm{X}}_{i+1}\stackrel{l_i{Q}_{ji}}{\to }{\mathrm{X}}_i+{\mathrm{X}}_j+{\mathrm{X}}_{i+1};i,j=1,...,n;i\mathrm{mod}n,$$

(1c)

$$\mathrm{A}\stackrel{r}{\to }Ø,\mathrm{and}$$

(1d)

$${\mathrm{X}}_i\stackrel{r}{\to }Ø;i=1,...,n.$$

(1e)

Reaction (1a) supplies the material required for reproduction. A solution with A at concentration $a_0$ flows into a continuously stirred tank reactor (CSTR) with a flow rate r [35] (p. 87ff.) We remark that the deterministic kinetic Equation (4) were extensively studied under the simplifying assumption of constant population size ${\sum }_{i=1}^n{x}_i=c_0=\mathrm{const}$ [29,31,36]. The solution curves formulated in relative concentrations ${\mathsf{\xi }}_i(t)=x_i(t)/{\sum }_{i=1}^n{x}_i(t)$ are identical for the CSTR and for constant population size but the stochastic system is unstable in the latter case [37]. The reactor operates at constant volume and this implies that the volume per time unit [t] of solution flowing into the reactor is compensated exactly by an outflow, which is described by the Equation (1d) and (1e) and concerns all $(n+1)$ molecular species, A and ${\mathrm{X}}_i,i=1,\cdots ,n$. More sophisticated flow reactors called chemostats and other experimental devices are also used to monitor and regulate available resources [38,39,40,41,42].

The notation for the rate parameters is chosen such that ambiguity is avoided: $k_i$ and $l_i$ with dimensions [M${}^{-1}$t${}^{-1}$] and [M${}^{-2}$t${}^{-1}$] refer to reactions in the flow reactor. For fitness and the cooperation parameter we use $f_i=k_{i}a$ and $h_i=l_{i}a$, which have the dimensions [t${}^{-1}$] and [M${}^{-1}$t${}^{-1}$]. Reactions (1b) model the processes of correct copying and mutation of replicators (Figure 2) and reactions (1c) finally include reproduction catalyzed by other replicators: ${\mathrm{X}}_{i+1}$ (with $imodn$) supports the process of copying ${\mathrm{X}}_i$ yielding either a correct copy of ${\mathrm{X}}_i$ or a mutant ${\mathrm{X}}_j$. The terms (1c) introduce cooperation between otherwise competing subspecies. In the molecular interpretation ${\mathrm{X}}_{i+1}$ acts as a specific catalyst for the reproduction process $\mathrm{A}+{\mathrm{X}}_i\to 2{\mathrm{X}}_i$. In principle, we could have $n^2$ instead of n catalytic terms in (1c), because every molecule might act as catalyst in the reproduction of every other molecule but efficient and specific catalysis is rare and a meaningful model should not have more terms than required. Not all catalytic networks of replicators are stable in the sense that all subspecies survive in a cooperative organization. The presumably smallest stable functional organization is a hypercycle [36] and requires only n catalysts for n reactions: A sequence of catalytic interactions is closed to a loop,

$$\cdots \leftarrow {\mathrm{X}}_n\leftarrow {\mathrm{X}}_1\leftarrow {\mathrm{X}}_2\leftarrow \cdots \leftarrow {\mathrm{X}}_{n-1}\leftarrow {\mathrm{X}}_n\leftarrow \cdots ,$$

(2)

and this introduces mutual dependence of all members of the hypercycle. If one subspecies, ${\mathrm{X}}_k$, dies out, the time-derivative of the previous species in the cycle, ${\mathrm{X}}_{k-1}$, becomes negative, ${\mathrm{X}}_{k-1}$ vanishes and the entire organization collapses.

The molecular concept behind mutation in Equation (1b) is sketched in Figure 2: A template is bound to an enzyme initiating thereby a copying process that leads either to a correct copy, $\mathrm{A}+{\mathrm{X}}_i\to 2{\mathrm{X}}_i$, or a mutant, $\mathrm{A}+{\mathrm{X}}_i\to {\mathrm{X}}_j+{\mathrm{X}}_i$. Equation (1c) describes the analogous reactions with subspecies ${\mathrm{X}}_{i+1}$ ($imodn$) being a specific catalyst for the reaction. The mutation rates are summarized in the mutation matrix $\mathrm{Q}=\{Q_{ij};i,j=1,\cdots ,n\}$ where $Q_{ji}$ is the frequency with which subspecies ${\mathrm{X}}_j$ is obtained as an error copy of the template ${\mathrm{X}}_i$. Since the copying processes leading to the different subspecies are parallel reaction channels, which taken together exhaust all possibilities—each copy has to be either correct or incorrect—we have the conservation relation ${\sum }_{j=1}^n{Q}_{ji}=1$. Consistent with a minimal model is the assumption of uniform mutation rates: The mutation or error rate p is assumed to be independent of the nature of the mutated nucleotide residue as well as its position along the sequence. Then all elements of the mutation matrix $\mathrm{Q}$ are obtained from three parameters only,

$$Q_{ij}(p)=Q{\epsilon }^{d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)}\mathrm{with}Q={(1-p)}^{\nu }=Q_{ii}\forall i=1,\cdots ,n\mathrm{and}\epsilon =\frac{p}{1-p}.$$

(3)

The Hamming-distance $d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)$ counts the minimum number of consecutive point mutations required to produce ${\mathrm{X}}_i$ as an error copy of ${\mathrm{X}}_j$.

Mechanism (1) is readily translated into kinetic differential equations with the molecular concentrations, $a=[\mathrm{A}]$ and $x_i=[{\mathrm{X}}_i]\forall i=1,\cdots ,n$, as variables:

$$\begin{array}{c}\frac{\mathrm{d}a}{\mathrm{dt}}=-a\sum _{j=1}^n{x}_j\left(k_j+l_j{x}_{j+1}\right)+r(a_0-a)\mathrm{and}\hfill \end{array}$$

(4a)

$$\begin{array}{cc}& \frac{\mathrm{d}{x}_i}{\mathrm{dt}}=a\left(\sum _{j=1}^n{Q}_{ij}(k_j+l_j{x}_{j+1})x_j\right)-r{x}_i;i,j=1\cdots ,n;jmodn,\hfill \end{array}$$

(4b)

where we used the conservation relation ${\sum }_{j=1}^n{Q}_{ji}=1$ in Equation (4b).

The competition-cooperation-mutation (CCM) model in the flow reactor has three internal and two external parameters. According to Figure 1 the three internal parameters are the fitness difference $\Delta f\ge 0$, the cooperation parameter $h\ge 0$, and the mutation rate $p\ge 0$. The external parameters are the concentration of resources $\mathrm{A}$ in the stock solution, $a_0>0$, and the flow rate, $r\ge 0$. The concentration of the material for reproduction in the stock solution, $a_0$, is a measure of the available resources and determines the carrying capacity of the system. The flow rate r sets a limit to the time needed for reproduction events, because it is the reciprocal mean residence time of the solution in the reactor, $r={\tau }_R^{-1}$, and $r=0$ allows for the approach towards thermodynamic equilibrium.

3. Competition and Cooperation

The dynamics of competition and cooperation with zero mutation is described by trajectories in plane $\mathcal{A}$ (Figure 1). The mutation process is neglected simply by setting $p=0$ leading to $\mathrm{Q}=\mathbb{I}$ where $\mathbb{I}$ is the unit matrix. The equation for $a(t)$ (4a) remains unchanged whereas the second equation simplifies to

$$\frac{\mathrm{d}{x}_i}{\mathrm{dt}}=x_i\left(a(k_i+l_i{x}_{i+1})-r\right);i=1\cdots ,n;imodn,$$

(4b’)

For the low dimensional cases, $n=2,3$, the kinetic equations are particularly simple and allow for straightforward mathematical analysis of steady states and bifurcations. Higher dimensional systems, $n\ge 4$, can also be handled analytically but the expressions become clumsy with increasing numbers of stationary states.

Steady state analysis reveals the long-time behavior of the system that is determined exclusively by stationary states, which are related by simple bifurcations, in particular transcritical and saddle-node bifurcations. Properly we start from Equation (4b’) and find for the stationary states:

$$\frac{\mathrm{d}{x}_i}{\mathrm{dt}}=0={\overline{x}}_i\left(\overline{a}(k_i+l_i{\overline{x}}_{i+1})-r\right).$$

(4c)

Equation (4c) sustains two solutions for the steady concentrations of every subspecies: (i) extinction of ${\mathrm{X}}_i$ given by ${\overline{x}}_i^{(i)}=0$; and (ii) survival of ${\mathrm{X}}_i$ is expressed by ${\overline{x}}_i^{(ii)}>0$ and requires

$$\overline{a}(k_i+l_i{\overline{x}}_{i+1}^{(ii)})-r=0\mathrm{and}{\overline{x}}_i^{(ii)}=\frac{1}{l_{i-1}}\left(\frac{r}{\overline{a}}-k_{i-1}\right),imodn.$$

These two solutions for each stationary replicator concentration can be combined in principle to $2^n$ different steady states unless there are cases of incompatibility. As a matter of fact we have $2^n+1$ solutions since the concentration $\overline{a}$ is obtained by means of a quadratic equation with two solutions. One of the two roots is always unstable [8,44]. For $n=2$ and $k_2>k_1$ three out of five states have regions of asymptotic stability in the $(a_0,r)$-plane: (i) the state of extinction $S_0$, $\overline{a}=a_0,{\overline{x}}_1={\overline{x}}_2=0$; (ii) the state of selection of ${\mathrm{X}}_2$, $S_1^{(2)}$ with $\overline{a}=r/k_2,{\overline{x}}_1=0,{\overline{x}}_2=a_0-r/k_2$; and (iii) the state of cooperation $S_2$, $\overline{a}=\alpha ,{\overline{x}}_1=(r-k_2\alpha )/l_2\alpha ,{\overline{x}}_2=(r-k_1\alpha )/l_1\alpha$. Here, the concentration value $\alpha =\overline{a}$ is the root with the minus sign of the quadratic Equation (5)

$$\overline{a}=\alpha =\frac{1}{2}\left(a_0+\psi -\sqrt{{(a_0+\psi )}^2-4r\varphi }\right),$$

with $\varphi =l_1^{-1}+l_2^{-1}$ and $\psi =k_1/l_1+k_2/l_2$. The second root corresponds to an unstable state $S_3$. The state at which the only the less fit variant is present, $S_1^{(1)}$, is unstable too.

For constant resources expressed by $a_0=\mathrm{const}$ and initial conditions $x_i(0)>0\forall i=1,\cdots ,n$ the long-time state of the system is uniquely defined. In

Table 1. Asymptotically stable stationary states of Equation (4) with $\mathit{n}\mathbf{=}\mathbf{2}$ [9]. The three stationary states are ordered with respect to increasing $a_0$-values of their asymptotically stable regimes. Validity of the relations $k_1<k_2$ and $l_1>l_2$ between rate parameters was assumed. For the cooperative state $S_2$ stationary concentration of A is obtained as the root with the negative sign of the quadratic Equation (5) with two combinations of the rate parameters, $\psi =k_1/l_1+k_2/l_2$ and $\varphi =1/l_1+1/l_2$. The existence of the cooperative state requires in addition a sufficiently small flow rate: $r\le {(a_0+\psi )}^2/4\varphi$.

Name	Symbol	Stationary Values			Stability Range
		$\overline{\mathit{a}}$	${\overline{\mathit{x}}}_1$	${\overline{\mathit{x}}}_2$
extinction	$S_0$	$a_0$	0	0	$0\le a_0\le \frac{r}{k_2}$
selection	$S_1^{(2)}$	$\frac{r}{k_2}$	0	$a_0-\frac{r}{k_2}$	$\frac{r}{k_2}\le a_0\le \frac{r}{k_2}+\frac{k_2-k_1}{l_1}$
cooperation	$S_2$	α	$\frac{r-k_2\alpha }{l_2\alpha }$	$\frac{r-k_1\alpha }{l_1\alpha }$	$\frac{r}{k_2}+\frac{k_2-k_1}{l_1}\le a_0$

and

Table 2. Asymptotically stable stationary states of Equation (4) with $\mathit{n}\mathbf{=}\mathbf{3}$ [9]. The four stationary states are ordered with respect to increasing $a_0$-values of their asymptotically stable regime. The relations $k_1<k_2<k_3$ and $l_1>l_2>l_3$ between the rate parameters were assumed. For the cooperative state $S_3$ the stationary concentration of A is obtained as one root of a quadratic Equation (5) with two combinations of the rate parameters, $\psi ={\sum }_{i=1}^3{k}_i/l_i$ and $\varphi ={\sum }_{i=1}^{3}1/l_i$. The existence of the cooperative state requires a sufficiently small flow rate: $r\le {(a_0+\psi )}^2/4\varphi$.

Name	Symbol	Stationary Values				Stability Range
		$\overline{\mathit{a}}$	${\overline{\mathit{x}}}_1$	${\overline{\mathit{x}}}_2$	${\overline{\mathit{x}}}_3$
extinction	$S_0$	$a_0$	0	0	0	$0\le a_0\le \frac{r}{k_3}$
selection	$S_1^{(3)}$	$\frac{r}{k_3}$	0	0	$a_0-\frac{r}{k_3}$	$\frac{r}{k_3}\le a_0\le \frac{r}{k_3}+\frac{k_3-k_2}{l_2}$
exclusion	$S_2^{(1)}$	$\frac{r}{k_3}$	0	$a_0-\frac{r}{k_3}-\frac{k_3-k_2}{l_2}$	$\frac{k_3-k_2}{l_2}$	$\frac{r}{k_3}+\frac{k_3-k_2}{l_2}\le a_0\le \frac{r}{k_3}+\frac{k_3-k_2}{l_2}+\frac{k_3-k_1}{l_1}$
cooperation	$S_3$	α	$\frac{r-k_3\alpha }{l_3\alpha }$	$\frac{r-k_1\alpha }{l_1\alpha }$	$\frac{r-k_2\alpha }{l_2\alpha }$	$\frac{r}{k_3}+\frac{k_3-k_2}{l_2}+\frac{k_3-k_1}{l_1}\le a_0$

we list all stable stationary states for the cases $n=2$ and $n=3$. The stationary coordinates of the states, $\mathbf{S}=(\overline{a},{\overline{x}}_1,{\overline{x}}_2)$ or $\mathbf{S}=(\overline{a},{\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3)$ with $\overline{a}(a_0,r)$, ${\overline{x}}_1(a_0,r)$, ${\overline{x}}_2(a_0,r)$ and ${\overline{x}}_3(a_0,r)$, are either constant or depend linearly on $a_0$ and/or r for all states except the cooperative state. In other words, the stationary states move along straight lines or are fixed in concentration space. Transitions between the individual solutions occur through transcritical bifurcations: The straight lines of two states, ${\mathbf{S}}_i$ and ${\mathbf{S}}_j$, cross at the bifurcation points defined by ${\mathbf{S}}_i\left({(a_0)}_{\mathrm{cr}},r_{\mathrm{cr}}\right)={\mathbf{S}}_j\left({(a_0)}_{\mathrm{cr}},r_{\mathrm{cr}}\right)$ and the stationary states exchange stability properties. The cooperative state $S_n$ is accompanied by an unstable state $S_{n+1}$. For increasing flow rates r the expression under the square root vanishes and both states disappear at the critical value $r_{\mathrm{cr}}={(a_0+\psi )}^2/4\varphi$.

The system with $n=3$ sustains an additional stable stationary state between selection and cooperation where two species are present. We choose for the kinetic parameters: $k_1<k_2<k_3$ and $l_1>l_2>l_3$. This choice simplifies the discussion to some extent because no distinction of subcases is required. We obtain ${\overline{x}}_1=0$ and denote the state $S_2^{(1)}$ as state of exclusion of ${\mathrm{X}}_1$. For higher dimensionality, $n>3$, further states with three and more coexisting species and other species missing will be observed. Eventually we find the following ordering of stability ranges of stationary or quasi-stationary states with increasing resources expressed as $a_0$-values:

$$\mathrm{extinction}\leftrightarrow \mathrm{selection}\leftrightarrow \mathrm{exclusion}\leftrightarrow \cdots \leftrightarrow \mathrm{cooperation}.$$

The state of cooperation is an asymptotically stable stationary state for $n=3$ as it is for $n=2$ provided the flow rate is below the critical value $r_{\mathrm{cr}}={(a_0+\psi )}^2/4\varphi$.

In the pure hypercycle equation—Equation (4) with $k_i=0\forall i=1,\cdots ,n$, $a=\mathrm{const}$, and ${\sum }_{i=1}^n{x}_i=c_0=\mathrm{const}$—weakly damped oscillations around an asymptotically stable stationary state $S_4$ are observed for $n=4$. For $n\ge 5$ the cooperative state $S_n$ is unstable [45] and the ODE (1) sustains a stable limit cycle [46]. As we shall show in the next section stochasticity introduces instability into hypercycle dynamics of cooperative systems with $n\ge 4$. In the CSTR the stationary concentration of A at the state of cooperation is not constant but results from the quadratic equation

$${\overline{a}}_{1,2}=\frac{1}{2}\left(a_0+\psi \mp \sqrt{{(a_0+\psi )}^2-4r\varphi }\right)\mathrm{with}\psi =\sum _{i=1}^n\frac{l_i}{k_i}\mathrm{and}\varphi =\sum _{i=1}^n\frac{1}{k_i}.$$

(5)

Interestingly, Equation (5) depends on n only through the summations in ψ and ϕ, and hence the equation is valid for arbitrary n. The second solution belongs to an unstable state $S_{n+1}$ (see e.g., [44] (pp. 674–675)) and is important in the context of a saddle-node bifurcation at which the states $S_n$ and $S_{n+1}$ disappear above the critical flow rate $r>r_{\mathrm{cr}}={(a_0+\psi )}^2/4\varphi$ [9].

The model for cooperation in the CSTR by means of hypercycles deserves additional attention because it is accessible to full analytical treatment. As in the case of constant concentrations of a and $c_0$ the nature of the solution curves and stationary states depends crucially on the number of molecular species n. Since the dimensions of the rate parameters are $[k_j]=$ [M${}^{-1}$t${}^{-1}$] and $[l_j]=$ [M${}^{-2}$t${}^{-1}$], respectively, the choice of sufficiently high concentrations $a_0$ allows for the neglect of terms containing $k_j$ against those with $l_j·x_{j+1}$ in Equation (4b’). Then we are dealing with pure hypercycle dynamics in the flow reactor as it has been studied analytically before [44]. The mechanism and the kinetic differential equations are obtained from Equations (1) and (4) by putting $k_j=0\forall j=1,\cdots ,n$. With a simple transformation of variables, $\mathbf{x}\Rightarrow \mathbf{z}$, with $z_{j+1}=l_j{x}_{j+1}\forall j,(j,j+1modn)$, which has been denoted characteristically as barycentric transformation [47,48,49], the stationary point inside the concentration simplex ${\mathbb{S}}_n^{(c)}=\{x_j;j=1,\cdots ,n,{\sum }_{i=1}^n{x}_i=c\}$ can be shifted into the center. From $\overline{a}·{\overline{z}}_j-r=0$ follows that all stationary concentrations have the same value ${\overline{z}}_j=\overline{z}=r/\overline{a}$, and $\overline{a}$ can be calculated from Equation (5) by setting $\psi =0$. Then, for $\overline{z}$ we obtain two solutions of the quadratic equation:

$${\overline{z}}_{1,2}=\frac{1}{2\varphi }\left(a_0\pm \sqrt{a_0^2-4r\varphi }\right)$$

(6)

The solution with the plus sign corresponds to the solution with the minus sign for A, ${\overline{a}}_1$. The stability of the central stationary point follows from the eigenvalues of the Jacobian matrix:

$$\begin{array}{cc}\hfill {\left({\lambda }_0\right)}_{1,2}& =\overline{z}(a_0^2-2\overline{z}\varphi )=\mp \sqrt{a_0^2-4r\varphi }\mathrm{and}\hfill \\ \hfill {\lambda }_j& =rexp\left(\frac{2\pi i}{n}j\right),j=1,\cdots ,n-1.\hfill \end{array}$$

(7)

The two solutions ${\left({\lambda }_0\right)}_{1,2}$ correspond to the two steady states $S_n$ and $S_{n+1}$, respectively, and inspection of Equation (7) shows that $S_{n+1}$ is always unstable because of the plus sign, ${\left({\lambda }_0\right)}_2>0$. The eigenvalue ${\left({\lambda }_0\right)}_1$ is always negative and stability or instability is therefore determined by the roots ${\lambda }_j$, which represent the complex roots of one except ${\lambda }_n=1$ that does not exist as a solution here. Accordingly, $S_n$ is asymptotically stable for $n=2$ and $n=3$. The case $n=4$ is special: The linearization around the stationary point $S_4$ is a marginally stable center with two eigenvalues with zero real part, ${\lambda }_{1,3}=\pm \frac{i}{2}r$. There is however a higher order term, which renders the state $S_4$ asymptotically stable. Integration of the differential Equation (4) with $n=4$ and $k_j=0\forall j=1,\cdots ,4$ yields indeed weakly damped oscillations [45]. For $n\ge 5$ the central point $S_n$ is unstable and the system sustains a stable limit cycle [46], which manifests itself through undamped oscillations.

The ODEs (4) combine competition and cooperation kinetics and are appropriate candidates for modeling transitions from Darwinian selection to stable symbiosis. At the same time it seems that nothing in the model can be omitted without loosing the capability to describe such a transition and therefore we intend to conjecture that it is also the simplest model possible. Despite its simplicity the model provides a straightforward explanation for the emergence of cooperation in a world of competitors triggered by increasing resources.

4. Stochastic Kinetics in the Competition-Cooperation System

The most natural way for considering stochasticity in chemical reactions and reaction networks is to write down and analyze the chemical master equation corresponding to mechanism (1). In order to distinguish particle numbers and concentrations we use upper case letters for former: $[\mathrm{A}]=A(t)$ and $[{\mathrm{X}}_i]=X_i(t)(i=1,\cdots ,n)$. The probabilities are expressed by $P_m(t)=\mathrm{Prob}\left(A(t)=m\right)$ and $P_{s_i}(t)=\mathrm{Prob}\left(X_i(t)=s_i\right)$ with $m,s_i\in \mathbb{N}$, and finally the multivariate master equation of the competition-cooperation system is of the form

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{P}_{\mathbf{m}}}{\mathrm{dt}}& =a_{0}r{P}_{(\mathbf{m};m-1)}+r\left((m+1)P_{(\mathbf{m};m+1)}+\sum _{j=1}^n(s_j+1)P_{(\mathbf{m};s_j+1)}\right)+\hfill \\ & +(m+1)\sum _{j=1}^n\left((k_j+l_j{s}_{j+1})(s_j-1)P_{(\mathbf{m};m+1,s_j-1)}\right)-\hfill \\ & -\left(r\left(a_0+m+\sum _{j=1}^n{s}_j\right)+m\left(\sum _{j=1}^n\left(k_j+l_j{s}_{j+1}\right)s_j\right)\right)P_{\mathbf{m}}.\hfill \end{array}$$

(8)

The index set is denoted by $\mathbf{m}=(m,s_1,\cdots ,s_n)$ and $(\mathbf{m};m-1)=(m-1,s_1,\cdots ,s_n)$, etc. Since the derivation of the master equation assumes that chemical processes occur one at a time all jumps involve single elementary steps and the reactions in the mechanism (1) change particle numbers by $\pm 1$, we apply here a notation in shorthand for changes of the index set for reaction $\mathbf{m}\to {\mathbf{m}}^{\prime }$,

$$\begin{array}{cc}\hfill {\mathbf{m}}^{\prime }& =({\mathbf{m}}^{\prime }=m\pm 1,s_1,\cdots ,s_n)\equiv (\mathbf{m};m\pm 1)\mathrm{for}\mathrm{Equation}(1\mathrm{a},\mathrm{d}),\hfill \\ \hfill {\mathbf{m}}^{\prime }& =({\mathbf{m}}^{\prime }=m,s_1,\cdots ,s_k-1,\cdots ,s_n)\equiv (\mathbf{m};s_k-1)\mathrm{for}\mathrm{Equation}(1\mathrm{e})\mathrm{and}\hfill \\ \hfill {\mathbf{m}}^{\prime }& =({\mathbf{m}}^{\prime }=m-1,s_1,\cdots ,s_k^{\prime }=s_k+1,\cdots ,s_n)\equiv (\mathbf{m};m-1,s_k+1)\mathrm{for}\mathrm{Equation}(1\mathrm{b},\mathrm{c}).\hfill \end{array}$$

Although it is not difficult to write down a multivariate master equation, the derivation of analytical solutions is successful only in exceptional cases, for example for networks of monomolecular reactions [50,51]. An alternative strategy for studying chemical master equations is computer simulation through sampling of trajectories. The theoretical background for trajectory harvesting has been laid down by Andrey Kolmogorov [52], Willy Feller [53], and Joe Doob [54,55]. With the incoming of electronic computers, simulations of stochastic processes became possible. The conception, analysis, and implementation of a simple and highly efficient algorithm by Daniel Gillespie [56,57,58] provided a very useful tool for investigations of stochastic effects in chemical kinetics. Here we present some results, which illustrate the differences between deterministic and stochastic solutions of the transition model.

A typical trajectory with $n=2$ is shown in Figure 3. In order to facilitate the analysis we start with an empty reactor, $A(0)=0$, that contains only the autocatalytic molecules in seeding quantities. These initial conditions allow for the distinction of four phases of the stochastic process. The inflow of stock solution with $[\mathrm{A}]=a_0$ results in a rapid increase of the number density of A molecules (phase I). In contrast to deterministic reaction kinetics knowledge of parameter values and initial conditions are not sufficient for correct predictions of the final state: In phase II random events guide the system either towards the absorbing state of extinction ($S_0$) or to one of the three quasi-stationary states ($S_1^{(1)}$, $S_1^{(2)}$ or $S_2$). In phase III the system approaches the long-time state and eventually in phase IV we observe the trajectories fluctuating around the deterministic values of the corresponding stationary states ($S_0$, $S_1^{(1)}$, $S_1^{(2)}$ or $S_2$).

The existence and uniqueness conditions for the solutions of ODEs apply to the kinetic Equation (4b) and we are dealing with a unique long-time behavior for precisely defined parameters and initial conditions. In contrast to the solution of a kinetic differential equation the outcome of stochastic trajectories with identical parameters and initial conditions need not be unique. Indeed the example of the trajectory shown in Figure 3 shows a representative case of one class—the class converging to $S_2$—whereas trajectories may converge to any of the four stationary states, the absorbing state $S_0$ and the three quasi-stationary states $S_1^{(1)}$, $S_1^{(2)}$, and $S_2$. By trajectory sampling we calculate the probabilities for ending up in a particular state. In

Table 3. Probabilities to reach states in the cooperative regime with $\mathit{n}\mathbf{=}\mathbf{2}$ and different initial conditions. The table provides counts of approach towards the four final states: extinction $S_0$, selection of ${\mathrm{X}}_1$ in $S_1^{(1)}$, selection of ${\mathrm{X}}_2$ in $S_1^{(2)}$, or cooperation $S_2$. The given values are sample means and unbiased standard deviations calculated from ten packages, each of them containing 10,000 trajectories computed with identical parameters and initial conditions, and differing only in the sequence of random events determined by the seeds of the pseudorandom number generator (Extended CA, Mathematica). Choice of parameters: $k_1=0.09$ [M${}^{-1}$t${}^{-1}$], $k_2=0.11$ [M${}^{-1}$t${}^{-1}$], $l_1=0.0011$ [M${}^{-2}$t${}^{-1}$], $l_2=0.0009$ [M${}^{-2}$t${}^{-1}$], $a_0=200$, $r=0.5$ [V t${}^{-1}$]. Initial conditions: $A(0)=0$, $X_1(0)$ and $X_2(0)$ see table.

Initial Values		Counted Numbers of States in Final Outcomes
${\mathit{X}}_{\mathbf{1}}(\mathbf{0})$	${\mathit{X}}_{\mathbf{2}}(\mathbf{0})$	${\mathit{N}}_{{\mathit{S}}_{\mathbf{0}}}$	${\mathit{N}}_{{\mathit{S}}_{\mathbf{1}}^{(\mathbf{1})}}$	${\mathit{N}}_{{\mathit{S}}_{\mathbf{1}}^{(\mathbf{2})}}$	${\mathit{N}}_{{\mathit{S}}_{\mathbf{2}}}$
1	1	$385.1\pm 23.6$	$1481.0\pm 36.8$	$1719.6\pm 37.8$	$6414.3\pm 53.8$
2	1	$77.4\pm 9.1$	$1822.6\pm 41.6$	$367.6\pm 17.0$	$7733.3\pm 38.3$
1	2	$71.6\pm 8.5$	$280.6\pm 20.0$	$2075.8\pm 28.9$	$7572.0\pm 39.2$
3	1	$15.0\pm 2.9$	$1900.4\pm 30.9$	$74.6\pm 10.0$	$8009.0\pm 35.3$
1	3	$14.0\pm 3.7$	$53.1\pm 4.8$	$2180.5\pm 48.4$	$7752.3\pm 53.8$
2	2	$14.9\pm 2.6$	$303.7\pm 16.0$	$354.5\pm 23.8$	$9326.8\pm 44.9$
3	3	0	$70.2\pm 10.0$	$106.2\pm 10.9$	$9823.4\pm 15.7$
4	4	0	$12.1\pm 2.6$	$28.0\pm 5.0$	$9959.9\pm 6.4$
5	5	0	$2.5\pm 1.1$	$6.3\pm 2.6$	$9991.2\pm 3.0$

we present the results for counting final states for different initial conditions and sets of parameter values. Seeding values of five molecules ${\mathrm{X}}_1$ and ${\mathrm{X}}_2$ each and a sufficiently large value of $a_0$ are enough to obtain almost certain convergence of the stochastic trajectories to the cooperative state $\left(P(S_2)>0.999\right)$.

The interpretation of the individual values in Table 3 is straightforward. The more seeding molecules of a given species we have, the less likely it is eliminated during the random phase II. We observe selection of both replicators ${\mathrm{X}}_1$ and ${\mathrm{X}}_2$, and the presence in more copies overrules the difference in fitness values—at least in the example shown here. One might argue: “When particle numbers of five are already sufficient to come close to the deterministic results why worry about stochastic effects at all?” This argument were valid, were there not the fact that every mutant inevitably has to begin with a single copy, and in the early evolutionary phase of a newly created molecular species stochasticity is always important. Interestingly, the probability to become extinct $P(S_0)=N_{S_0}/\mathrm{10,000}$ depends primarily on the sum of the seeding molecules, ${\mathrm{X}}_1(0)+{\mathrm{X}}_2(0)$, and not so much on the distribution over the two initial states. The interpretation is straightforward: The probability that all ${\mathrm{X}}_j$ molecules ($j=1,2$) are diluted out of the reactor before they are replicated does not depend on the particular subspecies, ${\mathrm{X}}_1$ or ${\mathrm{X}}_2$.

Stochasticity in competition-cooperation dynamics of systems with $n\ge 4$ shows additional features that are not present in the deterministic solutions of the kinetic equations. The interpretation is easier if one considers first pure hypercycle dynamics, $k_j=0\forall i=1,\cdots ,n$, for $n\ge 5$. The kinetic equations sustains two states: (i) extinction $S_0$ with $[\mathrm{A}]=\overline{a}=a_0$ and $[{\mathrm{X}}_j]={\overline{x}}_j=0$ ($j=1,\cdots ,n$); and (ii) cooperation consisting of an unstable state $S_n$ surrounded by a stable closed orbit resulting in oscillating variables ${\mathrm{X}}_j(t)$ [46]. The cooperative state is unstable in the stochastic process as well. Strong uncompensated fluctuations undergo autocatalytic self-enhancement for $n\ge 5$, the stochastic variables ${\mathrm{X}}_j(t)$ exhibit noisy oscillations and the amplitudes of the oscillations grow with time until one subspecies, say ${\mathrm{X}}_k$, dies out. Then we have $\mathrm{d}{x}_{k-1}/\mathrm{dt}<0$ and accordingly ${\mathrm{X}}_{k-1}$ vanishes followed by decreasing concentrations of ${\mathrm{X}}_{k-2}$ until it dies out, then ${\mathrm{X}}_{k-3}$ vanishes and so on. In the pure hypercyclic system the whole population dies out with the disappearance of ${\mathrm{X}}_{k+1}$ ($kmodn$). The sequence of extinction obviously is the same as the sequence of hypercyclic catalysis (2). The four dimensional system ($n=4$) is more subtle: At small population sizes the situation is the same just as described for $n\ge 5$ but in large populations it may take very long time before the concentration of one subspecies becomes zero and then we observe oscillation of the concentrations of all subspecies with fluctuating frequencies and amplitudes. Since every new species in biological evolution has to start ultimately from very small population sizes ($N=1$), we may expect that in reality cooperative systems with four members are unstable like the larger ones.

In essence, the addition of first order autocatalytic reaction terms, $k_{j}a{x}_j\forall j=1,\cdots ,n$, changes the dynamics only in one aspect (Figure 4): Frequent possible outcomes are not only the two states extinction $S_0$ and the cooperative state $S_n$ but also the n selection states $S_1^{(j)}$ ($j=1,\cdots ,n$). The last remaining subspecies in the sequence of extinction from above, ${\mathrm{X}}_{k+1},kmodn$, and A together form a quasi-stationary state. Since every subspecies may be wiped out first, every possible selection state $S_1^{(j)},j=1,\cdots ,n$, can result as final outcome.

In summary, stochasticity introduces two major changes into the competition-cooperation system: (i) Fluctuations allow for the choice between several final states of the stochastic process with identical parameters and initial conditions; (ii) cooperative systems with five and more partners are unstable and lead to random selection of one subspecies or extinction; and (iii) the case $n=4$ is intermediate in the sense that the quasi-stationary cooperative state $S_4$ is long-lived as with $n=2$ and $n=3$ for sufficiently large population sizes but unstable in small populations like the systems with $n\ge 5$.

5. Competition, Mutation and Quasispecies

Optimization of mean fitness may occur through competition and mutation in evolving populations. The competition-mutation system (coordinate plane $\mathcal{B}$) is studied in the CSTR and the equations result from mechanism (1) through neglect of second order autocatalysis (1c), $l_j=0\forall j=1,\cdots ,n$. In order to illustrate the problem with the simplest possible case we begin with correct replication, $p=0$. Then we are left with the kinetic differential equations

$$\begin{array}{cc}& \frac{\mathrm{d}a}{\mathrm{dt}}=-a\sum _{i=1}^n{f}_i{x}_i+r(a_0-a)\hfill \end{array}$$

(9a)

$$\begin{array}{cc}& \frac{\mathrm{d}{x}_j}{\mathrm{dt}}=f_{j}a{x}_j-r{x}_j=x_j(f_{j}a-r),j=1,\cdots ,n.\hfill \end{array}$$

(9b)

Equation (9) sustains two stationary states: (i) the state of extinction, $S_0$, characterized by $\overline{a}=a_0$ and ${\overline{x}}_j=0\forall j=1,\cdots ,n$; and (ii) the state resulting from natural selection, $S_1^{(m)}$, of the uniquely defined subspecies with highest fitness, ${\mathrm{X}}_m$ with $f_m=max\{f_i;i=1,\cdots ,n\}$, since we do not consider neutral evolution [34] here where two or more subspecies have identical largest fitness values. The state $S_1^{(m)}$ fulfils the conditions

$$f_m-r=k_m\overline{a}-r=0\mathrm{or}\overline{a}=\frac{r}{k_m}\mathrm{and}{\overline{x}}_m=a_0-\frac{r}{k_m}.$$

(9c)

In order to reveal the relation between selection and fitness optimization we consider the mean replication rate $\widehat{k}(t)={\sum }_{i=1}^n{k}_i{x}_i(t)/c(t)$ with $c(t)={\sum }_{i=1}^n{x}_i(t)$ and its time derivative,

$$\frac{\mathrm{d}\widehat{k}}{\mathrm{dt}}=\sum _{i=1}^n\frac{\partial \widehat{k}}{\partial x_i}·\frac{\mathrm{d}{x}_i}{\mathrm{dt}},\frac{\partial \widehat{k}}{\partial x_i}=\frac{1}{c}(k_i-\widehat{k})\mathrm{leading}\mathrm{to}\frac{\mathrm{d}\widehat{k}}{\mathrm{dt}}=a(t)(\widehat{k^2}-\widehat{k}{}^2)=a(t)\mathrm{var}(k).$$

(9d)

The concentration of A is positive, $a(t)>0$ and the variance of any distribution is a nonnegative quantity too, $\mathrm{var}(k)\ge 0$, where the equal sign requires a homogeneous population, $x_m=c$ and $x_i=0\forall i=1,\cdots ,n;i\ne m$, and hence $\mathrm{d}\widehat{k}/\mathrm{dt}\ge 0$ is always fulfilled. It is, of course, possible to choose $a(0)=0$ but then we have $\mathrm{d}a/\mathrm{dt}=r{a}_0>0$ and in the next instant $a(t)>0$ would be fulfilled. Equation (9d) states that the mean fitness is always increasing except at its maximum $\widehat{k}=k_m$, which is the stationary state $S_1^{(m)}$. Hence the mean fitness of populations is optimized during the process of selection in the flow reactor.

Equation (9b) can be rewritten and interpreted differently [59] (pp. 29–32): For every replicator ${\mathrm{X}}_j$ we define a specific growth function ${\mathrm{\Gamma }}_j(\mathbf{x})$ that describes unconstrained growth, $\mathrm{d}{x}_j/\mathrm{dt}={\mathrm{\Gamma }}_j(\mathbf{x})$, and an unspecific constraint $\mathsf{\varphi }(t)$, which is the same for all replicators:

$$\frac{\mathrm{d}{x}_j}{\mathrm{dt}}={\mathrm{\Gamma }}_j\left(\mathbf{x}(t)\right)-\frac{x_j(t)}{c(t)}\mathsf{\varphi }(t).$$

(9e)

In the flow reactor, for example, the constraint fulfils: $\mathsf{\varphi }=c(t)r$. In more elaborate reactors as far as the experimental implementation is concerned other constraints lead to expressions that are more complicated to write down but often easier to analyze. Summation over all subspecies yields an ODE for the total concentration of replicators:

$$\frac{\mathrm{d}c}{\mathrm{dt}}=\sum _{i=1}^n{\mathrm{\Gamma }}_j\left(\mathbf{x}(t)\right)-\mathsf{\varphi }(t),$$

(9f)

which can be used to calculate $c(t)$ from known $\mathsf{\varphi }(t)$ and vice versa:

$$c(t)=c(0)+{\int }_0^t\mathrm{d}\tau \left(\sum _{i=1}^n{\mathrm{\Gamma }}_i(\mathbf{x})-\mathsf{\varphi }(\tau )\right)\mathrm{or}\mathsf{\varphi }(t)=\sum _{i=1}^n{\mathrm{\Gamma }}_i(\mathbf{x})-\frac{\mathrm{d}c}{\mathrm{dt}}.$$

Stationarity in the total concentration, $\mathrm{d}c/\mathrm{dt}=0$, allows for an expression of the constraint in terms of the growth functions: $\mathsf{\varphi }(t)={\sum }_{i=1}^n{\mathrm{\Gamma }}_i\left(\mathbf{x}(t)\right)$. This constraint—often called constant concentration or constant organization—is particularly useful for the analysis of replicator equations [36,45,60]. Insertion into Equation (9e) yields

$$\frac{\mathrm{d}{x}_j}{\mathrm{dt}}={\mathrm{\Gamma }}_j\left(\mathbf{x}(t)\right)-\frac{x_j(t)}{{\sum }_{i=1}^n{x}_i(t)}\sum _{i=1}^n{\mathrm{\Gamma }}_i\left(\mathbf{x}(t)\right).$$

(9e’)

Next we introduce normalized concentrations, ${\mathsf{\xi }}_j=x_j/{\sum }_{i=1}^n{x}_i=x_j/c$, with ${\sum }_{i=1}^n{\mathsf{\xi }}_i=1$ and $\mathbf{x}=c·\mathsf{\xi }$ into Equation (9e) and the ODE for $\mathsf{\xi }$ becomes

$$\frac{\mathrm{d}{\mathsf{\xi }}_j}{\mathrm{dt}}=\frac{1}{c(t)}\left({\mathrm{\Gamma }}_j\left((c(t)\mathsf{\xi }\right)-{\mathsf{\xi }}_j\mathsf{\varphi }(t)\right).$$

Summation of normalized concentrations yields ${\sum }_{i=1}^n\mathrm{d}\mathsf{\xi }/\mathrm{dt}=0$. Accordingly we find for the constraint $\mathsf{\varphi }(t)={\sum }_{i=1}^n\mathrm{\Gamma }\left(c(t)\mathsf{\xi }\right)$ and obtain

$$\frac{\mathrm{d}{\mathsf{\xi }}_j}{\mathrm{dt}}=\frac{1}{c(t)}\left({\mathrm{\Gamma }}_j\left(c(t)\mathsf{\xi }\right)-{\mathsf{\xi }}_j\sum _{i=1}^n\mathrm{\Gamma }\left(c(t)\mathsf{\xi }\right)\right).$$

(9g)

It is worth noticing that Equation (9g) does not contain explicitly the constraint $\mathsf{\varphi }(t)$ and hence it is valid for almost all cases: The evolution of the population described in normalized variables is independent of the constraint as long as the population size does neither become zero nor approach infinity. An implicit dependence of the rate Equation (9g) on the constraint $\mathsf{\varphi }(t)$ is nevertheless given through the concentration $c(t)$. In case the growth functions ${\mathrm{\Gamma }}_j$ are homogeneous functions of degree λ in x, we find

$${\mathrm{\Gamma }}_j(\mathbf{x})={\mathrm{\Gamma }}_j(c\mathsf{\xi })=c^{\lambda }{\mathrm{\Gamma }}_j(\mathsf{\xi }),$$

and the ODE in normalized concentrations takes on the form:

$$\frac{\mathrm{d}{\mathsf{\xi }}_j}{\mathrm{dt}}=c^{\lambda -1}\left({\mathrm{\Gamma }}_j\left(\mathsf{\xi }\right)-{\mathsf{\xi }}_j\sum _{i=1}^n\mathrm{\Gamma }\left(\mathsf{\xi }\right)\right);j=1,\cdots ,n.$$

(9h)

Since ${\sum }_{i=1}^n{\mathsf{\xi }}_i=1$ Equation (9h) and (9e’) are identical apart from the factor $c{(t)}^{\lambda -1}$, which can be absorbed in the time axis as long as $c>0$ and $c<\infty$ are fulfilled. Two cases can be distinguished: (i) $\lambda \ne 1$, the stationary states of both equations are the same and so are the trajectories on the concentration simplex but the solution curves differ by the time factor $c{(t)}^{\lambda -1}$; and (ii) $\lambda =1$, the factor containing time and total concentration is one and thus time independent, and the two equations do not only have the same stationary states but also identical solution curves. In other words the course of competition and selection is the same in stationary and growing populations. Darwinian optimization is based on first order autocatalysis characterized by homogeneous growth functions with $\lambda =1$. Hypercycle dynamics uses homogeneous growth functions with $\lambda =2$ and the internal equilibrium in the growing system is identical with the stationary state approached at constant concentration. We remark that the growth functions in the competition-cooperation system, ${\mathrm{\Gamma }}_j(\mathbf{x})=(k_j+l_j{x}_{j+1})a{x}_j$ with $jmodn$, are not homogeneous and the regularities reported in this paragraph do not apply therefore.

Evolutionary optimization requires more than selection. Variation of genotypes is a conditio sine qua non for the Darwinian mechanism and in the spirit of the as simple as possible CCM-model we introduce point mutation (Figure 2) as source of variation. In order to study competition-mutation dynamics we set $l_j=0\forall j=1,\cdots ,n$ in the mechanism (1) and obtain the kinetic differential equations

$$\begin{array}{cc}& \frac{\mathrm{d}a}{\mathrm{dt}}=-a\sum _{i=1}^n{k}_i{x}_i\sum _{j=1}^n{Q}_{ji}+r(a_0-a)=-a\sum _{i=1}^n{k}_i{x}_i+r(a_0-a)\hfill \end{array}$$

(10a)

$$\begin{array}{cc}& \frac{\mathrm{d}{x}_j}{\mathrm{dt}}=a\sum _{i=1}^n{Q}_{ji}{k}_i{x}_i-r{x}_j,j=1,\cdots ,n.\hfill \end{array}$$

(10b)

The growth functions ${\mathrm{\Gamma }}_j=a{\sum }_{i=1}^n{Q}_{ji}{k}_i{x}_i$ are homogenous and linear in $\mathbf{x}$ and hence the equations in conventional and normalized concentrations, $\mathbf{x}(t)$ and $\mathsf{\xi }(t)$, are the same, and have identical solutions. The stationarity condition, $\mathrm{d}c/\mathrm{dt}=0$, yields again two solutions (i) the state of extinction $S_0=(\overline{a}=a_0,\overline{c}=0)$ and (ii) the state of quasispecies selection at $S_1=\left(\overline{a}=r/\overline{k},\overline{c}=a_0-r/\overline{k}\right)$ with $\overline{k}={lim}_{t\to \infty }\widehat{k}(t)={lim}_{t\to \infty }{\sum }_{i=1}^n{k}_j{x}_j(t)/c(t)$ [30]. The two states are related by a transcritical bifurcation at $a_0\overline{k}=r$. The state of selection, $S_1$, is asymptotically stable for small flow rates, or sufficiently large resources $a_0>r/\overline{k}$. Extinction occurs at high flow rates, $r>a_0\overline{k}$, and small resources $a_0<r/\overline{k}$. The problem is not yet solved completely because the stationary mutant distribution, the quasispecies $\overline{\mathsf{{\rm Y}}}=({\overline{x}}_i;i=1,\cdots ,n)$, is needed for the calculation of $\overline{k}$. The quasispecies is conventionally obtained through the solution of an eigenvalue problem and we shall briefly sketch this procedure here. Other approaches make use of techniques developed in the statistical mechanics of magnetic systems [61,62,63,64] (for a recent update see [65]).

Because of Equation (9h) with $\lambda =1$ the solutions approaching the state of quasispecies selection are the same at constant or variable population sizes and without loosing generality the replication-mutation problem can be and has been solved at constant resources $a=\mathrm{const}=\overline{a}$ as well as constant population sizes $c=\mathrm{const}=\overline{c}$ [29,30,31] (for a recent review see [32]). The constant concentration of A is absorbed in the fitness value $f_j=k_j\overline{a}$. We mention here only the most prominent results, which are relevant for the forthcoming discussion of changes in the information content of a population. The quasispecies is obtained in form of the largest eigenvector of the value matrix $\mathrm{W}$,

$$\mathrm{W}=\{W_{ji}=Q_{ji}{f}_i\}=\mathrm{Q}·\mathrm{F}\mathrm{with}\mathrm{Q}=\left\{Q_{ji}\right\}\mathrm{and}\mathrm{F}=\{F_{ij}=f_i{\delta }_{ij}\};i,j=1,\cdots ,n,$$

(10c)

which is a product of the mutation matrix $\mathrm{Q}$ and a diagonal matrix $\mathrm{F}$ containing the fitness values of all subspecies. The eigenvalue problem is of the form

$$\mathrm{W}·{{\mathsf{\zeta }}_k}^{\mathrm{t}}={\lambda }_k{{\mathsf{\zeta }}_k}^{\mathrm{t}}\mathrm{with}{\mathsf{\zeta }}_k=({\mathsf{\zeta }}_1^{(k)},\cdots ,{\mathsf{\zeta }}_n^{(k)});k=0,1,\cdots ,n-1\mathrm{and}\sum _{i=1}^n{\mathsf{\zeta }}_i^{(k)}=1.$$

(10d)

Vectors are assumed to be row vectors here and, in particular, the symbol ‘${}^{\mathrm{t}}$’ means transposed and indicates a column vector.

The matrix $\mathrm{W}$ is either positive or at least nonnegative and irreducible, the conditions for the applicability of the Perron-Frobenius theorem [66] are fulfilled, and hence the largest eigenvalue ${\lambda }_0$ is a non-degenerate eigenvalue and the corresponding eigenvector, the quasispecies $\overline{\mathsf{{\rm Y}}}={\mathsf{\zeta }}_0$, which represents the long-time solution of the competition-mutation problem, has exclusively positive components. In other words, all subspecies are present at the stationary state and no mutant vanishes in the process of natural selection combined with error-prone reproduction. The most frequent subspecies in the quasispecies is called the master sequence, ${\mathrm{X}}_m$. Often but not always it is the variant with the highest selective value, $W_{mm}=Q_{mm}{f}_m=max\{W_{ii};i=1,\cdots ,n\}$. In case the replication accuracy is the same for all subspecies, $Q_{ii}=Q\forall i=1,\cdots ,n$—as it occurs, for example, with the uniform error rate assumption—the sequence with the largest selective value is identical to the sequence with highest fitness. As we shall show in Section 6 this need not necessarily mean that the fittest sequence or the sequence with the largest selective value is the master sequence. The quasispecies consists of the master sequence and a mutant cloud surrounding it where the width of the cloud depends on the distribution of fitness values and the mutation rates.

How is Darwinian natural selection changed by the inclusion of mutations? The answer is readily given in mathematical terms. Frequent mutations couple individual subspecies to clans that are selected together and the clan, which reproduces with maximal efficiency is the dominant eigenvector ${\mathsf{\zeta }}_0$ of the value matrix $\mathrm{W}$. Eigenvectors corresponding to the other eigenvalues are less efficient since we have

$${\lambda }_0>{\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_n,$$

and after sufficiently long time only the term containing ${\lambda }_0$ is important, since the time dependent weighting factors of the contribution of the eigenvectors ${\mathsf{\zeta }}_k$ are: $z_k(t)\propto exp(-{\lambda }_{k}t)$. All eigenvectors except ${\mathsf{\zeta }}_0$ have positive and negative components and are positioned outside the physically reachable concentration space defined by nonnegative concentrations $x_i\ge 0\forall i=1\cdots ,n$. Population dynamics can be visualized as a process in the space spanned by the eigenvectors ${\mathsf{\zeta }}_k\forall k=0,\cdots ,n-1$:

$$\mathsf{{\rm Y}}(t)=\left(x_1(t),\cdots ,x_n(t)\right)=\sum _{i=1}^n{x}_i(t){\mathit{e}}_i=\sum _{k=0}^{n-1}{z}_k(t){\mathsf{\zeta }}_k,$$

(10e)

where ${\mathit{e}}_i$ is a Cartesian eigenvector in the direction of ${\mathrm{X}}_i$. Indeed rewriting replication-mutation dynamics in terms of the eigenvalues and eigenvectors of $\mathrm{W}$ yields a kinetic differential equation that looks identical to the mutation-free case with the fittest subspecies ${\mathrm{X}}_m$ being replaced by the quasispecies $\overline{\mathsf{{\rm Y}}}={\mathsf{\zeta }}_0$. The stationary solutions are defined by ${lim}_{t\to \infty }{x}_j(t)/c(t)={\zeta }_j^{(0)}$ and are—as required—independent of time and initial conditions. The stationary mean fitness, $\overline{f}={\sum }_{i=1}^n{f}_i{\zeta }_i^{(0)}$, is the maximal mean fitness, which the population can achieve at mutational equilibrium. It fulfils necessarily $\overline{f}\le f_m$, where the equal sign corresponds to no mutation and holds for a homogeneous population of master sequence. In other words, the maximal fitness $f_m$ can be obtained only for vanishing mutation rates $p\to 0$. Obviously, the optimization theorem of the mean fitness derived for error-free replication is not valid any more and trajectories along which mean fitness is decreasing or non-monotonously changing are easily found. According to Equation (10e) the time dependence of the population is given by a superposition of exponential functions with the eigenvalues ${\lambda }_k$. After sufficiently long time—when the system is close to the stationary state—only the largest eigenvalue ${\lambda }_0$ and the corresponding eigenvector ${\mathsf{\zeta }}_0$ are important.

Exact solutions in closed form are not available but phenomenological expressions for the purpose of illustration can be derived through three simplifying assumptions [29,32]:

(i)

A single-peak fitness landscape is applied that assumes equal fitness for all mutants and a higher fitness for the master sequence (Section 6).

(ii)

A uniform mutation rate per site and replication event, p, is assumed. In other words the frequency of mutation is assumed to be independent of the nature and the position of the mutated nucleotide. The mutation matrix is largely simplified by the uniform error rate assumption

$$Q_{ij}(p)=Q{\epsilon }^{d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)}\mathrm{with}Q={(1-p)}^{\nu }=Q_{ii}\forall i=1,\cdots ,n\mathrm{and}\epsilon =\frac{p}{1-p}.$$

The Hamming-distance $d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)$ counts the minimum number of point mutations needed to produce ${\mathrm{X}}_i$ as an error copy of ${\mathrm{X}}_j$.

(iii)

Mutational backflow in the kinetic differential Equation (10) is neglected. We rewrite Equation (10b) for $j=m$ and partition in two contributions coming from correct copying of the template ${\mathrm{X}}_m$ and from incorrect copying of all other ${\mathrm{X}}_i$ with $i\ne m$, and neglect the second term:

$$\frac{\mathrm{d}{x}_m}{\mathrm{dt}}=a\left(Q_{mm}{f}_m{x}_m+\sum _{i=1,i\ne m}^n{Q}_{mi}{f}_i{x}_i\right)-r{x}_m\approx x_m(Q_{mm}{f}_{m}a-r).$$

(10b’)

For small mutation rates p, ignoring backflow is an appropriate approximation. Neglect of backflow, in other words, means that we solve the ODEs without the terms in which the master sequence is produced from a mutant by mutation, and obtain for the stationary solutions [29]:

$${\overline{\mathsf{\xi }}}_m\approx \frac{Q-{\sigma }_m^{-1}}{1-{\sigma }_m^{-1}}\mathrm{and}{\overline{\mathsf{\xi }}}_j\approx {\epsilon }^{d_{\mathrm{H}}({\mathrm{X}}_j,{\mathrm{X}}_m)}\frac{Q-{\sigma }_m^{-1}}{{(1-{\sigma }_m^{-1})}^2}.$$

(11a)

The fact that individual fitness values $f_j$ do not enter the Equation (11a) is a results of the assumption of a single peak fitness landscape.

The stationary frequency at which a given subspecies ${\mathrm{X}}_j$ is present in the quasispecies is a function the mutation rate p and the degree of relatedness to the master sequence expressed by the Hamming distance. There is, of course, also a dependence on the difference in fitness values, $f_m-f_j$, which is encapsulated here in the superiority of the master sequence,

$${\sigma }_m=\frac{f_m}{{\widehat{f}}_{-m}}\mathrm{with}{\widehat{f}}_{-m}=\frac{{\sum }_{i=1,i\ne m}^n{f}_i{\mathsf{\xi }}_i}{1-{\mathsf{\xi }}_m},$$

(11b)

and which will be discussed in detail in Section 6. The quantity ${\widehat{f}}_{-m}$ is the mean fitness of all sequences except the master. Since $\epsilon \ll 1$ the stationary concentration of mutants decreases fast with increasing Hamming distance from the master sequence and the width of the distribution increases with increasing mutation rate p. The most prominent result of quasispecies theory is the existence of a sharply defined error threshold and follows directly from Equation (11a). All components of the quasispecies contain a common factor $Q-{\sigma }_m^{-1}$, which becomes zero at the critical mutation rate $p_{\mathrm{cr}}=1-{\sigma }_m^{-1/\nu }\approx ln{\sigma }_m/\nu$, where the chain length of the polynucleotide template is denoted by ν and ${\sigma }_m$ is the superiority of the master sequence defined above. The approximation applied in Equation (11a) assumes equal fitness of all mutants. This assumption can be relaxed in the sense that all fitness values are different without loosing the existence of a sharp error threshold (see Section 6 and [32]).

The existence of a critical error rate or error threshold $p_{\mathrm{cr}}$ can be interpreted heuristically in straightforward way: Mutation has the consequence that a certain fraction of the copies of the master sequence, $1-Q_{mm}\approx 1-{(1-p)}^{\nu }$, are less fit than the parent subspecies. This fraction apparently increases with increasing mutation rate p whereas the mean fitness of the stationary population decreases. There is a—very high—mutation rate $p=\widehat{p}$ at which the incorporations of correct or incorrect digits are equally probable—for binary sequences this occurs at an error rate of $\widehat{p}=1/2$: At this point all mutants are produced with equal probability and the stationary distribution of subspecies is the uniform distribution. Properly one can characterize reproduction at such a high level of inaccuracy as random replication. If an error threshold exists, the transition from the ordered regime of a structured quasispecies to a uniform distribution produced by random replication occurs to a very good approximation already at much lower error rates: $p_{\mathrm{cr}}\ll \widehat{p}$. In other words, in systems with error thresholds the point of random replication $\widehat{p}=1/2$ is widened to a broad zone $p_{\mathrm{cr}}<p\le \widehat{p}$.

6. Sequence Space, Fitness Landscape, and Population Dynamics

Understanding evolution is facilitated enormously by the application of two fundamental concepts: (i) sequence space; and (ii) fitness landscape. Nucleic acids are visualized as carriers of genetic information, the sequence or genotype space is a point space, and every nucleic acid sequence is represented by a point. The distance between two sequences ${\mathrm{X}}_i$ and ${\mathrm{X}}_j$ is the Hamming metric $d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)$ [29,67,68]. It should be mentioned that the idea of a sequence space for proteins originated about the same time [69].), which represents the minimal number of point mutations that are required in order to interconvert the two sequences. The notion of fitness landscape goes back to the population geneticist Sewall Wright [70,71]: A fitness value is assigned to every point in sequence space and the resulting object is a landscape over a high-dimensional support. Sequence spaces are huge—the number of possible sequences of chain lengths ν over an alphabet with κ digits is ${\kappa }^{\nu }$ and this amounts to $1.3\times {10}^{30}$ for small natural nucleotide sequences of length $\nu =50$. In reality sequence spaces can never be fully covered by populations, which only in exceptional cases can be as large as ${10}^{15}$ individuals. Thermodynamics on the other hand assumes equal distribution of molecular species over all degrees of freedom and the deterministic approach is based on the assumption that concentrations may become arbitrarily small. In chemical kinetics of replication and mutation (Equation (10) and Figure 2) qualitatively the same situation arises by the assumption that all reaction channels, $n={\kappa }^{\nu }$, are populated according to their weighting factors $Q_{ji}$. This is far away from any real situation where we have either one molecule or none and the usage of a discrete model is indispensable. A real population covers only a tiny part of sequence space—a typical distribution of virus genotypes, for example, consists of (i) a master sequence; (ii) a core of frequent mutants, which are present almost all the time; and (iii) mutants at the periphery, which “come and go”. Master sequence and its frequent mutants may be described by the deterministic quasispecies equation restricted to the area of sequence space that is covered by the core. We shall use the term local quasispecies here in order to express the fact that optimization of fitness is restricted to a small region in sequence space. The periphery, accordingly, cannot be modeled properly unless fluctuations are taken into account.

A snapshot of an evolving population will most probably catch a local quasispecies. New mutants are formed and the selective value of the majority of them does not exceed that of the current master sequence ${\mathrm{X}}_m^{(1)}$, $W_{mm}^{(1)}=W_{max}^{(1)}=max\{W_{jj}^{(1)};j=1,\cdots ,n_1\}$. Once in a while a mutant ${\mathrm{X}}_m^{(2)}$ is formed, which has a higher selective value than the current master, $W_{mm}^{(2)}>W_{max}^{(1)}$. Provided the new sequence is not lost during a stochastic initial period, a new local quasispecies centered around ${\mathrm{X}}_m^{(2)}$ with $W_{mm}^{(2)}=W_{max}^{(2)}=max\{W_{jj}^{(2)};j=1,\cdots ,n_2\}$. Computer simulations [72,73,74,75,76] have shown that quasispecies evolution as sketched in Figure 5 follows regularly a two phase process: (i) During the quasi-stationary phases along fitness plateaus the master sequence ${\mathrm{X}}_m^{(k)}$ stays the same, the mean fitness is approximately constant and the population broadens in sequence space; and (ii) a plateau phase ends by the advent of a fitter mutant, the width of the population shrinks instantaneously, the mean fitness increases and a new mutant cloud builds up around the new master sequence ${\mathrm{X}}_m^{(k+1)}$. Such two phase processes—broadening of the population though the formation of a mutant cloud around a master through spreading by stochastic drift, and narrowing of the population as a consequence of the transition to a new master—is repeated over and over again as long as the global optimum $W_{\mathrm{opt}}$ has not been reached. Because of the enormous size of sequence space a typical trajectory will not be able to come even close to the optimum and therefore the evolutionary process is practically open ended in most cases. The difference between the fitness and the selective value of the master sequence, $\Delta f_m=f_m-W_{mm}=f_m(1-Q_{mm})\approx f_m\left(1-{(1-p)}^{\nu }\right)$, is a quantity that depends exclusively on the properties of the master sequence and hence is different from the mutation load that measures the fitness of the master sequence relative to the mean fitness of the population: $L=(f_m-\widehat{f})/f_m$. Both quantities become zero in a homogeneous population of the master sequence. The evolutionary sequence of selective values is determined by the inequality

$${\mathrm{X}}_m^{(1)}\to {\mathrm{X}}_m^{(2)}\to \cdots :W_{mm}^{(1)}<W_{mm}^{(2)}<\cdots <W_{\mathrm{opt}}\le f_{opt}.$$

(12)

Fulfilling the equals sign in the last inequality requires vanishing mutation rates, $limp\to 0$.

How do realistic fitness landscapes look like [32] (pp. 62–75)? The investigation of biopolymer structures and functions as well as extensive works on pathogenic virus populations revealed three generic features of fitness landscapes in recent years [77,78,79,80]: (i) Realistic fitness landscapes are high dimensional—in a polynucleotide sequence of length ν all ν positions can be varied independently and the sequence space is a discrete ν-dimensional object with κ points in every direction; (ii) fitness landscapes are rugged in the sense that a small change in the sequence may cause dramatic fitness changes or no change at all; and (iii) neutrality implying that several or many sequences have the same selective value. Figure 6 sketches an adaptive walk and shows how a neutral segment of the walk may bridge an otherwise unsurmountable obstacle for the walk on a rugged landscape. The sketch at the same time suggests that both neutrality and ruggedness are required for evolution: Without sufficient neutrality adaptive walks would be trapped at some low-lying nearby peaks. If the peak is a member of a neutral network, however, the population can circumvent the trap by random drift in another direction and reach a point from which the adaptive walk can be continued. Ruggedness is a consequence of the generic relations between biopolymer sequences, structures and functions (For details of sequence-structure-function mappings of RNA molecules see [76,81,82]). It creates a multitude of local fitness optima and provides the basis for diversity in nature and adaptation to environmental changes.

Despite spectacular successes in the experimental determination and empirical modeling of fitness landscapes (see, for example, [78]) detailed information on sufficiently large parts of fitness landscapes is still missing. Accordingly, almost all studies were made with more or less simple model landscapes. Many results of quasispecies theory were derived by means of largely simplified landscapes and an important question concerns the general validity of these findings. As representative examples we mention here three landscapes. A very simple but nevertheless frequently used example is the single peak landscape:

$$f_j=\left\{\begin{array}{cc}{f}_0\hfill & \mathrm{for}j=m,\hfill \\ f_1\hfill & \mathrm{for}j=1,\cdots ,{\kappa }^{\nu };j\ne m,\hfill \end{array}\right.$$

(13a)

with $f_0>f_1$. It has only two fitness values, $f_0$ the fitness of the master sequence and $f_1$ the fitness value shared by all sequences ${\kappa }^{\nu }$ except the master. It was used in the calculation of the analytical expression for the error threshold $p_{\mathrm{cr}}$ and we ask now whether or not error thresholds are also found on more general landscapes.

For this goal we construct more realistic landscapes that allow to apply different fitness values for different individual sequences. The lack of sufficiently detailed empirical data forces us to use some random input, which we create by superposition of random scatter upon single peak landscapes. Because we aim at mimicking landscapes that are ultimately derived from biopolymers, such landscapes are called realistic random landscapes (For another class of empirical random landscapes see, for example, the Nk-model [83,84]):

$$f_j=\left\{\begin{array}{cc}{f}_0\hfill & \mathrm{for}j=m,\hfill \\ f_1+2d(f_0-f_1)({\eta }_j^{(s)}-0.5)\hfill & \mathrm{for}j=1,\cdots ,{\kappa }^{\nu };j\ne m,\hfill \end{array}\right.$$

(13b)

where ${\eta }_j(s)$ is the output of a pseudorandom number generator drawing numbers from a uniform distribution on the interval $[0,1]$, $\mathcal{U}:0\le {\eta }_j^{(s)}\le 1$, which had been started with the seeds s. The parameter d allows for tuning the degree of randomness. The value $d=0$ implies the single peak landscape with no random scatter at all and $d=1$ yields a band with fully developed random scatter covering the entire range $2f_1-f_0\le f_j\le f_0$.

The third example is a realistic random landscape with a tunable degree of neutrality, λ. Neutrality is incorporated into realistic random landscapes in straightforward manner: The fitness value $f_0$ is not only assigned to the master sequence but to all sequences ${\mathrm{X}}_j$ with pseudorandom numbers $1\ge {\eta }_j^{(s)}\ge 1-\lambda$ where $0<\lambda <1$ is a tunable degree of neutrality:

$$f_j=\left\{\begin{array}{cc}{f}_0\hfill & \mathrm{if}j=m,\hfill \\ f_0\hfill & \mathrm{if}{\eta }_j^{(s)}\ge 1-\lambda ,\hfill \\ f_1+\frac{2d}{1-\lambda }(f_0-f_1)({\eta }_j^{(s)}-0.5)\hfill & \mathrm{if}{\eta }_j^{(s)}<1-\lambda ,j=1,\cdots ,{\kappa }^{\nu };j\ne m.\hfill \end{array}\right.$$

(13c)

It is easy to see that $\lambda =0$ yields the non-neutral landscape and $\lambda =1$ results in completely flat landscape. Evolution on the flat landscape is described by the neutral theory of evolution developed by Motoo Kimura [34].

One general result derived from rugged fitness landscapes with resolution to individual sequences concerns the existence of an error threshold. Compared to the single peak landscape the position of the threshold is shifted towards smaller p-values with increasing random scatter d and this observation is readily explained by the fact that the difference between $f_0$ and the next highest fitness values is reduced with increasing d [32] (pp. 98–114). There are, however, rather smooth simple landscapes, which do not sustain error thresholds [32,85]. Neutrality, in general, does not prevent the existence of error thresholds. Single master sequences may be replaced by cluster of closely related neutral sequences [32]. Another relevant finding concerns landscapes with high degree of ruggedness, $d>0.9$: Depending on the distribution of fitness values controlled by the seeds of the pseudorandom number generator we observe two cases: (i) strong quasispecies where the master sequence ${\mathrm{X}}_m$ stays the same in the entire range $0\le p<p_{\mathrm{cr}}$ for $0\le d\le 1$, and (ii) standard quasispecies, which undergo one or more transitions ${\overline{\mathsf{{\rm Y}}}}_m\leftrightarrow {\overline{\mathsf{{\rm Y}}}}_k$ where the master sequence changes ${\mathrm{X}}_m\leftrightarrow {\mathrm{X}}_k$ at certain critical mutation rates $p_{\mathrm{tr}}$ [32].

The necessary and sufficient condition for the occurrence of a transition between quasispecies with different master sequences is crossing or avoided crossing of two eigenvalues at the transition point $p_{\mathrm{tr}}$ [86]. The eigenvalues as functions of the mutation rate p are accessible only by numerical calculation but the existence of a transition can be made plausible by inspection of the two kinetic differential equations for the two potential master sequences ${\mathrm{X}}_m$ and ${\mathrm{X}}_k$ (A rigorous derivation of the condition for transitions between quasispecies is found in [87]). We mention also that transitions between quasispecies were rediscovered thirteen years later by Claus Wilke et al. in numerical simulation by means of digital organisms [88] and characterized by the catchphrase survival of the flattest where flatness refers to the fitness landscape:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{x}_m}{\mathrm{dt}}& =\left(x_m{f}_m(d)Q_{mm}(p)+\sum _{i=1,i\ne m}^n{Q}_{mi}(p)f_i(d)x_i\right)-r{x}_m=\hfill \\ \text{(10b′)}\hfill & & =\left(x_m{W}_{mm}(p,d)+{\mathsf{\Phi }}_{m\leftarrow }(p,d)\right)-r{x}_m,\hfill \hfill \frac{\mathrm{d}{x}_k}{\mathrm{dt}}& =\left(x_k{f}_k(d)Q_{kk}(p)+\sum _{i=1,i\ne k}^n{Q}_{ki}(p)f_i(d)x_i\right)-r{x}_k=\hfill \\ \text{(10b″)}\hfill & & =\left(x_k{W}_{kk}(p,d)+{\mathsf{\Phi }}_{k\leftarrow }(p,d)\right)-r{x}_k.\hfill \end{array}$$

Which of the two candidates, ${\mathrm{X}}_m$ or ${\mathrm{X}}_k$, becomes the master sequence depends on the difference between the two differential quotients at the point $x_m=x_k$, as expressed by the two differences $\Delta {\mathsf{\Psi }}_{mk}=W_{mm}(p,d)-W_{kk}(p,d)$ and $\Delta {\mathsf{\Phi }}_{mk}={\mathsf{\Phi }}_{m\leftarrow }(p,d)-{\mathsf{\Phi }}_{k\leftarrow }(p,d)$. If $D_{mk}=\Delta {\mathsf{\Psi }}_{mk}-\Delta {\mathsf{\Phi }}_{mk}$ is positive, ${\mathrm{X}}_m$ is (very likely to become) the master sequence (As said the transition occurs at a crossing or avoided crossing of two eigenvalues ${\lambda }_0$ and ${\lambda }_1$ and the difference discussed here is the leading term in the difference of the eigenvalues, which determines almost always (but not always) the exact behavior.), whereas $D_{mk}<0$ is a very strong indication for ${\mathrm{X}}_k$ being the master sequence. Within the uniform error rate assumption we have $Q_{mm}=Q_{kk}=Q={(1-p)}^{\nu }$. At $p=0$ the two differences have the values $\Delta {\mathsf{\Psi }}_{mk}=Q(f_m-f_k)>0$, and $\Delta {\mathsf{\Phi }}_{mk}=0$ since $Q_{ij}(p)=0$ iff $i\ne j$, $D_{mk}>0$ and hence ${\mathrm{X}}_m$ is the master sequence. Next we increase p at constant d, and Q as well as $\Delta {\mathsf{\Psi }}_{mk}$ become smaller, whereas at the same time the terms containing $Q_{im}=Q{\epsilon }^{d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_m)}$ increase in absolute value. We need to consider only cases where $\Delta {\mathsf{\Phi }}_{mk}<0$ for $p>0$ because otherwise, for $\Delta {\mathsf{\Phi }}_{mk}>0$, $D_{mk}$ will be always positive and no transitions between quasispecies, ${\overline{\mathsf{{\rm Y}}}}_m\leftrightarrow {\overline{\mathsf{{\rm Y}}}}_k$, can occur. For $\Delta {\mathsf{\Phi }}_{mk}<0$ the difference $D_{mk}$ becomes smaller with increasing mutation rates p and it may become negative before the population reaches the error threshold at $p=p_{\mathrm{cr}}$ and then we observe a transition between two different quasispecies.

The role of the intensity parameter d of random scatter of fitness values is readily analyzed. As a reference we consider the case $d=0$, the single peak landscape. Straightforward calculations yield $\Delta {\mathsf{\Psi }}_{mk}=Q(f_m-f_k)=Q(f_0-f_1)$ and $\Delta {\mathsf{\Phi }}_{mk}=Q(f_0-f_1){\epsilon }^{d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)}$, and accordingly we obtain $D_{mk}=Q(f_0-f_1)\left(1-{\epsilon }^{d_{\mathrm{H}}({\mathrm{X}}_i,{\mathrm{X}}_j)}\right)>0$, and no transitions are possible on single peak landscapes. The influence of a distribution of fitness values with ${\overline{f}}_{-m}=f_1$ instead of the single value $f_1$ of the single-peak landscapes can be predicted straightforwardly: Since $f_m=f_0$ in Equation (13b) is independent of the fitness scatter d, and $f_k$, which evidently has to lie above the average ${\overline{f}}_{-m}=f_1$, is increasing with increasing scatter, the difference $f_m-f_k$ will decrease with increasing d. Consequently, the condition for a transition between quasispecies can be fulfilled at lower p-values the larger d is and we expect to find one or more transitions below the error threshold $p_{\mathrm{cr}}$. Numerical calculations show that $d>0.9$ is commonly needed for the occurrence of transitions [32].

7. Cooperation and Mutation

The question to be handled in this section concerns the influence of mutation on a dynamical system showing cooperation between subspecies. The influence of mutation is particularly important in oscillating stochastic systems with high probability of extinction. Intuitively one would guess that high mutation rates could compensate for extinction through reintroduction of the missing subspecies. Whether or not this is the case has been be investigated by computer simulation (Schuster, P., unpublished results, 2016), because mathematical analysis of the hypercycle equation with mutation is rather involved. As we have seen in Section 4 four membered hypercycles with small population size and all hypercycles with five or more members ($n\ge 5$) are endangered by extinction. Self-enhancing stochastic oscillations increase in amplitude until a subspecies vanishes and the whole system dies out. In the two examples shown in Figure 7 mutations prevent the system from extinction several times but the mutation rate parameter p is not large enough to sustain the population with probability one for arbitrarily long time. Figure 8 shows an episode where two subspecies are successfully replaced by mutation.

Cooperation-mutation dynamics is based on the mechanism (1) with $k_j=0\forall j=1,\cdots ,n$ leading to the differential equations

$$\begin{array}{cc}& \frac{\mathrm{d}a}{\mathrm{dt}}=-a\sum _{i=1}^n{l}_i{x}_i{x}_{i-1}\sum _{j=1}^n{Q}_{ji}+r(a_0-a)=-a\sum _{i=1}^n{l}_i{x}_i{x}_{i-1}+r(a_0-a)\hfill \end{array}$$

(14a)

$$\begin{array}{cc}& \frac{\mathrm{d}{x}_j}{\mathrm{dt}}=a\sum _{i=1}^n{Q}_{ji}{l}_i{x}_i{x}_{i-1}-r{x}_j;j=1,\cdots ,n;imodn.\hfill \end{array}$$

(14b)

A mutation matrix $\mathrm{Q}$ with the elements $Q_{ij}$ calculated from the mutation rate parameter p and a distance between the subspecies has to be defined. In case of the quasispecies it has been natural to choose the conventional sequence space as reference and to apply the uniform error rate model (3). Accordingly the appropriate choice for $n=4$ is the binary sequence space with chain length $\nu =2$. For $n=5$ there exists no analogue of a sequence space and we choose a mutation matrix that is built upon a pentagram (see caption to Figure 7). As before in the quasispecies equation the ODE for A does not depend explicitly on the mutation rate p but there is, of course, an implicit dependence via the concentration variables $x_j(t)$. The internal parameters on the cooperation-mutation plane (Figure 1) are the rate parameters $l_j$ or $h_j=l_{j}a$ and the mutation rate p. Since the dynamics of hypercycles remains essentially unchanged in a barycentric transformation we assume equal rate parameters, $l_j=l$ or $h_j=h=la$ ($\forall j=1,\cdots ,n$), without loosing generality. Then we obtain two stationary states: (i) the state of extinction $S_0$ with ${\overline{x}}_j=0\forall j=1,\cdots ,n$ and $\overline{a}=a_0$; and (ii) the cooperative state $S_n$ with ${\overline{x}}_j\ne 0$.

Stability analysis is straightforward for the systems without long-time oscillations ($n<5$). Assuming the existence of a stationary state with equal concentrations $\overline{x}$, which is true for $n=2,3,4$ but not for $n\ge 5$ allows for a straightforward calculation of the stationary concentration of A

$${\overline{a}}^2-a_0\overline{a}+\frac{n}{l}r=0\mathrm{and}{\overline{a}}_{n,n+1}=\frac{1}{2}(a_0\mp \sqrt{a_0^2-{\left(\frac{nr}{l}\right)}^2}.$$

(14c)

This quadratic equation has two solutions, one solution is the cooperative state $S_n$ and the second solution is a saddle point $S_{n+1}$ that separates the basins of attraction of $S_0$ and $S_n$. The situation is in complete analogy to the mutation-free system handled in Equation (5). The state of extinction $S_0$ is always stable, $S_n$ and its satellite $S_{n+1}$ exist for sufficiently small flow rates. If the flow rate exceeds a critical value, $r>r_{\mathrm{cr}}=l{a}_0/n$, the cooperative state $S_n$ and the saddle point $S_{n+1}$ do not exist and the state of extinction $S_0$ provides the only long-time solution. We point again to the fact that the systems with stable closed orbits are more complex and the assumptions made in the derivation of Equation (14c) are not fulfilled.

Stochasticity in the cooperation-mutation system is studied by means of computer simulation with Gillespie’s algorithm. In view of the enormous scatter of extinction times it is quite time consuming to do proper statistics. Therefore we were using an approach needing less computer resources. A set of ten seeds for the pseudorandom number generator was chosen, the corresponding ten trajectories were recorded as functions of mutation rates p, and mean extinction times were calculated for this rather very small sample. Accordingly, the mean values provide only hints on the regularities that would be derivable on the basis more accurate data obtained from larger samples. Nevertheless, the effects are sufficiently large and significant conclusions can be drawn. At first we show for a single trajectory that proper mutations can prevent a population from extinction (Figure 8). If a subspecies dies out and is replaced by mutation during the period of hypercycle decay the population can be saved from extinction. The figure shows on the basis of events seen with a single trajectory that the reintroduction of missing subspecies may be quite complicated and multiple mutations are often required before the oscillatory regime is reestablished. In

Table 4. Extinction times of cooperative populations with $\mathit{n}\mathbf{=}\mathbf{4}$ as a function of the mutation rate parameter $\mathit{p}$. A limit for recording trajectories is set with ${10}^6$ steps. Choice of parameters and initial conditions: $a_0=200$, $r=0.5$ [Vt${}^{-1}$], $l_1=l_2=l_3=l_4=0.1$ [M${}^{-2}$t${}^{-1}$], $A(0)=0$, $X_1(0)=X_2(0)=X_3(0)=X_4(0)=4$. Pseudorandom number generator: Extended CA (Mathematica 10).

	$\mathit{n}=4$: Mutation Rate Parameter p
Seed	0.0	0.001	0.002	0.0025	0.003	0.004	0.005
491d	1212.4	654.5	313.0	636.0	>3000	>3000	1939.4
919	292.4	1069.7	1764.8	2920.1	>3000	>3000	>3000
521	25.3	25.3	23.8	30.0	27.1	43.7	28.4
877	1920.5	367.1	1863.8	308.4	>3000	1940.9	1146.9
233	85.3	71.2	1532.6	>3000	1617.9	84.9	2398.7.4
373	1017.7	482.2	1668.3	525.1	1146.2	1231.5	>3000
089	319.4	338.5	75.3	82.4	81.0	1239.8	81.0
773	161.2	411.2	361.8	2297.2	>3000	>3000	>3000
131	936.2	1059.8	78.0	1362.4	1742.9	>3000	>3000
631	135.0	60.4	670.3	>3000	87.8	884.4	978.4
mean	600	500	800	-	-	-	-

and

Table 5. Extinction times of cooperative populations with $\mathit{n}\mathbf{=}\mathbf{5}$ as a function of the mutation rate parameter $\mathit{p}$. A limit for recording trajectories is set with ${10}^6$ steps. Choice of parameters and initial conditions: $a_0=400$, $r=0.5$ [Vt${}^{-1}$], $l_1=l_2=l_3=l_4=l_5=0.01$ [M${}^{-2}$t${}^{-1}$], $A(0)=0$, $X_1(0)=X_2(0)=X_3(0)=X_4(0)=X_5(0)=5$. Pseudorandom number generator: Extended CA (Mathematica 10).

	$\mathit{n}=5$: Mutation Rate Parameter p
Seed	0.0	0.0005	0.0010	0.0015	0.0020	0.0025
491	36.1	31.3	56.2	118.9	48.1	652.9
919	52.9	42.9	81.6	229.9	546.9	>1650
521	39.0	80.9	96.2	116.1	366.3	210.0
877	34.4	40.6	70.2	310.6	1611.6	1541.6
233	23.8	85.2	34.6	310.2	90.9	>1650
373	41.2	45.5	105.8	680.4	50.1	>1650
089	29.5	28.3	73.3	97.7	831.3	>1650
773	50.7	73.8	85.7	485.2	682.3	>1650
131	21.5	21.5	80.0	104.0	953.1	>1650
631	33.0	78.6	68.9	265.7	298.7	248.0
mean	37	53	75	272	548	-

we put this findings on a more quantitative basis. Times of extinction for oscillating systems with $n=4$ and $n=5$ are calculated by computer simulation as functions of the mutation rate parameter p. In both cases the step sizes in the mutation rate were chosen so small that at least in a single case no mutational effect was observed. In both tables the examples are indicated in bold-face—$s=521$ for $n=4$ and $s=131$ for $n=5$—and then the continuation to larger p-values was done with this step size, which was $\Delta p=0.001$ and $\Delta p=0.0005$ for $n=4$ and $n=5$, respectively. The series of p-values for $n=5$ shows gradual increase from $p=0$ to $p=0.002$. At the next higher p-value, $p=0.0025$, 60% of the trajectories did not reach extinction with the predefined time interval, $0\le t_{\mathrm{ext}}<1650$. For $n=4$ the situation seems to be more complex: A kind of low mutation rate zone is followed by the range of long-living trajectories, which again appear above $p=0.0002$ and amount roughly to the same fraction of trajectories—40%—for higher mutation rates up to $p=0.005$. Conclusions on a mutation threshold phenomenon are presumably premature in view of the small samples.

An interesting detail can be observed by inspection of trajectories in Figure 4 and Figure 7, which concerns the final states after the oscillations have died out. Oscillatory dynamics of competition-cooperation systems can pass over into every selection state, $S_1^{(j)}$ with $j=1,\cdots ,n$, in the first case, whereas pure hypercycles or populations sustaining cooperation and mutation always end up in the empty reactor containing exclusively A. The answer is trivial: All corners of the concentration simplex are unstable in hypercycle dynamics.

8. Information, Transitions and the Competition-Cooperation-Mutation Model

Biological information is a heavily discussed topic and has been the subject of many papers and books (For a recent survey of information and genomic sequences see the special volume of Philosophical Transactions of the Royal Society [89] and in particular the two papers [90,91]). One major issue is to reconcile the enormous complexity of processes in cells or organisms with a simple physical or technical concept resulting originally from the theory of communication via encoded messages [92]. The CCM model introduced here in its simplest form is free from most of the characteristic biological complications. To give an example: For genetic DNA sequences it is important to distinguish between coding, potentially coding and other stretches in the context of information. The model does not deal with translation from DNA to protein and hence no distinction in the sense mentioned above is required. The variation part does neither involve recombination nor more complex mechanisms of mutation that are changing the length of the information carriers. Another related burning issue deals with the origin of complexity in genomes [93], which can be identified with the amount of information on the environment that is stored in the genomic sequence [94]. The CCM model again has an exceedingly simple environment that is fully characterized by the two external parameters $a_0$ and r providing quantitative measures for resources and time constraints. In this special issue of Entropy the relations between information and self-organization are in the focus (See, e.g., [95,96]) and we shall investigate the role of information in the simple model for transitions presented and analyzed here. Consciously, all subtle aspects of information and the question whether or not it is meaningful to define a specific biological information are left aside. Only two crude features of information are considered: (i) syntactic or Shannon information related to the total coding capacity of an idealized genome; and (ii) semantic information understood as the evolutionary value of information as expressed by fitness.

Shannon information is dealing with a discrete random variable Ξ that is defined for a sample space $\mathrm{\Omega }=\{X_1,\cdots ,X_n\}$ with the probability mass function ${\mathsf{\xi }}_i=P(\Xi =X_i)$ for $i=1,\cdots ,n$ and ${\sum }_{i=1}^n{\mathsf{\xi }}_i=1$. Two quantities are of interest here: the content of information $I(\Xi )=-{log}_{2}P(\Xi )$ and its expectation value, the information entropy $H(\Xi )=E\left(I(\Xi )\right)$ [92,97]. For the application to the evolution of a population in the sense of the competition-cooperation-mutation (CCM) model we identify the elements $X_i$ with the subspecies and the variables ${\mathsf{\xi }}_i$ are the corresponding probabilities that a randomly picked individual belongs to class $X_i$. The calculation of the expectation value yields

$$H=\sum _{i=1}^n{\mathsf{\xi }}_i{log}_2\frac{1}{{\mathsf{\xi }}_i}=-\sum _{i=1}^n{\mathsf{\xi }}_i{log}_2{\mathsf{\xi }}_i\mathrm{with}\sum _{i=1}^n{\mathsf{\xi }}_i=1.$$

(15)

The formal identity of Equation (15) with Ludwig Boltzmann’s expression for the thermodynamic entropy explains Shannon’s choice of notion. The subspecies in the CSTR experiment are understood as polynucleotide sequences, RNA or DNA, of chain lengths ν. In our simple model the content of information is completely determined by the chain length ν: Since we are dealing with ${\kappa }^{\nu }$ different sequences of chain length ν the probability to draw one particular sequence from a uniform distribution is $P(\mathsf{\xi })={\kappa }^{-\nu }$ and $I=\nu lo{g}_2\kappa$ [bits]. We are not engaging in coding details in our simple model and the coding capacity of a subspecies is the same as its information content, and hence it is proportional to the chain length ν.

The Shannon entropy (15) is a measure of the broadness of probability distribution of the population being tantamount to the distribution of mutants. For a homogeneous population the sequence that is drawn from it is completely determined a priory and the information we gain by learning which sequence it was is zero as is the entropy: $I=-lo{g}_{2}1=0$ and $H=0$. For the uniform distribution the entropy attains its maximum value $H=\nu lo{g}_2\kappa$ [bits]. In Figure 9 we show the population entropy as a function of the mutation rate p and the bandwidth parameter d of the fitness landscape (13b). On the single peak landscape, $d=0$, and in absence of mutation, $p=0$, the stationary distribution consists of the master sequence ${\mathrm{X}}_m$ only and the information entropy H is zero. The entropy H increases with increasing p until it reaches the maximum value ${log}_{2}n$ corresponding to $\nu lo{g}_2\kappa$ [bits] near the error threshold. With increasing random scatter d or bandwidth of the fitness values $f_j$ the increase in the width of the quasispecies moves towards lower p-values and the error threshold occurs at smaller mutation rates $p_{\mathrm{cr}}$. Individual differences between landscapes with the same random scatter intensity d but different distribution of fitness values $f_j$ become visible only in the range $0.8<d\le 1.0$ (Figure 9). Transitions between quasispecies, ${\mathsf{{\rm Y}}}_m\leftrightarrow {\mathsf{{\rm Y}}}_k$, become detectable as discontinuities in the derivative $\mathrm{d}H/\mathrm{d}p$ on landscapes with large band width of fitness values, $0.9<d\le 1.0$. Eventually we mention that the entropy of the population is also an appropriate tool to illustrate the two-phase process of quasispecies evolution with shrinking and widening population diversity (Section 6, Figure 5, and [72,73,74,75,76]).

The evaluation of the meaning of information as done in semantics is context dependent and commonly very sophisticated. In order to be able to understand and evaluate the message, the receiver requires information that is at least as complex as the information in the message [98]. In the spirit of the model discussed here we propose an approach that is based on simplified RNA structures called secondary structures [76]. The idea is based on the conventional paradigm of structural biology: sequence determines structure and structure determines function. Accordingly, and as sketched in Figure 10 the evaluation of subspecies as carriers of semantic information is done in two steps: (i) RNA sequences ${\mathrm{X}}_j$ are folded into RNA secondary structures ${\mathrm{S}}_j$; and (ii) fitness values $f_j$ are derived from the structures in form of RNA replication rate parameters. Sequence space is discrete by definition and so is structure space. The mapping into function therefore yields also a discrete spectrum of fitness values, which are ready to enter population dynamics. Considering a population of subspecies ${\mathrm{X}}_j$ we can define an expectation value of the evaluation, which in the case considered here is the mean fitness of the population $\mathsf{{\rm Y}}$:

$$W(\mathsf{{\rm Y}})=\widehat{f}(\mathsf{{\rm Y}})=\sum _{i=1}^n{\mathsf{\xi }}_i{f}_i\mathrm{with}\sum _{i=1}^n{\mathsf{\xi }}_i=1.$$

(16)

As an expectation value the semantic value of the population has the same structure as Shannon’s entropy when the information content is replaced by a functional value.

How does information change in evolutionary optimization? With Figure 10 and Equation (16) in mind we recognize that sequence length and coding capacity have no direct influence. Intuitively we might expect that more coding capacity will allow for the construction of more efficient enzymes and thereby increase fitness. Indeed examples for a positive correlation of genome size with reproductive success were found in nature (for a recent example see [99]), but there is no correlation in general. Cases where the genome length decreased with increasing fitness are also well known. The most popular examples are genome size reduction in experimental evolution found with bacteria [100] and the spectacular loss of genes in the RNA of the bacteriophage Qβ in the first experiment on evolution in the test-tube by Sol Spiegelman and his group [2,13].

Because of the very high percentage of non-protein-coding DNA in higher organism the genome length is no measure for the number of genes. This is not so in bacteria with rather small fractions of non-coding DNA where the number of genes correlates strongly with genome length [101]. On a rough scale, however, genome size in biological evolution increases with the complexity of organisms or, perhaps expressed better, there are no complex organisms with genomes as small as those of bacteria [7]. Apparently there were evolutionary events where Shannon information or coding capacity underwent substantial increase. It is suggestive to identify such events with major transitions where several carriers of genetic information come together to form a cooperative unit.

Symbiosis is indeed one straightforward way among others to increase substantially the amount of genetic information: The genomes of two or more organisms become available for the superorganism and the same is true for semantic information. Eukaryotic cells formed by endosymbiosis are presumably the best understood examples: The nucleus and the organelles share their genetic information but processing of information reveals hierarchical control. In the toy model discussed here the Shannon information is—redundancy neglected—essentially doubled or tripled and the semantic information receives new qualities through the evaluation of skills and abilities of the functional organization leading to cooperation. One important prerequisite for the transitions can be seen already in the simple model: The basic molecular function, which has the capacity to induce cooperation, needs to be present prior to the transition. In our case this is the capability of ${\mathrm{X}}_j$ to perform second order catalysis in the reaction $\mathrm{A}+{\mathrm{X}}_i+{\mathrm{X}}_{j}0.25ex\begin{array}{c}\hfill c{l}_{ij}\\ \hfill \to \\ \hfill \end{array}2{\mathrm{X}}_i+{\mathrm{X}}_j$. This reaction must exist in the repertoire of possible catalytic processes but it plays no major role before the transition, since $a_0$ is too small under conditions leading to competition. To give an example, which is often mentioned in the context of prebiotic evolution: RNA is not only capable of being a template for replication it is also a catalyst and has the capacity to form catalytic networks and to act as universal replicase [20,21]. One caveat must eventually be repeated: What the toy model sketches is—at best—the origin and the beginning of a major transition and what we observe at present is the result of a long lasting evolutionary process with a plethora of steps most of them optimizing cellular performance. The organization of a eukaryotic cell, for example, is entirely different and much more complex than that of a prokaryotic cell. The change in cellular organization has not come about in a single step leading to endosymbiosis.

In order to stress the existence of a diversity of mechanisms leading to increases in information we mention an alternative process that can lead to substantial gains in syntactic and semantic information. Gene duplications are important evolutionary events, which increase the functional repertoire of organisms [102,103,104,105]. The initiation of gene duplication increases the genome lengths but neither the syntactic nor the semantic information because in essence only redundancy in the genome had been increased. Then the function of one copy of the duplicated genes is altered by mutation and this creates new function. In fortunate cases the new function is integrated into the reaction network of the organism and the new gene is stably and permanently integrated. Whole genome duplication is a rare event but happens and has a decisive influence on further evolution. The best studied example is a genome duplication in yeast [106,107]. Out of 16 genes in one particular segment 14 genes are eliminated and only two genes stay integrated and are ready for adaptation to new functions.

9. Discussion

The analysis and the discussion of the results were focussing on the dynamics taking place on the three faces of the Cartesian evolution space shown in Figure 1. Although it is not difficult to write down the general equations as we did in (1) and (4), analytical results for the dynamics in the interior of the evolution space are very hard to derive. On the other hand solution curves are easily obtained by numerical integration or stochastic simulation. Near the three surfaces the solution curves commonly resemble those discussed for the cases of the third parameter being zero. Sufficiently far away from all three surfaces it is not risky to conjecture that an asymptotically stable state with no extinct species in the sense of a quasispecies somewhat distorted by cooperation between the partners will be found. Depending on the intensity of the cooperation parameter and the number of species the trajectories will either converge monotonously to the stable state or spiral into it. The threshold phenomena are likely to become smoother and will eventually disappear inside the evolution space.

Despite its simplicity the model presented here is perfect for the description of in vitro evolution of RNA [108]. Needless to say there is ample room for making the model, in particular the reactions more complicated. Template induced RNA replication by means of virus specific replicases, for example, follows a complex many step polymerization mechanism [14,15,109], which under suitable conditions can be approximated properly by the simple autocatalytic reaction $\mathrm{A}+\mathrm{X}\to 2\mathrm{X}$. It is not difficult to introduce more elaborate many step mechanisms into the kinetic equations but the attractive possibility to do analytical mathematics is lost. The same is true for the embedding of the reactions into a flow reactor: More complex environments are readily conceived and implemented but the analytical approach becomes very tedious or even impossible for practical purposes.

What has been sketched here is, in essence, the evolution of populations by competition, cooperation, and mutation, as it occurs in laboratory assays under controlled experimental conditions [3] or with viruses in nature. Bacteria follow similar regularities in experimental evolution [11,12,110]. Bacteria in nature and eukaryotes have a much richer repertoire of variation including occasional and regular recombination of genetic information and require more involved models. In particular, the simultaneous consideration of mutation and recombination is anything but trivial.

Major transitions in evolution are, of course, much more complex than any simple model can describe. Needless to say, realistic and detailed molecular models for transitions are sophisticated and involve high-dimensional functional, metabolic and genetic networks. The formation of a cooperative collective of otherwise competing entities is only the first step of the transition and as said before further logical steps creating a new unit of selection at the higher hierarchical level must inevitably follow [6,111,112]. An important issue in most major transitions is the loss of autonomous replication capabilities of the individual elements in favor of the reproduction of the entire new unit. Another basic issue is the avoidance of exploitation by nonmembers of the new organization. Delimitation is required and this can be achieved spatially through cell membranes, cell walls, skins, shells, or other structures forming borderlines. An indispensable requirement for growth in compartments is coordinated synthesis of genetic material and cell division. Separation of societies occurs through signals or languages that are only understood by the members. The delimitation of the new system from the rest of the world leads to the new unit of selection and reintroduces the Darwinian optimization mechanism among variants at the higher hierarchical level. There are other features of major transitions but they are not universal in the sense that they occur in all transitions. The transition from independent to coordinated replication is mediated through catalyzed replication: The reproduction of a replicator is supported by the presence of another replicator. In case replicators are polynucleotides, RNA or DNA, catalysis in the form of ribozymes or deoxyribozymes is well known and discussed as an important feature of the hypothetical RNA world [113]. Catalysis by acting on polynucleotides has one very attractive feature: The basis for the catalysis is given by the molecular structure whereas the specificity can be tuned by sequence complementarity. The catalytic terms lead to suppression of competition and to the formation of groups of cooperating replicators provided sufficient resources are available. A new form of organization integrating several replicators into a functional unit is created and after other important steps stabilizing the new unit, the transition is completed.

Finally, we come to the role of information in the model presented here. We consider both Shannon information and semantic information in their simplest form: Shannon information as the coding capacity of an information carrier and sematic information as fitness derived through evaluation of the phenotype associated with the information carrier (Figure 10), and distinguish (i) optimization through mutation (Figure 2) and selection; and (ii) the transition from competition to cooperation. Optimization based on pure selection is not reflected by any systematic change in Shannon information. In simple systems semantic information deals with the evaluation of mutants with respect to fitness and the selection Equation (9d) guarantees increase—or more precisely no decrease—of semantic information in evolution. Mutation complicates the picture through formation of mutants clouds. Quasi-stationary populations in the sense of quasispecies, however, follow a similar scenario since a new quasispecies can replace a previous one only if the mean population fitness is larger.

The role of genetic information in major transitions is much more involved, and we observe different processes in different transitions. In the case of symbiosis—the example upon which we were focussing here—the symbiontic unit contains the genetic information of all partners at the beginning, and during the consolidation period various processes may take place. Genes no more needed by all partners can be eliminated from the genomes of some subspecies or genes may be transferred from one subspecies to another as it has happened, for example, in eukaryotic endosymbiosis in from of gene transmission from the mitochondria to the nucleus. The crucial steps required for reproduction at the higher hierarchical need to come under central control. In the eukaryotic cell this is control by the nucleus. Shannon information content, in general, will be larger in the symbiontic unit than in an individual subspecies but in total the information content will decrease because redundances can and will be eliminated. Evaluation in the form of semantics was important before the transition and will occur again at the higher hierarchical level between the new units of selection which evolve by optimization through mutation and selection. Little more can be said in general since individual major transitions are very different in detail [6]. One feature, nevertheless, is worth being considered. The property that is operative for the coupling of subspecies must already exist before the transition. To give an example coupling competing RNA replicators to a cooperative unit requires a molecular reaction, which can suppress Darwin’s natural selection providing a specific bonus for replication in presence of a partner molecule. RNA catalysis of RNA replication is the required capacity for the transition from independent replicators to the replicating network.

Despite the fact that the model is exceedingly simple it makes three highly relevant predictions: (i) Cooperation between otherwise competing partners requires more resources than optimization of properties by means of a Darwinian mechanism; and (ii) stable symbiosis can be formed only by small numbers of replicating units. Indeed many examples of symbiontic units formed by two partners are known, three-way symbiosis is not uncommon, but four way symbiosis is very rare and higher forms of symbiontic organizations are unknown [27,114]; and (iii) mutation at sufficiently high rates can prevent oscillatory systems from extinction. The fact that effective cooperation occurs only above a certain limit in the abundance of resources is among other things a consequence of the molecularity of the two autocatalytic reactions: First order autocatalysis or template induced replication without catalysis by another replicator is, in essence, a second order process whereas the catalyzed process is third order in concentrations and hence dominates only at high concentration, which are tantamount to large resources.