Molecular evolution between chemistry and biology

Biological evolution is reduced to three fundamental processes in the spirit of a minimal model: (i) Competition caused by differential fitness, (ii) cooperation of competitors in the sense of symbiosis, and (iii) variation introduced by mutation understood as error-prone reproduction. The three combinations of two fundamental processes each, (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal A}$$\end{document}A) competition and mutation, (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal B}$$\end{document}B) cooperation and competition, and (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal C}$$\end{document}C) cooperation and mutation, are analyzed. Changes in population dynamics that are induced by bifurcations and threshold phenomena are discussed.


Introduction
Biology is centered around evolutionary thinking or understanding biology implies understanding evolution as Theodosius Dobzhansky pointed out clearly in his famous book on evolution: "Nothing in biology makes sense except in the light of evolution" (Dobzhansky et al. 1977). Pure Darwinian evolution is a simple process but its embedding in nature renders it complex: Natural selection would follow uncomplicated laws in a simple environment. In the light of current molecular biology, there is need for a simple but comprehensive mathematical model of evolution to be able to account for modern genetics. The various epigenetic mechanisms have to be part of any comprehensive model of evolution and to keep such a model amenable to analysis and handling, molecular details must be reduced to a coarse-grained level. This article deals with a flexible model of evolution under defined environmental conditions. We present a concise and comprehensive review of work that was published elsewhere (Schuster 2016a(Schuster , d, 2017a together with a few new results. The model focusses on three basic processes: (i) fitnessdriven competition through differential reproductive success, (ii) reproduction-relevant cooperation between competitors, and (iii) reproduction-induced variation. The model is conceived with a modular structure and allows for the implementation of different, more or less complicated mechanisms for the processes (i), (ii), and (iii). For example, variation may be implemented by mutation, by recombination or by both. Here, we shall apply the simplest conceivable mechanisms: reproduction as single enzyme mediated replication (Biebricher 1983), cooperation of competitors as hypercycle dynamics (Eigen and Schuster 1978a, b), and variation as single point mutations based on the uniform error rate assumption (Swetina 1982).
Usage of the notion evolution is often ambiguous and a precise definition is desirable. Here evolution is understood as a process based on reproduction of a genotype being a DNA or an RNA sequence that carries encoded information on the formation of a phenotype, which is evaluated with respect to success in reproduction. The evolutionary process is built upon two foundations: (i) the dynamics at the population level and (ii) the environment dependent encoding of phenotypes in genotype. At the molecular level, the latter boils down to sequence-structure-function relations . The simplest systems that are capable of evolution in the sense of the given definition are nothing but special autocatalytic reactions involving polynucleotides under suitable conditions. In the next "Preliminaries", some prerequisites are presented for the model, which is introduced in "A minimal mathematical model for the evolution of molecules". The model comes in a deterministic version based on kinetic differential equations or it is formulated as a stochastic process modeled by means of chemical master equations . Although almost no analytic solutions are available for master equations derived from nonlinear reaction kinetics 1 in two or more variables, the stochastic version of the model can be studied by efficient simulation methods for the calculations of trajectories (Gillespie 1977(Gillespie , 2007. In three sections, solution curves for the three two-dimensional subspaces, (  ) reproduction and mutation, (  ) reproduction & cooperation, and (  ) cooperation & mutation, are presented, analyzed and interpreted. Then follows a brief discussion of results obtained with the full three-dimensional model, reproduction & cooperation & mutation and in the final section we discuss the simple model in the context of the complex processes in nature.

Preliminaries
Darwinian evolution is often-and incorrectly-seen as a synonym for optimization mainly because Ronald Fisher's fundamental theorem of natural selection (Fisher 1930(Fisher , 1958Price 1972;Ewens and Lessard 2015) focusses on non-decreasing mean fitness (t) in evolving populations. In the simplest form derived from idealized models, the time derivative of the mean fitness is the variance of the fitness values and therefore a non-negative quantity. In reality, mean fitness is the result of many factors and as Fisher was been certainly aware, only a few cases-like single locus genetics-fulfil the theorem in pure form (for elaborate discussions see the more recent literature Plutynski 2006;Okasha 2008).
Evolution is intimately related to the environment in which it takes place and accordingly environment and environmental changes are major factors shaping evolutionary processes. Here, we are primarily interested in the internal dynamics and hence a well-defined and controllable environment is required. For this goal, we introduce a flow reactor in "Idealized environments", which represents a simple device that is not only useful for performing evolution experiments but also provides at the same time a suitable setup for theoretical modeling. More complex environments can be implemented as long as they can be cast in kinetic equations. In a separate "Basic processes", the processes that will be used to model evolution are introduced. The following section describes the deterministic and stochastic methods, which are applied to find solutions of the kinetic equations. Finally, we review some fundamental features of autocatalysis since this is the chemical counterpart of reproduction.

Idealized environments
Environments that allow for investigations of observations as functions of one or few parameters with everything else being constant require elaborate design in the form of sophisticated experiments since devices controlling environmental conditions may be quite involved. In theoretical approaches, often the silent assumption is made that there exists a hypothetical machinery, which takes care of fixing parameter values as needed for the mathematical analysis.
Environmental influences on phenotypes are commonly large, manifold, and easy to observe. In this study, however, we are interested in the intrinsic driving forces of evolution, which result from reproduction, symbiotic cooperation and variation, and therefore impacts on evolution caused by changes in the environment are intended to be kept to a minimum. To reduce the influence of the environment as much as possible we shall assume a control device in the form of a simple flow reactor ( Fig. 1; see also Schmidt 2005). More elaborate reactors, which keep, for example, the numbers of bacterial cells constant have been designed and implemented (Novick and Szillard 1950a, b;Bryson 1952). The flow reactors called chemostat, cellstat or turbidostat and other experimental devices for monitoring and controlling evolution in the laboratory may serve as examples (Husimi et al. 1982;Dykhuizen and Hartl 1983;Husimi 1989;Koltermann and Kettling 1997;Strunk and Ederhof 1997).
Implementation of a physical device rather than application of idealized assumptions like constant population size is required for the stochastic description of evolution. An illustrative example where the deterministic approach yields a stable solution whereas the corresponding stochastic system is unstable is provided by the linear birth-and-death process (Goel and Richter-Dyn 1974), which is described by the kinetic equations + f ⟶ 2 and d ⟶ ∅ . For equal birth and death rate parameters, f = d , the population number of the deterministic system stays constant whereas in the stochastic model the fluctuations are unregulated. With increasing amplitude of fluctuations, the populations will hit the death state, which is an absorbing boundary and where the system therefore remains caught forever. In other words, the system is unstable despite a (marginally) stable deterministic solution. The direct incorporation of constant population size into Fisher's selection (1930) or Eigen's equation (1971) leads to an instability of the same kind since fluctuations are not self-regulating (Jones and Leung 1981). Implementation of the system in the flow reactor provides stability due to balance control of inflow and outflow modeled by the pseudoreactions * a 0 ⋅r ⟶ and r ⟶ ∅ as well as i r ⟶ ∅ (i = 1, … , n) ( Fig. 1).

Basic processes
The minimal system for modeling evolution of molecules is based on the three classes of processes: (i) competitive reproduction, (ii) symbiotic cooperation, and (iii) reproduction based variation. For the minimal model we shall choose the simplest possible chemical reactions. Reproduction will be modeled by simple replication, + i f i ⟶ 2 i with f i being the fitness of species i and competition being introduced through living on the same resource A. Symbiontic cooperation is introduced as catalyzed replication, , where i represents the template and j is the catalyst. In its most general form-every molecule j has the potential to act as catalyst in the replication of every molecular species i -the number of catalytic terms is n 2 and unrealistically large since specific catalysis is a rare property. The simplest example of stable cooperative catalytic networks with fewer catalytic reactions known so far is the catalytic hypercycle (Eigen 1971;Eigen and Schuster 1977: The catalyzed reactions form a closed cycle with n members. Genetic variation occurs at the level of a DNA or RNA genotype in forms of mutation and recombination. The simplest form of variation is the point mutation that consists of the exchange of a single nucleotide in the sequence and caused by the incorporation of a wrong nucleotide during the replication process. Correct and error-prone replication are considered as parallel reaction channels within one and the same replication mechanism (Eigen 1971). In terms of a simple over-all replication kinetics of the multistep process, with E being a replicase enzyme, the reaction rate is obtained as the product of two parameters: the fitness of the template that depends on the availability of resources here the concentration [ A]. The dimensionless factors Q ji with i, j = 1, … , n are understood as the elements of a mutation matrix Q = {Q ij } and provide the probability to obtain j as a copy of the template i that can be either correct ( j ≡ i ) or error-prone ( j , j ≠ i ). Conservation of probabilities requires: ∑ n j=1 Q ji = 1 since each copy has to be either correct, Q ii , or incorrect, ∑ n j=1,j≠i Q ji = 1 − Q ii . For the sake of simplicity, binary sequences will be considered here and we distinguish only between correct and incorrect base pairs.
A useful simplifying approximation is made by the uniform error rate model (Swetina 1982): The error per nucleotide and replication event, p, is assumed to be independent of the position of the nucleotide along the sequence and the nature of the nucleotide to be complemented. Then all elements of the mutation matrix can be expressed by a simple formula, with only three parameters: (i) the sequence length of the RNA molecules l, (ii) the mutation rate p, and (iii) the Hamming distance (Hamming 1950(Hamming , 1986 between the two sequences interrelated by the mutation process, d ij = d H (X i , X j ) . Without changing important results for the purposes pursued here, the analysis of the model is substantially simplified by the assumption of constant chain lengths l, which is also consistent with the restriction to As an alternative to Eq. (1), mutation can be seen as a consequence of DNA change-damage and imperfect repair-during the whole life time of an organism, which is the idea in the Crow-Kimura mutation model (2009), pp 264-266]: Interestingly, both models (1) and (3) although different with respect to the underlying physics give rise to identical mathematical problems (see e.g. Baake andGabriel 1999 andSchuster 2016, pp 76-78).

Deterministic and stochastic approaches
Reaction mechanisms are commonly analyzed by deterministic and stochastic approaches. The former translate the chemical reaction equations into kinetic differential equations, which can be directly solved by mathematics, studied by means of qualitative analysis or investigated by numerical integration. The theorems of existence and uniqueness of solutions of differential equations are applicable and a single integration provides the complete information for a given input set. Competitive selection with nonzero mutation leads to one unique asymptotically stable stationary state (Thompson 1974;Jones et al. 1976), whereas the long-time dynamics of cooperative systems is much richer and multiple stationary states, oscillations or deterministic chaos may be observed Schnabl et al. 1991).
Stochastic analysis in general in based on searching for a stochastic process that fits the model to be studied as closely as possible. Chemical reaction kinetics prefers master equations although the repertoire of analytical solutions is very limited. It is not difficult to write down a multivariate master equation but the derivation of analytical solutions is successful only in exceptional cases, for example for networks of monomolecular reactions (Jahnke et al. 2007;Deuflhard et al. 2008. If no analytical solutions are available information on the stochastic system can be obtained by trajectory sampling. The theoretical background for trajectory harvesting has been laid down by Andrey Kolmogorov (1931), Willy Feller (1940, and Joe Doob (1942Doob ( , 1945. With electronic computers now being generally available elaborate simulations of stochastic processes became possible. The more recent conception, analysis, and implementation of a simple but highly efficient algorithm by Daniel Gillespie (1976, 1977 provides a very useful tool for investigations of stochastic effects in chemical kinetics. Distributions of trajectories are characterized by expectation values and higher moments, commonly only by variances or standard deviations. (3) Equilibrium fluctuations in conventional chemical reaction kinetics follow an approximate √ N-law and hence play almost no role in systems with particle numbers that are typical for chemical systems. In biology a different scenario is encountered. For example, every new variant originating from mutation has to start out from a single copy. The deterministic approach commonly uses continuous variables, which can only be an approximation to reality, since numbers of molecules or biological entities are integers by definition. Continuous concentrations can adopt arbitrarily small values whereas stochastic variables cannot pass low values beyond unity, because then molecular species go extinct, and deterministic solutions become unrealistic and differ strongly from the stochastic results. Large population size alone is not sufficient, each variable has to be sufficiently large at every instant to guarantee similarity between stochastic and deterministic solutions. Two more phenomena are relevant in the stochastic approach at low particle numbers and may lead to large standard deviations in the variables: (i) nucleation steps in reactions involving two or more molecules and (ii) multiple (quasi)stationary states 2 of the stochastic system. The stochastic system may approach any stationary state and the distribution of the possible outcomes is determined by probabilities. Then the calculations of expectation values and higher distribution moments build averages over different final states and may give rise to enormous standard deviations.
For the purpose of illustration, we consider the equilibration of the flow reactor as an example of a stochastic system that can be fully analyzed by analytic calculations (Schuster 2016, pp 436-441). The two pseudoreactions, are converted into the easy to solve kinetic differential equation where n 0 = a(0) . The master equation for the probability P n (t) = P A(t) = n with number density A(t) being the discrete pendant of concentration a(t), P n t = r a 0 P n−1 + (n + 1)P n+1 − (a 0 + n)P n , n ∈ ℕ , 1 3 has the analytical solution with n, n 0 , a 0 ∈ ℕ for the sharp initial condition P n (0) = n,n 0 . In the limit t → ∞ the distribution (4e) converges to a Poisson distribution First and second moment yield expectation values and standard deviations The standard deviation starts out from (t = 0) = 0 and approaches the equilibrium value = √ a 0 either from below for a 0 > n 0 or from above for a 0 < n 0 after having passed a maximum at In general the maximum if it exists is very flat (see Fig. 3) as can be easily checked from the analytical expressions through inspection: The flow reactor is one of the few rather exceptional cases where the stochastic approach can be done completely by mathematical analysis.

Autocatalysis in the batch reactor
A batch reactor is an elaborate device that allows for performing chemical reactions under controlled conditions without inflow and outflow (Schmidt 2005). Here, the term batch reactor is used to indicate that reactions are carried out in a well-mixed closed system under constant temperature and pressure, and in the long run approach a thermodynamic equilibrium state.
In conventional chemistry autocatalysis is a rather rare phenomenon but in biology it represents the most important process since multiplication is just a special form of autocatalysis that in simplest form can be expressed by the reversible reaction (4e) The forward reaction of (5) leads to molecular self-enhancement: The resource A is consumed in the synthesis of the replicator X and the process is catalyzed by the presence of further molecules X. Depending on the molecularity of the autocatalytic process as expressed by the value of m autocatalysis comes in different forms that, in essence, fall into two classes: (i) simple or first-order autocatalysis with m = 1 and (ii) second-and higher-order autocatalysis with m ≥ 2 [ (Phillipson and Schuster 2009, pp 18-27)]. For consistency, we add here also the uncatalyzed reaction and call it zeroth order autocatalytic ( m = 0). 3 All processes in closed systems converge to an equilibrium state and show mass conservation that implies a conservation relation The reversible autocatalytic reactions converge to the equilibrium state: Although the expressions for the equilibria are the same for the all reactions independently of the stoichiometric coefficient m, K = g → ∕g ← , there are subtle differences in the probability distributions, which become important at small concentrations: The particle numbers are discrete quantities, change by integers only, and the stoichiometric factors are 3 rather than X(t) 2 and X(t) 3 , respectively. 4 The states with X(t) = 1 for m = 1 , or the states X(t) = 1 and X(t) = 2 for m = 2 cannot react because two or three molecules X are needed for the conversion of X into A. The inaccessibility of the state with X(t) = 0 or the states with X(t) = 0 and X(t) = 1 in the first-or second-order autocatalytic reactions require different normalizations for the stationary probabilities of the three systems: where as before the stochastic variable is n = A . Mass conservation provides X = N − n . As expected the truncation of P (0) n is important for small values of N only. For firstorder autocatalysis with N = 10 , k → = k ← = 1.0 we obtain E(A) = 4.9951 , for second-order autocatalysis E(A) = 4.9605 compared to E(A) = 5 for the uncatalyzed process.
In principle, there are three major sources of randomness in autocatalytic reactions: (i) thermal fluctuations, (ii) delayed onset of autocatalytic reactions and (iii) multiple stationary states. (ad i) In the transients towards equilibrium the thermal fluctuation bands are essentially the same in autocatalytic and conventional reactions as can be seen best from the comparison of standard deviations at equilibrium (Fig. 2). The differences between conventional and autocatalytic equilibrium densities can be recognized numerically for very small particle numbers only (). (ad ii) The reaction rate for an autocatalytic reaction of order m is v(0) = k A(0) X(0) m − l X(0) m+1 (with m ≥ 1 ). At small t the factor A(t) is large and the factor(s) containing X(t) are small and only the first term in v(t) matters in the early phase of the reaction. Thus X is produced and X(t) increases, which in turn leads to an increase v(t) yielding the typical scenario of self-enhancement. Selfenhancement in chemical reactions is tantamount to an increase of the reaction rate with concentration in the early phase and together with the late saturation phase gives rise to "s"-shaped or sigmoid curves whereas the uncatalyzed reaction ( m = 0 ) follows a simple exponential decay (Fig. 2). Higher values of m lead to steeper curves, which approach a step function with increasing m. The maximum standard deviation in the approach towards equilibrium max = max{ (t)} is a measure of the random scatter in the delay in the onset of the autocatalytic reaction (Table 1). (ad iii) Multiple final states give rise to an additional stochastic component often called anomalous fluctuations (de Pasquale et al. 1980) since "Autocatalysis in theow reactor").
In Fig. 2 transients for the two processes + m ⇌ (m + 1) with m = 0 and 1 are compared by means of expectation values and fluctuation bands. As initial conditions we apply an empty reactor, A(0) = 0 , and the smallest possible values for X: X(0) = 0 and X(0) = 1 for the uncatalyzed and the autocatalytic reaction, respectively. A total population size of N = 1000 was chosen so that the one-standard-deviation fluctuation band of order √ N appears small and the deterministic solutions coincide with the expectation values in the reference process A ⇌ X ( m = 0 ). The transient of the autocatalytic process ( m = 1 ) is different: It shows substantial broadening of the fluctuation band (Table 1) because of delayed onset of the reaction before it narrows down to the equilibrium value. In the intermediate range the expectation values differ remarkably from the deterministic solution curves. As expected an increase in X 0 leads to a decrease in the width of the fluctuation band. Second order autocatalysis ( m = 2 ) differs from firstorder autocatalysis mainly by the size of the characteristic effects: Both autocatalytic reactions show broadening of the fluctuation bands but the band width in the second-order case is about four times as wide. In particular, the scatter in the waiting times until the first reaction events occurs is much larger because we are dealing with two small factors, The standard deviation in the course of the reactions, (t) , is shown in Fig. 3. Because of sharp initial conditions the fluctuation band starts out from zero-(0) = 0 , increases, becomes broad in the intermediate range and settles down at the equilibrium value (6). Substantial deviations between the deterministic solution and the stochastic expectation value, a(t) and the E A(t) or x(t) and E X(t) , respectively, are observed in the intermediate range. Accordingly, the standard deviation goes through a pronounced maximum that is (black) and E X(t) (red) together with the one standard deviation band E(t) ± (t) (gray and pink) obtained by sampling of 10,000 trajectories that were calculated by Gillespie's simulation method for the uncatalyzed ( m = 0 ; top plot) and the first-order autocatalytic reaction ( m = 1 ; bottom plot). Both reactions approach almost identical thermodynamic equilibrium states. Choice of parameters and initial conditions: . Solutions of the corresponding kinetic differential equations are shown as broken lines.
qualitatively different from the shallow maximum observed with the standard deviation of conventional chemical processes see Eq. (4h) .
The ir reversible processes, + m → (m + 1) obtained from Eq. (5) by putting g ← = 0 , are illustrative because they are closer to biology where replication or reproduction occur always under the conditions of irreversibility. The shape of the solution curves compared to those shown in Fig. 2 shows similarities except a twice as wide range for the values of the stochastic variables, (0) and expectation values and standard deviations approach zero at sufficiently long times. Again, we observe a sigmoid shape of the solution curves, for small initial values ( X(0) < 5 ) the standard deviation (t) becomes large in the intermediate range (Fig. 3), and the deterministic curve deviates substantially from the stochastic expectation values.

Autocatalysis in the flow reactor
Implementation of autocatalytic reactions in the flow reactor provides additional insights into the different forms of randomness. In particular we are interested in multiple stationary states as a source of stochasticity (item iii). The reaction equations for first order autocatalysis are: with the kinetic differential equations The simple first-order autocatalytic process in the flow reactor sustains two long time states: (i) The state of extinction S 0 with lim t→∞ a(t) = a = a 0 and lim t→∞ x(t) = x = 0 , and (ii) the quasistationary state S 1 with lim t→∞ a(t) = a = r∕k and

Fig. 3
Comparison of the standard deviation in the reversible first-order autocatalytic reaction + ⇌ 2 and the reversible uncatalyzed reaction ⇌ . The reactions are recorded for the closed system fulfilling The figure presents the results of statistical evaluation of 10 000 trajectories obtained by computer simulations with Gillespie's method (Gillespie 2007) for the autocatalytic reaction (1) X (t) = (1) (t) (black) and for the uncatalyzed reaction (0) The equilibrium value of the standard deviation is practically the same for both reactions: Table 1 Fluctuations in the two autocatalytic reactions A+ X ⇌ 2 X and A+2 X ⇌ 3 X a The table presents maximal standard deviations max computed from 1 000 trajectories of the two autocatalytic reactions for different initial conditions X 0 = X(0) . For the autocatalytic reactions the standard deviation (t) passes through the maximum max listed here whereas for the uncatalyzed process it increases monotonously from (0) = 0 to the equilibrium value (see Fig. 3

) a
The equilibrium fluctuations calculated from equations () are practically the same for all three reactions. Choice of parameters, N = 1000 , A ⇌ X: Depending on rate parameters and initial conditions the trajectories may pass a flat maximum before they decrease to the equilibrium value (see Eq. (4h) and Schuster 2016, pp. 445-449) b The accurate value obtained from the stationary master equation is = 15.8114

Reaction
Initial conditions Limit Starting from an empty reactor containing no A and the autocatalyst only in seeding quantities, A(0) = 0 and X(0) = 1, 2, 3, … , the trajectory shown in the upper part of Fig. 4 allows for the distinction of several phases: (I) the flow reactor is filled with the resource A in phase I, (II) in the random phase II the decision is made to which state the trajectory will converge, (III) the trajectory approaches the long-time state, and (IV) the trajectory fluctuates around the state (see Figs. 4 and 7). First-order autocatalysis sustains the two long-time solutions S 0 = (a (0) , 0) and S 1 = (a (1) , x (1) ) , stochastic trajectories approach either of the two states, and parameters and initial conditions determine the probabilities to end up here or there. For sufficiently large population sizes, the long-time expectation values of the stochastic variables can be well approximated by linear combination of the deterministic values: and N 1 are the counted numbers of trajectories ending up in S 0 or S 1 from a sufficiently large sample. Although only 3/100 trajectories go to extinction in the example shown in the lower plot of Fig. 4 the influence on the expectation values E A(t) and E X(t) and the standard deviations A(t) and X(t) is remarkable. This random component of processes has been intensively studied in the nineteen eighties by Paolo Tombesi, Francesco de Pasquale, Piero Tartaglia and the notion anomalous fluctuation caused by a chemical instability was coined for this kind of stochasticity (de Pasquale et al. 1980. It is advantageous to collect trajectories separately for the different final states, because then the anomalous fluctuations disappear. For example, the standard deviation of A at t = 30 is reduced from A(30) = 335.7 to A(30) = 7.12 if one changes from a sample of hundred trajectories with three extinction events (Fig. 4, Pseudorandom number generator: Extend-edCA, Mathematica, seeds s = 491 ) to one without extinction ( s = 919).
The major difference between the two classes of autocatalytic reactions lies in the repertoire of possible dynamic behaviors. First-order autocatalysis gives rise to exponential growth in unconstrained systems and to selection and optimization of mean fitness in multispecies cases with finite resources (see ''Competition, mutation and quasispecies"). Accordingly first-order autocatalysis leads to a Darwinian scenario of selection of the fittest. In contrast to first-order autocatalytic reaction networks, the dynamics in second-order systems is very rich and includes multiple stationary states, oscillations and deterministic chaos. The second-order autocatalytic elementary step, A+2 X ⇌ 3 X, represents a kind of generally used prototype for theoretical models, for example the Brusselator (Nicolis and Prigogine 1977). It provides a simple enough reaction step for studies by means of rigorous mathematics. Qualitative analysis of Brusselator dynamics is straightforward and numerical integration causes no problem provided the integration software can handle stiff differential equations. In reality, however, single-step autocatalytic reactions are extremely rare, instead we are commonly dealing with multistep-reaction networks (Noyes et al. 1972) (see also the review by Francesc Sagués and Irving Epstein (2003).
In biology, in particular in the theory of evolution, the process A+2 X → 3 X plays a special role since in the simple form of hypercycles it is the basis for suppression of The sample size of the bottom plot was 100 trajectories, 3 led to S 0 (extinction) and 97 reached the pseudo-stationary state S 1 . Solutions of the corresponding kinetic differential equations are shown as dashed lines competitive selection without destroying template-induced reproduction. It provides one fundamental mechanism for major transitions (Maynard Schuster 1996) and will be discussed extensively in "Competition and cooperation". The enormous flexibility of secondorder autocatalysis follows, for example, from Fisher's selection equation and the proof for the optimization of mean fitness in sexual reproduction under idealized condition. In a caricature model we may explain how the above mentioned reaction step could be related to sexual reproduction: 2 X on the l.h.s. are (at least stoichiometrically) related to the parents and 3 X on the r.h.s. model parents and one offspring. Apart from being illustrative toys autocatalytic processes set the stage for modeling evolution in the sense that we shall reencounter all special phenomena of autocatalysis in the more elaborate model for evolution to be presented and discussed in the next "A minimal mathematical model for the evolution of molecules".

A minimal mathematical model for the evolution of molecules
The minimal model is dealing with the time dependence of the distribution of genotypes in populations Π(t) and hence it is rooted in chemical kinetics of (i) competitive reproduction, (ii) symbiontic cooperation, and (iii) genetic variation. The quantification of these three properties yields three parameters or sets of parameters, which can be plotted on the three axes of a Cartesian coordinate system (see Fig. 5 and the previous publications (Schuster 2016a(Schuster , b, d, 2017b. We consider the simplest case here, where the three quantities are a fitness parameter f, a cooperation parameter h, and an mutation rate parameter p. For an implementation of the model in the flow reactor we need two additional external parameters measuring the accessible resources expressed, for example, as number density or concentration of a (hypothetical) compound A, A 0 or a 0 , respectively, and the mean residence time R = V∕r of the reaction mixture in the reactor where V is the volume of the reactor and r is the (volumetric) flow rate. The parameter R defines the time resolution of the reactor since slow reactions, which do not progress appreciably during the time interval R cannot be studied.
In the next step, the model is implemented by means of a suitable and simple reaction mechanism. Based on "Idealized environments" and "Basic processes" we consider the following set of 2n 2 + n + 2 reactions in the flow reactor: The process (8a) 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 parameter r (Schmidt 2005, p. 87ff). The reactor operates at constant volume and this implies that the volume of solution flowing into the reactor per time unit [t] is compensated exactly by an outflow, which is described by the Eqs. (8d) and (8e) and concerns all (n + 1) molecular species, A and i , i = 1, … , n . Inflow and outflow are often characterized (8c)

Fig. 5
A minimal model for modeling evolution. Evolution is considered as an interplay of three processes: (i) competition through reproduction, (ii) cooperation through symbiosis, and (iii) mutation through error-prone replication. In parameter space, the intensity parameters of all three processes, (i) fitness parameters f corresponding to reaction rates for competition, (ii) reaction rates h for catalyzed reproduction, and (iii) an error rate parameter p for mutation are plotted on the axes of a Cartesian coordinate system. On the three twodimensional faces of the coordinate system we are dealing with the three fundamental evolutionary processes: (  ) competitive reproduction and mutation are the basis of Darwinian optimization through natural selection and give rise to the formation of quasispecies and eventually to the occurrence of error thresholds (Eigen 1971;Eigen and Schuster 1977;  ) the interplay of competition and cooperation allows for the description of major transitions, which are seen as the consequences of crossing resource thresholds (Maynard Szathmáry and Maynard Smith 1995;Schuster 1996; (  ) the combination of cooperation and mutation enables reintroduction of extinguished species provided the error rate is sufficiently large such that a mutation threshold for stochastic survival can be recognized (Schuster 20160 1 3 as pseudoreactions because they are no chemical reactions in the strict sense, which are converting reactants are into products. The two classes of reactions producing progeny, template induced replication (8b) and catalyzed template induced replication (8c), represent the core of the evolution model. In agreement with the conditions in biology, both reproduction steps are implemented irreversibly in the direction of polynucleotide synthesis. A basic assumption for both reproduction steps is that correct reproduction and mutation are parallel chemical reaction channels ("Basic processes"). In other words, there is no mutation under conditions that do not sustain reproduction.
As an alternative to the Eigen model (1) mutation can be seen, for example, as the result of DNA damage and imperfect damage repair during the whole life span of an organism, which is the idea underlying the Crow-Kimura mutation model (Crow and Kimura 2009, pp 264-266). Then reproduction and mutation are completely independent processes, and in the kinetic differential equations they appear as additive terms. Interestingly, the Eigen and the Crow-Kimura model although being fundamentally different with respect to the assumptions about the nature of the mutation process give rise to the same mathematical problems [see, e.g., (Schuster 2016d, pp.76-78)]. The mutation matrix Q corresponds formally to the mutation matrix but there are non-negligible differences Q covers correct and error-prone replication but the process i → 2 i is handled separately in the Crow-Kimura model and hence all diagonal terms are zero μ ii = 0 ∀ i = 1, … , n.
The equations that will be applied in the analysis of the dynamics of the model (8) implement three processes along the coordinate axes: (i) Darwinian selection of the fittest on the competition axis, (ii) hypercycle dynamics on the cooperation axis, and (iii) neutral evolution on the mutation axis. The kinetic differential equations of the model mechanism (8) are of the form: In (9a) we made use of the conservation relation ∑ n i=1 Q ij = 1 . No analytical solutions are available for (9) in general but numerical integration is straightforward as long as n is not too large. In absence of cooperation terms, l i = 0 ∀ i = 1, … , n , Eq. (9) can be transformed into an eigenvalue problem of a symmetric matrix, which is readily diagonalized provided n is not very large ( n < 10 6 ; "Competition, mutation and quasispecies"). 5 In mutation-free systems, p = 0 ("Competition and cooperation"), qualitative analysis and determination of stationary states are straightforward, and the dynamics of the complete system can be derived by extrapolation from the error-free results to finite mutations rates. The cooperation system with mutation, (9) with k i = 0 ∀ i = 1, … , n is used here to study the relevance of mutation in symbiontic systems. It also serves as an example for the study of unconventional consequences of replication with frequent errors in the strong mutation scenario ("Cooperation and mutation").

Processes on individual coordinate axes
The processes along the individual coordinate axes are considered in order to verify the initial statement on the three basic processes. For this goal it is easiest to set certain parameters zero: (I) no natural selection implies k 1 = k 2 = … = k n = k (= 0) , (II) no cooperative coupling requires l 1 = l 2 = … = l n = l = 0 , and (III) no mutation leads to Q = where is the unit matrix. A process taking place on the selection axis is given by (II) and (III) being true leads to the ODE Competition between the reproducing elements i leads to survival of the fittest subspecies m , which is defined by k m a = f m = max{f i ;i = 1, … , n} . For constant [ ] = a 0 = const , the mean fitness of the population, , is a nondecreasing function of time: d ∕ dt = var{f } ≥ 0 . This result is the formal analogue of Fisher's fundamental theorem for asexual reproduction.
We remark that the deterministic kinetic Eqs. for (8b) and (8c) were extensively studied under the simplifying assumption of constant population size ∑ n i=1 x i (t) = c 0 = const (Eigen 1971;Eigen and McCaskill 1989). The deterministic solution curves formulated in relative concentrations are identical for the CSTR and for constant population size (Schuster and Sigmund 1985). As we mentioned already in "Idealized environments" the stochastic system is unstable for unregulated constant population size (Jones and Leung 1981) and a proper description has to account explicitly for the physical setup applied, here the flow reactor. assumption that chemical processes occur one at a time, all jumps involve single steps and the particle numbers change by ±1 unless the elementary steps involve more than a single molecule per class. The jumps S ′ → S or S → S ′ are denoted by the shorthand notation Then the master equation of mechanism (8) takes on the form Each reaction step involving S changes the probability to be in state S , P , in two ways: It increases the probability through reactions or pseudoreactions S ′ → S and decreases the probability through the reaction steps S → S ′ where ′ is summed over all states from which S can be reached or vice versa. The two terms in the first line, for example, describe the two pseudoreactions modeling inflow and outflow of the material A, and further each reaction is represented by two steps. It is also worth noticing that stoichiometry requires two slightly different replication terms depending on whether the copy is correct or incorrect.
Master equations are easily written down and stationary solutions can be derived by generally applicable techniques as it was shown for the flow equilibrium in "Deterministic and stochastic approaches" but explicit time-dependent solutions are very hard to obtain and known only in exceptionally simple cases [Schuster 2016e, pp 347-568]. Here, we shall use the simulation technique of sampling trajectories introduced already in "Autocatalysis in the batch reactor". Expectation values and second moments of variables can be computed through sampling of trajectories-an example is shown in Fig. 6-but often this approach exceeds the available computational facilities and it is necessary to interpret single trajectories. As examples we consider single trajectories in Fig. 4 and in Fig. 7. The process for convenience starting from an empty flow reactor is split into four phases: (i) establishment of the flow equilibrium of A, (ii) random decision on the (quasi)stationary state towards which the trajectory � = (m ± 1, s 1 , … , s n ) ≡ ( ; m ± 1) or � = (m, s 1 , … , s k ± 1, … , s n ) ≡ ( ; s k ± 1).
On the cooperation axis, the conditions (I) and (III) are fulfilled and we obtain the equations of hypercycle dynamics which were studied in great detail in a number of previous papers Schuster et al. 1979Schuster et al. , 1980Hofbauer et al. 1991).
The third case-with (I) and (II) being true-yields a degenerate deterministic solution: All distributions of subspecies with fixed population size ∑ n i=1 x i = c yield the same solutions and therefore constitute an (n − 1)-dimensional invariant manifold fulfilling In the absence of selection and cooperation neutral evolution in the sense of Motoo Kimura ( , 1983) is observed on the mutation axis. For an adequate description of the process a stochastic treatment is required. Random selection of a single arbitrarily chosen subspecies is observed in the no-mutation limit, lim p → 0 . For non-vanishing mutation rates once selected subspecies may be replaced by other subspecies and the mean time of replacement of one randomly selected subspecies is T rep = −1 where is the mutation rate per generation (Kimura , 1983. Translated into our model, we find for single point mutations ≈ p.

Master equation and simulation
Reaction (8) can be cast into chemical master equations. The particle numbers of the molecular species, [ ] = A(t) and [ i ] = X i (t) with i = 1, … , n , are integers and in the absence of flows they fulfil the conservation relations,

Competition, mutation and quasispecies
The bottom face of the three-dimensional parameter space- in Fig. 5-is dealing with selection provided by the combination of competition and mutation. Natural selection at zero mutation rate leads to a homogeneous population of the fittest subspecies and at non-zero mutation rates this scenario yields selection of a fittest ensemble of subspecies, which has been characterized as quasispecies (Eigen and Schuster 1977;Domingo and Schuster 2016). Precisely, the quasispecies is the stable stationary long-time distribution of a population of subspecies undergoing replication and mutation: A population that consists of several genotypes present in time-dependent concentrations, h a s t h e s t a t i o n a r y s o l u t i o n , lim t→∞ Π(t) = = x 1 1 ⊕ x 2 2 ⊕ … ⊕ x n n , which is called the quasispecies .
Deterministic quasispecies The deterministic or continuous quasispecies represents the unique deterministic longtime solution of the replication mutation problem, which is described in the flow reactor by the ODE (Schuster and Sigmund 1985): . 6 Quasispecies formation. The two plots show the formation of quasispecies from an initially empty reactor, A(0) = 0 , with replicators in seeding amounts. The system in the upper plot was initiated to be a single copy of the master sequence, X 1 (0) = 1 and X 2 (0) = X 3 (0) = X 4 (0) = 0 . At the beginning the system might be extinguished by a single dilution event X 1 → ∅ and the high probability of extinction gives rise to enormously broad and overlapping one-standard-deviation bands. The initial values in the lower plot were chosen such that the probability of extinction is zero for all practical purposes: X 1 (0) = 10 and X 2 (0) = X 3 (0) = X 4 (0) = 0 . Choice of parameters: , Color code: A black, 1 red, 2 yellow, 3 green, and 4 blue. Expectation values are shown as full lines, deterministic solutions as broken lines. Since x 2 (t) = x 3 (t) is fulfilled for the parameter values used here the curve is shown as a yellow-green broken line Fig. 7 Sequence of phases in the approach towards a quasistationary state for n = 2 . A stochastic trajectory simulating competition and cooperation ("Competition and cooperation") of two species in the flow reactor is shown in the plot above. The corresponding master equation is derived from (14) by putting Q = . The stochastic process is assumed to start with an empty reactor except seeds for the two autocatalysts 1 and 2 . It can be partitioned into four phases: (I) fast raise in the concentration of A, (II) a random phase where the decision is made into which final state-S 0 , S (1) 1 , S (2) 1 or S 2 -the trajectory progresses, (III) the approach towards the final statehere S 2 -and (IV) fluctuations around the values of the (quasi) stationary state. Comparison with the simple autocatalytic process in Fig. 4 reveals great similarity. Choice of parameter values: , pseudorandom number generator: Extend-edCA, Mathematica, seeds s = 631 . Initial conditions: A(0) = 0 , X 1 (0) = X 2 (0) = 1 . Color code: A(t) black, X 1 (t) red, and X 2 (t) yellow. The figure is adapted from Schuster (2017a) 1 3 All early works on quasispecies dynamics were performed with the constraint of constant population size: ∑ n i=1 x i (t) = c = const with a(t) = a 0 and f i = k i a 0 that gives rise to the differential equation (16) and (16') have identical solutions in normalized variables but the stability properties of the corresponding stochastic systems are different. [For more details see (Thompson 1974;Jones et al. 1976;Ebeling and Mahnke 1979;Jones and Leung 1981;Schuster and Sigmund 1985)].
The quasispecies as a function of the mutation rate parameter (p) begins at p = 0 as a homogeneous population containing exclusively the fittest subspecies m and becomes a distribution of subspecies at nonzero p. This distribution consists of a most frequent sequence called master sequence, which is surrounded by a mutant cloud. Commonly, the most frequent or master sequence is also the fittest one, but this is not necessarily so: In case of strong stationary mutation flow from the mutant cloud to the master sequence, a less fit master may outgrow a fitter sequence with less efficient mutational backflow. At further increase in p the distribution broadens, and eventually ends up as the uniform distribution  ∶ {x 1 =x 2 = ⋯ =x n } at p = 1 − ( − 1)p∕ = 1∕ where n = l and x i = 1∕ l (i = 1, … l ) for sequences of chain length l over an alphabet with letters. 6 At the mutation rate p =p , correct and wrong digits are incorporated with the same frequency. The uniform distribution is the dynamical answer to the absence of any preference for nucleotide assignment: All subspecies in the stationary distribution occur with the same frequency. The quasispecies in the intermediate range is determined by the fitness landscapethe distribution of fitness values f i in sequence space, by the move set of allowed mutations as well as the mutation rate p. Typically, sharp transitions occur at some critical mutation rates p = p tr : The distribution changes smoothly from p = 0 to p = p tr , where the distribution turns abruptly into another distribution. At the transition with the largest p-value p = 1∕ > p cr > p tr , the quasispecies becomes an approximate uniform distribution. An explanation is straightforward: Above this critical transition, mutations occur too often to sustain sufficiently accurate reproduction of the template sequence and the result is random replication: In the long run, every sequence is obtained with the same probability. The transition has been characterized as error threshold (Eigen 1971;Eigen and Schuster 1977) since evolutionary dynamics does not sustain a structured longtime population at higher error rates, p > p cr . On typical fitness landscapes, the error threshold sharpens with increasing chain length l and becomes a first-order phase transition in the limit l → ∞ (Tarazona 1992;Park et al. 2010;Huang et al. 2017). 7 Discrete quasispecies In the continuous description, the quasispecies contains all species at finite positive concentrations no matter how small the concentrations might be. Dealing with less than one molecule per reactor volume, however, is unrealistic. Numbers of molecules are non-negative integers, X i ∈ ℕ , subspecies distributions are truncated in reality and we call them stochastic or discrete quasispecies: The discreteness of the stochastic variables leads to a modification of the scenario for the mutation dependence of quasispecies (p) . In the mutation-free system, p = 0 , survival of the fittest is observed and the quasispecies consists of a single sequence: ̃ = X m m = N m with m indicating maximal fitness, f m = max{f i ;i = 1, … , n} , and N being the population size. At sufficiently small values of the mutation rate parameter p, no subspecies except the master exceeds the threshold x i ≥ 1 and the discrete quasispecies ̃ still consists of a single fittest genotype m . The scenario in the small mutation rate regime-populations are almost always homogeneous except short non-stationary periods during which advantageous mutations take over the population-is tantamount to the strong selection-weak mutation scenario discussed in evolutionary biology (Gillespie 1983;Joyce et al. 2008). 8 If the p-value is so large that for one or The value = 2 is used here and applies for binary sequences. For the natural AUGC -nucleotide alphabet we have = 4. 7 We mention that some model landscapes like the additive landscape or the multiplicative landscape sustain smooth rather than sharp transitions (Wiehe 1997). For a recent update and a review of the state of the art in the relation between fitness landscape and selection dynamics see (Schuster 2016d). 8 Biologists (Dean and Thronton 2007;Sniegowski and Gerrish 2010) and computer scientists commonly distinguish strong and weak mutation. The weak mutation scenario assumes that adaptive mutations are sufficiently rare and do not interfere with the selection process but initiate replacements of currently fittest genotypes by still fitter variants. The strong mutation scenario is characterized by sufficiently large values of p that give rise to the quasispecies dynamics described here [for details see (Domingo and Schuster 2016)] and to mutation induced cooperative dynamics ("Cooperation and mutation").
more subspecies x i ≥ 1 is fulfilled, selection of a family of genotypes is observed: ̃ consist of a master genotype m , together with a stationary distribution of sufficiently frequent mutants i ( i ≠ m ). The quasispecies becomes broader with increasing mutation parameter p until a threshold value p cr is reached above which error propagation does not sustain a stationary state and the population Π drifts randomly through sequence space (Huynen et al. 1996;Higgs and Derrida 1991). Interestingly, the critical mutation rate p cr can be derived from the continuous quasispecies theory.
As calculations of the time dependencies of the first two moments of the probability distribution of subspecies in the population, P (Π) (t) , reveal (Jones and Leung 1981) (16') is only marginally stable: The time derivative of the first moment vanishes as expected since the condition of a constant population size is fulfilled by the differential equation (16'). The variance, however, increases with time since the integrand in (18b) is always positive. After sufficiently long time, the fluctuation band becomes so large that the expectation value is irrelevant for the description of the system. The instability in Eigen's quasispecies equation (Eigen 1971) is well known from similar problems: (i) the neutral birth-and-death process with equal birth and death parameters, = , and (ii) the Wiener process. In both cases a constant expectation value is jeopardized by a variance that grows linearly with time. Stability against fluctuation is easily introduced into (16'): one needs only to give up the condition of strictly constant population size and to replace the denominator in (t) by a constant Θ, where Θ is the population size towards which the population converges after sufficiently long time. The approach chosen here, the implementation of a flow reactor as a physical device, yields a stable system as well.
Fluctuations at small particle numbers have different origins ("Autocatalysis in the batch reactor"): (i) conventional thermal fluctuations, (ii) enhanced fluctuations related to autocatalytic self-enhancement, and (iii) anomalous fluctuations in the stochastic variables arising from two or more quasistationary states. The standard deviations (t) fulfil the √ N-law for the resource A but are larger for the autocatalysts 1 , 2 , 3 , and 4 . The stationary states of the stochastic system are extinction S 0 and quasispecies S 1 . Fig. 6 shows a typical example: The anomalous fluctuations in the upper plot are in full analogy to first-order autocatalysis (Fig. 4). Since the initial condition X 1 (0) = 1 was chosen, one outflow step may extinguish the population and the probability of dying out is indeed as large as 20%. The fluctuations bands are extremely broad, and large differences between the deterministic solution and the corresponding expectation values are observed. In the lower plot initial conditions X 1 (0) = 10 were chosen, which are sufficient to reduce the contributions of anomalous fluctuations practically to zero. Then, for a population size of N = 2000 the concentration a(t) coincides with the expectation value E A(t) for all practical purposes and the fluctuations fall into a typical ± √ N -band. The autocatalysts, X 1 , … , X 4 , show broader than usual fluctuation bands because of self-enhancement as we saw in the intermediate range of the first-order autocatalytic reaction (Fig. 2). For long-times the standard deviation (t) stays large in the quasispecies because the flow reactor is an open system and does not approach an equilibrium state (see also Fig. 4).
It is worth recalling what means stochasticity for quasispecies: (i) continuous concentrations are replaced by discrete particle numbers, (ii) fluctuations replace single line trajectories by bands within which trajectories follow a probability distribution, (iii) subspecies can be diluted out of the flow reactor and if this happens for all subspecies the population goes extinct giving rise to anomalous fluctuations, and (iv) error thresholds introduce random reproduction that is closely related to Motoo Kimura's random drift. An increase in the error rate up to the error threshold leads to broadening of the mutant spectrum surrounding the master sequence. Above the thresholds, the populations migrate by random drift through sequence space.

Competition and cooperation
The kinetic equations for replication describing the template induced, uncatalyzed and catalyzed processes are obtained from Eq. () by neglect of mutation. Then the kinetic differential equations for competition and cooperation of competitors result from (9) by setting Q = : Both equations contain terms of different molecularity in and this has the consequence that the dynamic behavior depends strongly on the total concentration c = ∑ n i=1 x i , which in turn is determined by the amount of available resources, a 0 : At sufficiently low concentration, the firstorder terms, k i x i , dominate whereas the second-order terms, l i x i+1 x i become important at high concentrations. No direct analytical solutions are available but exhaustive qualitative analysis is possible and the concentrations of the stationary solutions, a, x i (i = 1, … , n) , can be computed analytically Schuster (2016a, d) ( Table 2).
The basis of the calculations of stationary states is the factorability of the r.h.s. of (19b). The relation x i ⋅ (k i + l i x i+1 )a − r = 0 with i mod n is compatible with stationary states that fulfil either and the conservation condition a + ∑ n i=1 x i = a 0 . In total there are 2 n possible states out of which a single one is asymptotically stable for a given set of parameters. In Table 2 the results for n = 3 are shown. The number of subspecies present at the stationary state, n S , is suitable for characterizing the state: n S = 0 is the state of extinction, n S = 1 is the state of selection, n S = 2 is the state of exclusion, etc., and finally n S = 3 = n occurs at the state of cooperation where all three subspecies are present. The numbers of long-time subspecies n S depend on the resource input a 0 . The relative size ordering of the parameters determines the identity of the selected, excluded, etc., subspecies. For the calculations shown in Table 2 the orderings k 1 < k 2 < k 3 and l 1 > l 2 > l 3 were applied and hence 3 is selected and 1 is excluded. Considering steady states as functions of the available resources, the state of extinction S 0 comes first at small a 0 and is stable for a 0 < r∕k 3 , the state of selection of 3 , S (3) 1 , is stable in the range r∕k 3 < a 0 < r∕k 3 + (k 3 − k 2 )∕l 2 , followed by exclusion of 1 , S (1) 2 , in the range k 3 + (k 3 − k 2 )∕l 2 < a 0 < k 3 + (k 3 − k 2 )∕l 2 + (k 3 − k 1 )∕l 1 . Above this value for a 0 cooperation of all three subspecies is observed provided the flow rate is not too large:

The free concentration of A is obtained as solution of a quadratic equation 9
where the minus sign corresponds to the cooperative state S 3 . The second solution belongs to an unstable state S ′ 3 , which separates the basins of attraction of the states of cooperation and extinction. Above the critical flow rate, r > r cr the states S 3 and S ′ 3 do not exist. For n > 3 , the situation becomes more complex since solutions may oscillate. Many systems with n = 4 have oscillations with very weak damping factors, n ≥ 5 commonly leads to undamped oscillations. The properties of these systems have been discussed extensively in previous publications to which we refer here (Schuster 2016e;Schuster and Sigmund 1985;. Stochasticity has common effects on competition-cooperation systems like thermal fluctuations and fluctuation through autocatalytic enhancement. In addition, there are many more quasistationary states than the asymptotically stable states of the deterministic system. For example, states in which less efficient subspecies are selected show up as quasistationary states as well. In case of the smallest possible system with n = 2 and k 2 > k 1 , we have four states: (i) the absorbing boundary as state of extinction S 0 , (ii) the state of natural selection S (2) 1 where the fittest variant 2 is  (2016) 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. The abbreviations 3 = r∕k 3 , 31 = (k 3 − k 1 )∕l 1 ) , and 32 = (k 3 − k 2 )∕l 2 ) are used for the combination of parameters. For the cooperative state S 3 the stationary concentration of A is obtained as one root of a quadratic equation (20) with two combinations of the rate parameters, = ∑ 3 i=1 k i ∕l i and = ∑ 3 i=1 1∕l i , and i = (r − k i )∕(l i ) for i = 1, 2, 3 and = a 1 from (20). The existence of the non-trivial stationary state requires a sufficiently small flow rate: r ≤ (a 0 + ) 2 ∕4

Symbol
Stationary values Stability range It is worth mentioning that Eq. (20) and the solutions for a 1,2 are valid for arbitrary n provided the summations in and are adjusted accordingly.

3
selected, (iii) the state of selection of the less fit subspecies 1 , and eventually (iv) the state of cooperation S 2 . The relative stabilities of the individual states are reflected by the probabilities to reach these states by randomly chosen trajectories (Table 3). Parameters for the calculations shown in the table were chosen such that the corresponding deterministic system is situated in the cooperative domain in parameter space, and indeed the approach towards the cooperative state S 2 has always the largest probability. Again, initial particle numbers around X 1 (0) + X 2 (0) = 4 are sufficient for strong dominance of the state corresponding to the deterministic solution. In case of n = 3 we present expectation values of the four stochastic variables, A(t) and X i (t) ( i = 1, 2, 3 ) at some predefined time of the simulation end, t end for different initial conditions and mutation rate parameters p. The values for p = 0.0 refer to the pure competition-cooperation case discussed here (Table 4). As expected extinction plays a major role at small initial values, X i (0) = 1, 2, 3 ( i = 1, 2, 3 ), but X i (0) = 4 is already sufficient for coming close to the expectation values obtained for large initial values, X i (0) = 10 is enough for reaching the deterministic values for practical purposes. Table 3 Probabilities to reach quasistationary states in the cooperative regime with n = 2 and different initial conditions The table provides probabilities of occurrence for all four possible long-term states: extinction S 0 , selection of 1 in S (1) 1 , selection of 2 in S (2) 1 , or cooperation S 2 . The counted numbers of events 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: Probabilities are obtained by multiplication by 10 −4

Initial values
Counted numbers of states in final outcomes 28.0 ± 5.0 9959.9 ± 6.4 5 5 0 2.5 ± 1.1 6.3 ± 2.6 9991.2 ± 3.0 Table 4 Long-time behavior in the cooperation-mutation system with n = 3 and different initial conditions The table provides expectation values at the time of the end of the simulation ( t end = 30 ) for different mutation rate parameters p and different initial conditions. Choice of parameters: l 1 = 0.011 , l 2 = 0.010 , In systems with n ≥ 4 deterministic and stochastic solution curves oscillate. The solutions of the ODE's are different for n = 4 where weakly damped oscillations occur and for n ≥ 5 showing undamped relaxation oscillations (Phillipson et al. 1984). In the stochastic approach, the systems die out after population numbers of individual subspecies went beyond X(t) = 1 , and for sufficiently large population sizes, four-membered systems may survive for very long time whereas systems with five or more members go extinct. Cases with n = 5 are well suited for studying the transition from selection to cooperative dynamics through increase of the parameter ration ratio h∕f = l∕k (Fig. 9). At dominant competition or small h-values the system approaches selection of the fittest as long-time solution. The upcoming role of cooperation in a series of systems with increasing parameters l can be nicely illustrated by a series of plots of trajectories from selection to somewhat chaotic looking intermediate scenarios and further to oscillatory hypercycle dynamics (see Fig. 9 where a nonzero mutation rate parameter p was chosen and where accordingly quasispecies formation instead of selection is observed).

Cooperation and mutation
The combination of cooperation and mutation reveals a less common role of mutation in addition to the creation of diversity through variation. In principle, mutation can reintroduce extinguished subspecies into the population. Here, we shall focus on this aspect and, in particular, study the influence of the mutation rate parameter p on extinction times. To study the role of mutation in low-dimensional cooperative systems ( n = 2, 3 ) expectation values of the stochastic variables A(t), X 1 (t) and X 2 (t) , or X 2 (t) and X 3 (t) were calculated at predefined times ( t end ) and compared with the probabilities of trajectories to end up in the cooperative state, S 2 or S 3 , respectively. For the initial conditions X 1 (0) = X 2 (0) = X 3 (0) > 4 the results are practically indistinguishable from the deterministic values. An increase in the mutation rate parameter p shows the expected influence: Extinguished subspecies can be reintroduced and this increases the probability of reaching the cooperative state, S 2 or S 3 . The case n = 3 is shown as an example in Table 4.
The oscillating systems are more difficult to investigate. Here, we consider the time of extinction of the entire population as a function of the mutation rate and the available resources, T 0 (p, A 0 ) . The results are shown in Fig 8: The extinction times T 0 show very strong scatter and their appearance is dependent on the resolution of the calculations. By resolution, we mean here the number of molecules A in A 0 between two neighboring points. The highest resolution is achieved when the calculations are performed with every (integer) number of A 0 molecules, e.g., ΔA 0 = 1 yields 100, 101, 102, … . Computations at somewhat lower resolutions are less time consuming and provide in essence the same results. The plots shown in Fig. 8 show enormous scatter but, nevertheless, allow for drawing two conclusions: (i) In the mutation-free case the extinction time T 0 is independent of the amount of available resources up to a value A 0 ≈ 1300 , and (ii) for non-zero mutation rates a kind of noisy or stochastic threshold phenomenon is observed. Considering the noisy function T 0 (p, A 0 ) and taking A 0 at the first value of T 0 (p, A 0 ) ≥ 1000 we find for the parameter values applied in Fig. 8: A (T 0 ≥1000) 0 = A (cr) 0 = 1130, 690, 360 for the mutation rates p = 0.0005, 0.0010, 0.0020 , respectively. As expected the threshold moves to lower A (cr) 0 -values with increasing mutation rate. The behavior of the extinction times T 0 (A 0 ) is Fig. 8 Times to extinction as a function of available resources in the five membered cooperative system ( n = 5 ). Extinction times T 0 of the populations Π are shown for different amounts of available resources measured as inflow concentrations a 0 or A 0 when expressed in numbers of molecules per unit volume. The upper diagram presents the data at a resolution of ten molecules ( ΔA 0 = 10 ; 100, 110, 120, … ) for four different values of the mutation rate parameter: p = 0.0000 (red), p = 0.0005 (yellow), p = 0.0010 (green), and p = 0.0020 (blue). T 0 -values above 1000 are truncated at this value. The lower diagram shows the two plots p = 0.0000 (red) and p = 0.0020 (blue) at the highest possible resolution ( ΔA 0 = 1 ). Choice of other parameters: r = 0.5[V −1 t −1 ], l 1 = l 2 = l 3 = l 4 = l 5 = 0.01[M −2 t −1 ]. Pseudorandom number generator: Extended CA, Mathematica, seed: s = 491 . Initial conditions: A(0) = 0 , X 1 (0) = X 2 (0) = X 3 (0) = X 4 (0) = X 5 (0) = 5 1 3 similar for n = 4 but the critical concentrations A (cr) 0 for the different p-values lie much closer together and the analysis is more difficult. Considering survival at constant resources A 0 reveals a mutation threshold above which the population survives to long times.
Hypercycle extinction is an example that reflects well the expected increase in lifetime with increasing mutation rate. One general remark nevertheless is important: This mechanism of reintroduction of extinguished subspecies requires that template and mutant are close relatives and that the Hamming distance between them is not too large. What we need in reality, however, is not a perfect revertant being genetically identical to the lost original, we need only a subspecies that can replace the original with respect to its phenotypic function. Suppression of deleterious mutations (Gorini and Beckwith 1966;Hartman and Roth 1973;Prelich 1999) as well as the relation between protein sequence, structure and functional efficiency Knowles 1976, 1977) have been extensively studied in the last decades of the twentieth century.

The complete model
Completion of the model brings together the three faces of the coordinate system in Fig. 5 and is concerned with an analysis of the dynamics in the interior. An appropriate strategy for analyzing the interior consists in choosing certain type of behavior on one of the three faces and increasing the third parameter from zero to the value of interest. Raising the third parameter will change the dynamic behavior either gradually or in threshold-like manner or stepwise through a cascade of bifurcations. Illustrative prototype examples are seen through rising the mutation rate in competitive or cooperative reproduction, or with the introduction of cooperation into Darwinian systems.
The characteristic dependence of the population dynamics on n, the number of subspecies, prevails also in the full model. For small numbers ( n = 2, 3 ) and p = 0 the transition from the competitive to the cooperative system has been discussed in "Competition and cooperation". An increase of the cooperation parameter h = l A(t) leads in steps from selection of the fittest to a cooperative state with all subspecies present. Oscillating systems ( n ≥ 4 ) are more spectacular since the hypercycles are unstable at p = 0 in the stochastic approach and raising the cooperation parameter h leads from selection to extinction. In the intermediate parameter range where the deterministic system shows stepwise increase in the number of coexisting subspecies ( 1, 2, … , n or, expressed phenomenologically, selection, exclusion, … , cooperation) the stochastic approach yields highly irregular dynamics with different numbers of non-extinguished subspecies whereby the number of species present increases with increasing values of the cooperation parameter h. The scenario in parameter space is completed by considering the series of states at different mutation rate parameters p. The case n = 5 , which for large h-values leads to undamped oscillations with a stochastic contribution to the amplitude, was chosen to facilitate the illustration ( Fig. 9): At zero mutation rate p = 0 the series with increasing h is selection → irregular dynamics with two species → irregular dynamics with three species → … → extinction. At intermediate p-values, we find quasispecies → irregular dynamics → oscillations with highly irregular spacings. At high mutation rates, we are dealing with quasispecies → irregular dynamics with increasing numbers of dominant subspecies → stochastic hypercycle dynamics (Fig. 9).

Concluding remarks
The model presented here has been conceived with modular structure in the sense that different mechanisms can be applied for each of the three basic components. Here, it has been presented in its simplest form. Each of the three modules, competition, cooperation and variation, can be made arbitrarily complex. Variation, for example, can be extended to include more elaborate mutation mechanisms and recombination as well as environmental influences. Even in case of viruses the reproduction mechanism is commonly much more elaborate and comprises a whole molecular machinery instead of a single enzyme. Virus reproduction may include also the defense system of the host, epigenetic phenomena could be taken into account through simultaneous consideration of several generations, and for higher organisms the real challenge in reproduction is to deal with the enormous complexity of development in a form that is simple enough for modeling. Cooperation at the molecular level could also involve reproductive autocatalytic networks whereas social phenomena in reproductive groups or societies represent the currently highest step in the open ended complexity increase of biological evolution. Cooperation has been frequently modeled by game theory Maynard Smith (1982); Hofbauer and Sigmund (1998). There is no limitation to make the model more complex, the problem evidently is to include the desired phenomena but to keep the model sufficiently simple for mathematical analysis or simulation.
In the simple form in which the model was introduced here, it has been tested experimentally by in vitro evolution experiments [For an overview of early works on this subject see (Spiegelman 1971;Biebricher 1983); as a recent review we mention (Joyce 2007)]. The kinetic equations describing replication and mutation were introduced 1971 by Manfred Eigen in his scholarly written paper on self-organization of biological macromolecules (Eigen 1971). Eigen's mutationselection equation describes the evolution of the distribution of asexually reproducing genotypes in a population of constant size N. Correct replication and mutation are seen as parallel chemical reactions leading to a uniquely defined stationary population called quasispecies (Eigen and Schuster 1977). RNA replication catalyzed by single virus specific enzymes from RNA bacteriophages provides a bridge from chemistry to biology: The mechanism of the replication process is well understood in all molecular details (Biebricher 1983;Biebricher et al. 1984Biebricher et al. , 1985 and an appropriate replication assay serves for in vitro evolution studies (Mills et al. 1967;Biebricher 1983). The mutation-selection scenario was found to provide an appropriate molecular basis for understanding also virus evolution [For a recent survey see the contributed volume (Domingo and Schuster 2016)]. More complex systems, for example, bacteria and populations of cancer cells, were found to be describable by quasispecies theory as well (Bertels et al. 2017;Covacci and Rappuoli 1998;Napoletani et al. 2013;Brumer et al. 2006).
In the strong mutation scenario 8 Darwin's view of evolution has to be modified. Not a single fittest genotype is selected but a uniquely defined distribution of genotypes, which is represented by the largest eigenvector of a value matrix that represents the long time or stationary solution of Eigen's mutation-selection equation. The mean fitness of populations is not always optimized since situations can be constructed in which the fitness is decreasing in the approach towards the stationary state. A trivial but illustrative example of decreasing fitness during evolution considers a homogeneous population consisting exclusively of fittest genotypes as initial condition: Mutations introduce mutants into the population and since they have lower fitness by definition the mean fitness is doomed to decrease. Such situations, however, are rather rare and Darwinian optimization still remains a very powerful heuristic that applies to almost all scenarios. For error rate parameters exceeding a critical value p cr , the largest eigenvector approaches the uniform distribution over the entire sequence space, which is the exact solution for the value p =p leading to incorporations of correct and incorrect nucleotides with equal probabilities-for binary sequences this happens at p = 1 −p = 1 2 . In realistic populations, the uniform distribution is incompatible with a discrete quasispecies (17). Instead populations are observed that migrate randomly through sequence space Higgs and Derrida (1991); Huynen et al. (1996).
In the second half of the twentieth century, most of the molecular insights into reproduction and inheritance came from viruses and bacteria and a high percentage of molecular biologists thought that the basic regulation mechanisms of gene activities are understood. Eukaryotic cells, however, are not "giant bacteria". Although the genetic code is the same, the gene expression and inheritance system of higher organisms are different from the prokaryotic one and much more complex. A true wealth of information on eukaryotic cells has been discovered in the past fifty years, but gene expression in animals, plants, and fungi is still a subject of current cutting-edge research. Most of the recently revealed gene expression regulation mechanisms are subsumed under the notion of epigenetics for which an operational definition has been proposed at the Cold Spring Harbor Meeting in the year 2008 (Berger et al. 2009): An epigenetic trait is a stably heritable phenotype resulting from changes in a chromosome without alterations in the DNA sequence.
The diversity of epigenetic effects on gene regulation is enormous. It ranges from specific methylation of DNA, in particular cytosine methylation in position 5 of CpG elements (Zemach et al. 2010), histone methylation and acetylation (Lawrence et al. 2016) to post-transcriptional RNAmethylation of adenine in position 6 (Barbosa Dogini et al. 2014;Yue et al. 2015) and small interfering RNAs (He and Hanon 2004). Epigenetics provides an extremely diverse, complex and flexible richness of regulatory actions on genes, which so far was not yet cast into a comprehensive theory and precisely this is one of the greatest challenges for the future of evolutionary biology.
There is neither a convincing theoretical model nor experimental evidence that Darwinian evolution leads to an obligatory increase in complexity. The combination of competitive selection and cooperation, however, may lead from one level of complexity to the next higher one by integration of competitors through cooperation (Szathmáry and Schuster 1996). The evolution model presented here proposes a mechanism for this integration of competitors and identifies the abundance of resources as one driving force towards higher complexity. This simple model distinguishes four steps (Schuster 1996): (i) Initially the systems consists of independent replicators competing for a single resource, (ii) the capability of cooperative interaction allows to form an autocatalytic network, which couples the replicators and suppresses competition but makes the network vulnerable to exploitation by parasites, which consume resources without contributing a share to the common properties, (iii) the members of the autocatalytic network are separated from the environment by means if a suitable boundary that prevents the system from exploitation and allows for the formation of a new unit at a hierarchically higher level, and (iv) the individual units at the higher level diversify by variation, compete for common resources, Darwinian selection sets in and takes place now a the higher level. The previously autonomous units at the lower level lost their autonomy at least in part when they were integrated into the higher unit of selection. Although modeling major transitions as shown in "Competition and cooperation" is not difficult, suitable experimental molecular models are very hard to conceive, because second-and higher-order autocatalytic systems consist almost always of complex reaction networks rather than single-step reactions [as examples for attempts to construct simple systems of this kind see (McCaskill 1997;Wlotzka and McCaskill 1997)]. John Maynard Smith and Eörs Szathmáry collected a true wealth of evidence for the historic occurrence of such major evolutionary transitions  in the evolution of life.
It is illustrative to think about transitions in terms of thresholds: Up to a certain value of the transition determining parameter the typical feature is hardly detectable and does not become evident before the parameter exceeds the transitions value. Accordingly, such thresholds correspond to sharp transitions. Nevertheless, it appears useful to be less fussy and to accept the notion of threshold also for smooth transitions. On the three faces of the coordinate system (Fig. 5) we observe an error threshold between selection and random replication, a cooperation threshold between selection and symbiosis, and a mutation threshold that separates the regime of independent replication of all subspecies from mutual support through frequent mutation.
Understanding evolution implies knowledge on the relation between genotypes being DNA or RNA sequences and phenotypes giving rise to all fitness relevant parameters. The metaphor of an abstract fitness landscape has been originally introduced by Sewall Wright for the purpose of illustration (Wright 1922). In formal mathematical terms such a relation can be modeled as a mapping from a genotype or sequence space into fitness values. In molecular structural biology such a mapping is split into two parts: (i) a mapping from sequences into structures or genotypes into phenotypes, and (ii) a second mapping assigning fitness values to structures or phenotypes . Computer models of RNA evolution with sequence-secondary structure-fitness mappings have been studied in the past (Fontana and Schuster 1987;Fontana et al. 1989;Fontana and Schuster 1998a, b and these studies provided the basis for a definition of evolutionary nearness of phenotypes in the presence of neutrality (Stadler et al. 2001). With more and more data on sequences and fitness values of mutants becoming currently available fitness landscapes may also be determined directly by experiment (Kouyos et al. 2012) and it is not risky to predict that genotype-phenotype relations will become an integral part of evolution models in the near future. Then evolution can be described in a self-contained manner where the genotype-phenotype mapping is a genuine part of the evolving system. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.