Molecular evolution between chemistry and biology

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, ( ) competition and mutation, (  ) cooperation and competition, and (  ) cooperation and mutation, are analyzed. Changes in population dynam- ics that are induced by bifurcations and threshold phenomena are discussed. Keywords Competition · Cooperation · Darwinian optimization · Error threshold · Mutation · Neutral evolution · Selection Introduction The model focusses on three basic processes: (i) fitness- driven competition through differential reproductive success, Biology is centered around evolutionary thinking or under- (ii) reproduction-relevant cooperation between competitors, standing biology implies understanding evolution as Theo- and (iii) reproduction-induced variation. The model is con- dosius Dobzhansky pointed out clearly in his famous book ceived with a modular structure and allows for the imple- on evolution: “Nothing in biology makes sense except in the mentation of different, more or less complicated mechanisms light of evolution” (Dobzhansky et al. 1977). Pure Darwin- for the processes (i), (ii), and (iii). For example, variation ian evolution is a simple process but its embedding in nature may be implemented by mutation, by recombination or by renders it complex: Natural selection would follow uncom- both. Here, we shall apply the simplest conceivable mecha- plicated laws in a simple environment. In the light of current nisms: reproduction as single enzyme mediated replication molecular biology, there is need for a simple but comprehen- (Biebricher 1983), cooperation of competitors as hypercy- sive mathematical model of evolution to be able to account cle dynamics (Eigen and Schuster 1978a, b), and variation for modern genetics. The various epigenetic mechanisms as single point mutations based on the uniform error rate have to be part of any comprehensive model of evolution assumption (Swetina 1982). and to keep such a model amenable to analysis and handling, Usage of the notion evolution is often ambiguous and a molecular details must be reduced to a coarse-grained level. precise definition is desirable. Here evolution is understood This article deals with a flexible model of evolution under as a process based on reproduction of a genotype being a defined environmental conditions. We present a concise and DNA or an RNA sequence that carries encoded information comprehensive review of work that was published elsewhere on the formation of a phenotype, which is evaluated with (Schuster 2016a, d, 2017a) together with a few new results. respect to success in reproduction. The evolutionary pro- cess is built upon two foundations: (i) the dynamics at the population level and (ii) the environment dependent encod- Special Issue: Chemical Kinetics, Biological Mechanisms and ing of phenotypes in genotype. At the molecular level, the Molecular Evolution. latter boils down to sequence–structure–function relations * Peter Schuster (Schuster 2016). The simplest systems that are capable of pks@tbi.univie.ac.at evolution in the sense of the given definition are nothing but special autocatalytic reactions involving polynucleotides Institut für Theoretische Chemie, Universität Wien, under suitable conditions. Währingerstraße 17, 1090 Wien, Austria Vol.:(0123456789) 1 3 404 European Biophysics Journal (2018) 47:403–425 In the next "Preliminaries", some prerequisites are pre- be implemented as long as they can be cast in kinetic equa- sented for the model, which is introduced in "A minimal tions. In a separate "Basic processes", the processes that will mathematical model for the evolution of molecules". The be used to model evolution are introduced. The following model comes in a deterministic version based on kinetic dif- section describes the deterministic and stochastic methods, ferential equations or it is formulated as a stochastic process which are applied to find solutions of the kinetic equations. modeled by means of chemical master equations (Schuster Finally, we review some fundamental features of autocataly- 2016). Although almost no analytic solutions are available sis since this is the chemical counterpart of reproduction. for master equations derived from nonlinear reaction kinet- ics in two or more variables, the stochastic version of the Idealized environments model can be studied by efficient simulation methods for the calculations of trajectories (Gillespie 1977, 2007). In Environments that allow for investigations of observa- three sections, solution curves for the three two-dimensional tions as functions of one or few parameters with everything subspaces, (  ) reproduction and mutation, (  ) reproduc- else being constant require elaborate design in the form of tion & cooperation, and ( ) cooperation & mutation, are pre- sophisticated experiments since devices controlling envi- sented, analyzed and interpreted. Then follows a brief dis- ronmental conditions may be quite involved. In theoretical cussion of results obtained with the full three-dimensional approaches, often the silent assumption is made that there model, reproduction & cooperation & mutation and in the exists a hypothetical machinery, which takes care of fixing final section we discuss the simple model in the context of parameter values as needed for the mathematical analysis. the complex processes in nature. 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 evolu- Preliminaries tion, which result from reproduction, symbiotic cooperation and variation, and therefore impacts on evolution caused by Darwinian evolution is often—and incorrectly—seen as a changes in the environment are intended to be kept to a mini- synonym for optimization mainly because Ronald Fisher’s mum. To reduce the influence of the environment as much fundamental theorem of natural selection (Fisher 1930, as possible we shall assume a control device in the form of 1958; Price 1972; Ewens and Lessard 2015) focusses on a simple flow reactor (Fig.  1; see also Schmidt 2005). More non-decreasing mean fitness (t) in evolving populations. In elaborate reactors, which keep, for example, the numbers of the simplest form derived from idealized models, the time bacterial cells constant have been designed and implemented derivative of the mean fitness is the variance of the fitness (Novick and Szillard 1950a, b; Bryson 1952). The flow reac- values and therefore a non-negative quantity. In reality, mean tors called chemostat, cellstat or turbidostat and other exper- fitness is the result of many factors and as Fisher was been imental devices for monitoring and controlling evolution in certainly aware, only a few cases—like single locus genet- the laboratory may serve as examples (Husimi et al. 1982; ics—fulfil the theorem in pure form (for elaborate discus- Dykhuizen and Hartl 1983; Husimi 1989; Koltermann and sions see the more recent literature Plutynski 2006; Okasha Kettling 1997; Strunk and Ederhof 1997). 2008). Implementation of a physical device rather than appli- Evolution is intimately related to the environment in cation of idealized assumptions like constant population which it takes place and accordingly environment and envi- size is required for the stochastic description of evolution. ronmental changes are major factors shaping evolutionary An illustrative example where the deterministic approach processes. Here, we are primarily interested in the internal yields a stable solution whereas the corresponding stochastic dynamics and hence a well-defined and controllable environ- system is unstable is provided by the linear birth-and-death ment is required. For this goal, we introduce a flow reactor in process (Goel and Richter-Dyn 1974), which is described by f d "Idealized environments", which represents a simple device the kinetic equations  + ⟶ 2 and ⟶∅ . For equal that is not only useful for performing evolution experi- birth and death rate parameters, f = d , the population num- ments but also provides at the same time a suitable setup ber of the deterministic system stays constant whereas in for theoretical modeling. More complex environments can 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 The term linear reaction kinetics is used in this contribution for words, the system is unstable despite a (marginally) stable reactions and flow terms of zero and first order, which give rise to deterministic solution. The direct incorporation of constant constants or linear functions in the kinetic ODEs. First-order kinetics population size into Fisher’s selection (1930) or Eigen’s gives rise to exponential functions, which become linear in ln[]∕t -plots. equation (1971) leads to an instability of the same kind since 1 3 European Biophysics Journal (2018) 47:403–425 405 in the replication of every molecular species  —the number of catalytic terms is n 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, 1978, 1978): The catalyzed reactions 2 +  ⟵ +  +  (i = 1,… , n; i mod n) i i+1 i i+1 ⋯ ←  ←  ←  ← ⋯ ←  ←  ← … n 1 2 n−1 n 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 sim- plest 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 replica- tion kinetics of the multistep process, Q ⋅k [] ji i (1) +  + →→→ +  +  , i j i Fig. 1 The continuous-flow stirred-tank reactor (CSTR). The fig- ure sketches a device for controlling the environmental condition with E being a replicase enzyme, the reaction rate is obtained of evolution experiments. The material needed for reproduction is as the product of two parameters: Q ⋅ f where f = k [] is ji i i i subsumed by A, it flows into the reactor with a (volumetric) flow the fitness of the template that depends on the availability rate r   [V  /  t]≡[volume  /  time] in form of a solution with concentra- tion []= a or number density []= A . In the reactor molecules of resources here the concentration [ A]. The dimensionless 0 0 (i = 1, … , n) are reproduced and A is consumed. The volume V factors Q with i, j = 1,… , n are understood as the elements ji of the reactor is constant and hence reaction mixture compensating of a mutation matrix Q ={Q } and provide the probability ij the volume increase through influx of stock solution has to flow out to obtain  as a copy of the template  that can be either of the reactor. The mean residence time of a volume element in the j i CSTR is  = V∕r [t]. Inflow and outflow of materials are handled as correct (  ≡  ) or error-prone (  ,  j ≠ i ). Conservation of a ⋅r j i j virtual chemical reactions or pseudoreactions: ∗⟶ for inflow and, r r probabilities requires: Q = 1 since each copy has to be ji j=1 ⟶∅ and  ⟶∅ (i = 1, … , n) for outflow either correct, Q , or incorrect, Q = 1 − Q . For ii ji ii j=1,j≠i the sake of simplicity, binary sequences will be considered here and we distinguish only between correct and incorrect fluctuations are not self-regulating (Jones and Leung 1981). base pairs. Implementation of the system in the flow reactor provides A useful simplifying approximation is made by the uniform stability due to balance control of inflow and outflow mod- a ⋅r 0 error rate model (Swetina 1982): The error per nucleotide and eled by the pseudoreactions and ⟶∅ as well as ∗⟶ replication event, p, is assumed to be independent of the posi- ⟶∅ (i = 1, … , n) (Fig. 1). tion of the nucleotide along the sequence and the nature of the nucleotide to be complemented. Then all elements of the Basic processes mutation matrix can be expressed by a simple formula, The minimal system for modeling evolution of molecules d l−d l d ij ij ij Q = p (1 − p) = p  with  = , ji (2) is based on the three classes of processes: (i) competitive 1 − p reproduction, (ii) symbiotic cooperation, and (iii) repro- with only three parameters: (i) the sequence length of duction based variation. For the minimal model we shall the RNA molecules l, (ii) the mutation rate p, and (iii) choose the simplest possible chemical reactions. Reproduc- the Hamming distance (Hamming 1950, 1986) between tion will be modeled by simple replication,  +  ⟶ 2 i i the two sequences interrelated by the mutation process, with f being the fitness of species  and competition being i i d = d (X , X ) . Without changing important results for the ij H i j introduced through living on the same resource A. Symbi- purposes pursued here, the analysis of the model is sub- ontic cooperation is introduced as catalyzed replication, ij stantially simplified by the assumption of constant chain +  +  ⟶ 2 +   (i, j = 1, … , n) , where  repre- i j i j i lengths l, which is also consistent with the restriction to sents the template and  is the catalyst. In its most general form—every molecule  has the potential to act as catalyst 1 3 406 European Biophysics Journal (2018) 47:403–425 point mutations since point mutations do not change chain Equilibrium fluctuations in conventional chemical reac- lengths by definition. tion kinetics follow an approximate N -law and hence play As an alternative to Eq. (1), mutation can be seen as almost no role in systems with particle numbers that are a consequence of DNA change—damage and imperfect typical for chemical systems. In biology a different scenario repair—during the whole life time of an organism, which is encountered. For example, every new variant originat- is the idea in the Crow–Kimura mutation model (2009), ing from mutation has to start out from a single copy. The pp 264–266]: deterministic approach commonly uses continuous vari- ables, which can only be an approximation to reality, since k [] ji numbers of molecules or biological entities are integers by (3) +  + →→→ 2 +  and  →→→ . i i i j definition. Continuous concentrations can adopt arbitrarily Interestingly, both models (1) and (3) although different with small values whereas stochastic variables cannot pass low respect to the underlying physics give rise to identical math- values beyond unity, because then molecular species go ematical problems (see e.g. Baake and Gabriel 1999 and extinct, and deterministic solutions become unrealistic and Schuster 2016, pp 76–78). differ strongly from the stochastic results. Large population size alone is not sufficient, each variable has to be suffi- Deterministic and stochastic approaches ciently large at every instant to guarantee similarity between stochastic and deterministic solutions. Two more phenom- Reaction mechanisms are commonly analyzed by determin- ena are relevant in the stochastic approach at low particle istic and stochastic approaches. The former translate the numbers and may lead to large standard deviations in the chemical reaction equations into kinetic differential equa- variables: (i) nucleation steps in reactions involving two or tions, which can be directly solved by mathematics, studied more molecules and (ii) multiple (quasi)stationary states of by means of qualitative analysis or investigated by numeri- the stochastic system. The stochastic system may approach cal integration. The theorems of existence and uniqueness any stationary state and the distribution of the possible out- of solutions of differential equations are applicable and a comes is determined by probabilities. Then the calculations single integration provides the complete information for a of expectation values and higher distribution moments build given input set. Competitive selection with nonzero mutation averages over different final states and may give rise to enor - leads to one unique asymptotically stable stationary state mous standard deviations. (Thompson 1974; Jones et al. 1976), whereas the long-time For the purpose of illustration, we consider the equilibra- dynamics of cooperative systems is much richer and multiple tion of the flow reactor as an example of a stochastic system stationary states, oscillations or deterministic chaos may be that can be fully analyzed by analytic calculations (Schuster observed (Schuster 1978; Schnabl et al. 1991). 2016, pp 436–441). The two pseudoreactions, Stochastic analysis in general in based on searching for a a r (4a) ∗ ⟶  and stochastic process that fits the model to be studied as closely as possible. Chemical reaction kinetics prefers master equa- r (4b) ⟶ ∅ , tions although the repertoire of analytical solutions is very are converted into the easy to solve kinetic differential limited. It is not difficult to write down a multivariate master equation equation but the derivation of analytical solutions is suc- cessful only in exceptional cases, for example for networks da −r t = a ⋅ r − a ⋅ r =(a − a) r with a(t)= n + n − a e 0 0 0 0 0 of monomolecular reactions (Jahnke et al. 2007; Deuflhard dt (4c) et al. 2008. If no analytical solutions are available informa- where n = a(0) . The master equation for the probability tion on the stochastic system can be obtained by trajectory sampling. The theoretical background for trajectory har- P (t)= P A(t)= n with number density A(t) being the dis- crete pendant of concentration a(t), vesting has been laid down by Andrey Kolmogorov (1931), Willy Feller (1940), and Joe Doob (1942, 1945). With elec- tronic computers now being generally available elaborate = r a P +(n + 1)P −(a + n)P , n ∈ ℕ , 0 n−1 n+1 0 n simulations of stochastic processes became possible. The (4d) more recent conception, analysis, and implementation of a simple but highly efficient algorithm by Daniel Gillespie (1976, 1977, 2007) provides a very useful tool for investiga- A stochastic quasistationary state is a state towards which the sys- tions of stochastic effects in chemical kinetics. Distributions tem converges stochastically in the long-time limit and around which it fluctuates. It is not an absorbing boundary, and if true asymptoti- of trajectories are characterized by expectation values and cally stable stationary states exist the system converges to one of higher moments, commonly only by variances or standard them in the limit t → ∞ although the mean time of convergence may deviations. be extremely long. 1 3 European Biophysics Journal (2018) 47:403–425 407 has the analytical solution The forward reaction of (5) leads to molecular self-enhance- ment: The resource A is consumed in the synthesis of the min{n ,n} −krt −rt n +n−2k replicator X and the process is catalyzed by the presence n e (1 − e ) −rt n−k a (1−e ) P (t)= a e of further molecules X. Depending on the molecularity k (n − k)! k=0 of the autocatalytic process as expressed by the value of (4e) m autocatalysis comes in different forms that, in essence, with n, n , a ∈ ℕ for the sharp initial condition P (0)=  . 0 0 n n,n fall into two classes: (i) simple or first-order autocatalysis In the limit t → ∞ the distribution (4e) converges to a Pois- with m = 1 and (ii) second- and higher-order autocatalysis son distribution with m ≥ 2 [(Phillipson and Schuster 2009, pp 18–27)]. For consistency, we add here also the uncatalyzed reaction and lim P (t)= P = exp(−a ) . (4f) n n 0 t→∞ n! call it zeroth order autocatalytic ( m = 0). All processes in closed systems converge to an equilibrium state and show First and second moment yield expectation values and stand- mass conservation that implies a conservation relation ard deviations A(0)+ X(0)= A(t)+ X(t)= N . The reversible autocatalytic � � −rt reactions converge to the equilibrium state: E A(t) = a +(n − a ) e and 0 0 0 � � (4g) −rt −rt A(t) = (a + n e )(1 − e ) . S ∶ a = N∕(1 + K), x = N ⋅ K∕(1 + K). 0 0 Although the expressions for the equilibria are the same The standard deviation starts out from (t = 0)= 0 and for the all reactions independently of the stoichiomet- approaches the equilibrium value  = a either from ric coefficient m, K = g ∕g , there are subtle differ- → ← below for a > n or from above for a < n after having 0 0 0 0 ences in the probability distributions, which become passed a maximum at important at small concentrations: The particle numbers are discrete quantities, change by integers only, and the 2n t = ln . (4h) max stoichiometric factors are X(t) X(t)− 1 ≡ X(t) or r n − a 2 0 0 X(t) X(t)− 1 X(t)− 2 ≡ X(t) rather than X(t) and In general the maximum if it exists is very flat (see Fig.  3) X(t) , respectively . The states with X(t)= 1 for m = 1 , or the as can be easily checked from the analytical expressions states X(t)= 1 and X(t)= 2 for m = 2 cannot react because through inspection: 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 � � (n − a ) � � n − a 0 0 0 0 states with X(t)= 0 and X(t)= 1 in the first- or second-order E A(t ) = a + and  A(t ) = . max 0 max 2n 2 n autocatalytic reactions require different normalizations for the stationary probabilities of the three systems: The flow reactor is one of the few rather exceptional cases (0) where the stochastic approach can be done completely by P = ; n ∈[0, N] , (6a) (1 + K) mathematical analysis. (1) Autocatalysis in the batch reactor P = ; n ∈[0, N − 1] , (6b) N N (1 + K) − K A batch reactor is an elaborate device that allows for per- (2) forming chemical reactions under controlled conditions N K P = ; n ∈[1, N − 2] , N N N−1 (1 + K) − K − N K − 1 without inflow and outflow (Schmidt 2005). Here, the term batch reactor is used to indicate that reactions are carried out (6c) in a well-mixed closed system under constant temperature where as before the stochastic variable is n = A . Mass con- and pressure, and in the long run approach a thermodynamic servation provides X = N − n . As expected the truncation (0) equilibrium state. of P is important for small values of N only. For first- In conventional chemistry autocatalysis is a rather rare order autocatalysis with N = 10 , k = k = 1.0 we obtain → ← phenomenon but in biology it represents the most impor- tant process since multiplication is just a special form of autocatalysis that in simplest form can be expressed by the We shall use different letters for the rate parameters: For m = 0 we reversible reaction use g ≡ h , for m = 1 g ≡ k , and for m = 2 g ≡ l . It is worth recalling −1 −1 −1 −2 −1 the different dimensions: h [t ], k [M t ], and l [M t ] that are the same in both directions. −−→ + m ←−− (m + 1) . (5) g Here the expression (x) = x!∕(x − n)! denotes the falling Pochham- → n mer symbol. 1 3 408 European Biophysics Journal (2018) 47:403–425 small t the factor A(t) is large and the factor(s) contain- ing 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. Self- enhancement 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 uncata- lyzed 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{(t)} is a measure of the random scatter in the max delay in the onset of the autocatalytic reaction (Table 1). (ad iii) Multiple final states give rise to an additional sto- chastic 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 ini- tial 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, respec- tively. A total population size of N = 1000 was chosen so Fig. 2 Comparison of fluctuations in the reversible zero- and first- that the one-standard-deviation fluctuation band of order order autocatalytic reaction  + m ⇌ (m + 1) with m = 0, 1 bf in N appears small and the deterministic solutions coincide the batch reactor. The two plots show the expectation values E A(t) (black) and E X(t) (red) together with the one standard deviation with the expectation values in the reference process A ⇌ X band E(t)± (t) (gray and pink) obtained by sampling of 10,000 tra- ( m = 0 ). The transient of the autocatalytic process ( m = 1 ) jectories that were calculated by Gillespie’s simulation method for the is different: It shows substantial broadening of the fluctua- uncatalyzed ( m = 0 ; top plot) and the first-order autocatalytic reac- tion band (Table 1) because of delayed onset of the reaction tion ( m = 1 ; bottom plot). Both reactions approach almost identical thermodynamic equilibrium states. Choice of parameters and initial before it narrows down to the equilibrium value. In the inter- −1 conditions: N = 1000 ; h = h = 1.5[t ], A(0)=1000, X(0)=0 (top mediate range the expectation values differ remarkably from → ← −1 −1 plot); k = k = 0.01[M   t ], A(0)=999, X(0)=1 (bottom). Solu- → ← the deterministic solution curves. As expected an increase tions of the corresponding kinetic differential equations are shown as in X leads to a decrease in the width of the fluctuation broken lines. 0 band. Second order autocatalysis ( m = 2 ) differs from first- order autocatalysis mainly by the size of the characteristic effects: Both autocatalytic reactions show broadening of the E(A)= 4.9951 , for second-order autocatalysis E(A)= 4.9605 fluctuation bands but the band width in the second-order compared to E(A)= 5 for the uncatalyzed process. case is about four times as wide. In particular, the scatter In principle, there are three major sources of random- in the waiting times until the first reaction events occurs is ness in autocatalytic reactions: (i) thermal fluctuations, (ii) much larger because we are dealing with two small factors, delayed onset of autocatalytic reactions and (iii) multiple (X − 1) ⋅ (X − 2) , in the expression for v(0). stationary states. (ad i) In the transients towards equilib- The standard deviation in the course of the reactions, (t) , rium the thermal fluctuation bands are essentially the same is shown in Fig. 3. Because of sharp initial conditions the in autocatalytic and conventional reactions as can be seen fluctuation band starts out from zero— (0)= 0 , increases, best from the comparison of standard deviations at equi- becomes broad in the intermediate range and settles down at librium (Fig.  2). The differences between conventional the equilibrium value (6). Substantial deviations between the and autocatalytic equilibrium densities can be recognized deterministic solution and the stochastic expectation value, numerically for very small particle numbers only (). (ad a(t) and the E A(t) or x(t) and E X(t) , respectively, are ii) The reaction rate for an autocatalytic reaction of order observed in the intermediate range. Accordingly, the stand- m is v(0)= kA(0) X(0) − l X(0) (with m ≥ 1 ). At m m+1 ard deviation goes through a pronounced maximum that is 1 3 European Biophysics Journal (2018) 47:403–425 409 those shown in Fig.  2 shows similarities except a twice as wide range for the values of the stochastic variables, A(0)+ X(0) ≥ A(t) ≥ 0 and 1 ≤ X(t) ≤ A(0)+ X(0) and expectation values and standard deviations approach zero at sufficiently long times. Again, we observe a sig- moid shape of the solution curves, for small initial val- ues ( X(0) < 5 ) the standard deviation (t) becomes large in the intermediate range (Fig. 3), and the deterministic curve deviates substantially from the stochastic expecta- tion values. Autocatalysis in the flow reactor Fig. 3 Comparison of the standard deviation in the revers- Implementation of autocatalytic reactions in the flow reac - ible first-order autocatalytic reaction  +  ⇌ 2 and the revers- tor provides additional insights into the different forms of ible uncatalyzed reaction  ⇌  . The reactions are recorded for the closed system fulfilling A(t)+ X(t)= A + X = const = N where randomness. In particular we are interested in multiple 0 0 A(0)= A and X(0)= X for both reactions. The figure presents 0 0 stationary states as a source of stochasticity (item iii). The the results of statistical evaluation of 10  000 trajectories obtained reaction equations for first order autocatalysis are: by computer simulations with Gillespie’s method (Gillespie 2007) (1) (1) a r for the autocatalytic reaction  (t)=  (t) (black) and for the (0) (7a) (0) ∗ ⟶ , uncatalyzed reaction  (t)=  (t) (red). Choice of parameters: −1 −1 −1 k = k = 0.01[M   t ] and h = h = 1.5[t ]; equilibr ium → ← → ← parameters K = k ∕k = h ∕h = 1 ; initial conditions: N = 1000 , → ← → ← (7b) +  ⟶ 2, X = 1 , A = 999 and N = 1000 , X = 0 , A = 1000 yielding the 0 0 0 0 numerical equilibrium values A = X = 500 . The equilibrium value of the standard deviation is practically the same for both reactions: (7c) ⟶ ∅, and (1) (0) lim  (t)≈ lim  (t)≈ 15.8114 t→∞ t→∞ X X (7d) ⟶ ∅, with the kinetic differential equations qualitatively different from the shallow maximum observed with the standard deviation of conventional chemical pro- da dx =− k ax +(a − a) r and = x (ka − r). cesses see Eq. (4h) . (7e) dt dt The irreversible processes,  + m → (m + 1) The simple first-order autocatalytic process in the flow reac- obtained from Eq. (5) by putting g = 0 , are illustrative tor sustains two long time states: (i) The state of extinction S because they are closer to biology where replication or with lim a(t)= a = a and lim x(t)= x = 0 , and (ii) reproduction occur always under the conditions of irre- t→∞ 0 t→∞ the quasistationary state S with lim a(t)= a = r∕k and versibility. The shape of the solution curves compared to 1 t→∞ Table 1 Fluctuations in the two Autocatalysis Reaction Initial conditions Limit autocatalytic reactions A+ X order ⇌ 2 X and A+2 X ⇌ 3 X X = 1 X = 2 X = 5 X = 10 t → ∞ 0 0 0 0 Zero A ⇌ X – – – – 15.8 First A+ X ⇌ 2 X 123.2 89.8 59.6 44.4 15.8 Second A+2 X ⇌ 3 X – 245.1 226.7 189.4 15.8 The table presents maximal standard deviations  computed from 1 000 trajectories of the two autocata- max lytic reactions for different initial conditions X = X(0) . For the autocatalytic reactions the standard devia- tion (t) passes through the maximum  listed here whereas for the uncatalyzed process it increases max monotonously from (0)= 0 to the equilibrium value (see Fig. 3) The equilibrium fluctuations calculated from equations  () are practically the same for all three reactions. −1 −1 −1 Choice of parameters, N = 1000 , A ⇌ X: h = h = 1.5[t ]; A+ X ⇌ 2 X: k = k = 0.01[M  t ]; A+2 → ← → ← −2 −1 X ⇌ 3 X: l = l = 0.00001[M  t ] → ← 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) The accurate value obtained from the stationary master equation is  = 15.8114 1 3 410 European Biophysics Journal (2018) 47:403–425 determine the probabilities to end up here or there. For suf- ficiently large population sizes, the long-time expectation values of the stochastic variables can be well approximated by linear combination of the deterministic values: (0) (1) (0) (1) E A =  a +  a andE X =  x +  x 0 1 0 1 with  = N ∕(N + N ) and  = N ∕(N + N ) where N 0 0 0 1 1 1 0 1 0 and N are the counted numbers of trajectories ending up in S or S from a sufficiently large sample. Although only 0 1 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, 1982, 1982). 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 extinc- tion ( s = 919). Fig. 4 The first-order autocatalytic reaction.  +  ⇌ 2 in the flow The major difference between the two classes of auto- reactor. The two plots show a single trajectory (top plot) and statis- catalytic reactions lies in the repertoire of possible dynamic tics consisting of expectation value within the one-standard deviation band, E A(t) ±  A(t) and E X(t) ±  X(t) taken form sample of behaviors. First-order autocatalysis gives rise to exponen- 100 trajectories (bottom plot). The four phases of the approach to the tial growth in unconstrained systems and to selection and long-time solution are indicated (top plot; see text). Choice of param- optimization of mean fitness in multispecies cases with −1 −1 eters and initial conditions: N = 2000 ; k = 0.01[M   t ], r = 0.5 −1 finite resources (see ’’Competition, mutation and qua- [V  t ], A(0)= 0 , X(0)= 1 (top) and X(0)= 3 (bottom). Pseudoran- dom number generator: ExtendedCA, Mathematica, seeds s = 491 . sispecies"). Accordingly first-order autocatalysis leads to a The sample size of the bottom plot was 100 trajectories, 3 led to S Darwinian scenario of selection of the fittest. In contrast (extinction) and 97 reached the pseudo-stationary state S . Solu- to first-order autocatalytic reaction networks, the dynamics tions of the corresponding kinetic differential equations are shown as in second-order systems is very rich and includes multiple dashed lines 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 lim x(t)= x = a − r∕k . Stability of S requires that the t→∞ 0 1 models, for example the Brusselator (Nicolis and Prigogine condition r < a k is fulfilled. 1977). It provides a simple enough reaction step for studies Starting from an empty reactor containing no A and by means of rigorous mathematics. Qualitative analysis of the autocatalyst only in seeding quantities, A(0)= 0 and Brusselator dynamics is straightforward and numerical inte- X(0)= 1, 2, 3, … , the trajectory shown in the upper part of gration causes no problem provided the integration software Fig. 4 allows for the distinction of several phases: (I) the can handle stiff differential equations. In reality, however, flow reactor is filled with the resource A in phase I, (II) in single-step autocatalytic reactions are extremely rare, instead the random phase II the decision is made to which state we are commonly dealing with multistep-reaction networks the trajectory will converge, (III) the trajectory approaches (Noyes et al. 1972) (see also the review by Francesc Sagués the long-time state, and (IV) the trajectory fluctuates and Irving Epstein (2003). around the state (see Figs. 4 and 7). First-order autoca- (0) In biology, in particular in the theory of evolution, the talysis sustains the two long-time solutions S =(a ,0) (1) (1) process A+2 X → 3 X plays a special role since in the sim- and S =(a , x ) , stochastic trajectories approach either ple form of hypercycles it is the basis for suppression of of the two states, and parameters and initial conditions 1 3 European Biophysics Journal (2018) 47:403–425 411 competitive selection without destroying template-induced reproduction. It provides one fundamental mechanism for major transitions (Maynard Smith and Szathmáry 1995; Schuster 1996) and will be discussed extensively in "Compe- tition and cooperation". The enormous flexibility of second- order autocatalysis follows, for example, from Fisher’s selec- tion 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 par- ents 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". Fig. 5 A minimal model for modeling evolution. Evolution is con- sidered 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 A minimal mathematical model parameters of all three processes, (i) fitness parameters f correspond- for the evolution of molecules ing to reaction rates for competition, (ii) reaction rates h for catalyzed reproduction, and (iii) an error rate parameter p for mutation are plot- ted on the axes of a Cartesian coordinate system. On the three two- The minimal model is dealing with the time dependence dimensional faces of the coordinate system we are dealing with the of the distribution of genotypes in populations Π(t) and three fundamental evolutionary processes: ( ) competitive reproduc- hence it is rooted in chemical kinetics of (i) competitive tion and mutation are the basis of Darwinian optimization through reproduction, (ii) symbiontic cooperation, and (iii) genetic natural selection and give rise to the formation of quasispecies and eventually to the occurrence of error thresholds (Eigen 1971; Eigen variation. The quantification of these three properties yields and Schuster 1977;  ) the interplay of competition and cooperation three parameters or sets of parameters, which can be plot- allows for the description of major transitions, which are seen as the ted on the three axes of a Cartesian coordinate system (see consequences of crossing resource thresholds (Maynard Smith and Fig. 5 and the previous publications (Schuster 2016a, b, d, Szathmáry 1995; Szathmáry and Maynard Smith 1995; Schuster 1996, 2016); ( ) the combination of cooperation and mutation ena- 2017b). We consider the simplest case here, where the three bles reintroduction of extinguished species provided the error rate is quantities are a fitness parameter f , a cooperation parameter sufficiently large such that a mutation threshold for stochastic survival h, and an mutation rate parameter p. For an implementa- can be recognized (Schuster 20160 tion 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 l Q i ji of a (hypothetical) compound A, A or a , respectively, and +  +  ⟶  +  +  , i, j = 1,… , n ; i mod n ; 0 0 i i+1 i j i+1 the mean residence time  = V∕r of the reaction mixture (8c) in the reactor where V is the volume of the reactor and r is (8d) ⟶ ∅ ; and the (volumetric) flow rate. The parameter  defines the time resolution of the reactor since slow reactions, which do not progress appreciably during the time interval  cannot be ⟶ ∅ , i = 1,… , n. (8e) R i studied. The process (8a) supplies the material required for repro- In the next step, the model is implemented by means of duction. A solution with A at concentration a flows into a a suitable and simple reaction mechanism. Based on "Ideal- continuously stirred tank reactor (CSTR) with a flow rate ized environments" and "Basic processes" we consider the 2 parameter r (Schmidt 2005, p. 87ff). The reactor operates at following set of 2n + n + 2 reactions in the flow reactor: constant volume and this implies that the volume of solu- a r (8a) ∗ ⟶  ; tion flowing into the reactor per time unit [ t] is compensated exactly by an outflow, which is described by the Eqs. (8d) k Q i ji and (8e) and concerns all (n + 1) molecular species, A and (8b) +  ⟶  +  , i, j = 1,… , n ; i i j , i = 1,… , n . Inflow and outflow are often characterized 1 3 412 European Biophysics Journal (2018) 47:403–425 as pseudoreactions because they are no chemical reactions long as n is not too large. In absence of cooperation terms, in the strict sense, which are converting reactants are into l = 0∀ i = 1,… , n , Eq.  (9) can be transformed into an products. The two classes of reactions producing progeny, eigenvalue problem of a symmetric matrix, which is readily template induced replication (8b) and catalyzed template diagonalized provided n is not very large ( n < 10 ; "Com- induced replication (8c), represent the core of the evolu- petition, mutation and quasispecies"). In mutation-free tion model. In agreement with the conditions in biology, systems, p = 0 ("Competition and cooperation"), qualitative both reproduction steps are implemented irreversibly in the analysis and determination of stationary states are straight- direction of polynucleotide synthesis. A basic assumption forward, and the dynamics of the complete system can be for both reproduction steps is that correct reproduction and derived by extrapolation from the error-free results to finite mutation are parallel chemical reaction channels ("Basic mutations rates. The cooperation system with mutation, (9) processes"). In other words, there is no mutation under con- with k = 0∀ i = 1,… , n is used here to study the relevance ditions that do not sustain reproduction. of mutation in symbiontic systems. It also serves as an exam- As an alternative to the Eigen model (1) mutation can ple for the study of unconventional consequences of repli- be seen, for example, as the result of DNA damage and cation with frequent errors in the strong mutation scenario imperfect damage repair during the whole life span of an ("Cooperation and mutation"). organism, which is the idea underlying the Crow–Kimura mutation model (Crow and Kimura 2009, pp  264–266). Processes on individual coordinate axes Then reproduction and mutation are completely independ- ent processes, The processes along the individual coordinate axes are considered in order to verify the initial statement on ji the three basic processes. For this goal it is easiest to set +  ⟶ 2 and  ⟶  , (8b’) i i i j certain parameters zero: (I) no natural selection implies and in the kinetic differential equations they appear as addi- k = k =…= k = k (= 0) , (II) no cooperative coupling 1 2 n tive terms. Interestingly, the Eigen and the Crow–Kimura requires l = l =…= l = l = 0 , and (III) no mutation 1 2 n model although being fundamentally different with respect leads to Q =  where  is the unit matrix. to the assumptions about the nature of the mutation pro- A process taking place on the selection axis is given by cess give rise to the same mathematical problems [see, e.g., (II) and (III) being true leads to the ODE (Schuster 2016d, pp.76-78)]. The mutation matrix Q cor- responds formally to the mutation matrix  but there are da =−a k x + r (a − a) and (10a) i i 0 non-negligible differences Q covers correct and error-prone dt i=1 replication but the process  → 2 is handled separately i i in the Crow–Kimura model and hence all diagonal terms are dx zero μ = 0∀ i = 1,… , n. ii =(k a − r) x ; i = 1, 2, … , n . (10b) i i dt The equations that will be applied in the analysis of the dynamics of the model (8) implement three processes along Competition between the reproducing elements the coordinate axes: (i) Darwinian selection of the fittest leads to survival of the fittest subspecies  , which is on the competition axis, (ii) hypercycle dynamics on the defined by k a = f = max{f ;i = 1, … , n} . For constant m m i cooperation axis, and (iii) neutral evolution on the mutation []= a = const , the mean fitness of the population, ∑ ∑ n n axis. The kinetic differential equations of the model mecha- (t)= f x (t) c(t) with c(t)= x (t) , is a non- i i i i=1 i=1 nism (8) are of the form: decreasing function of time: d∕ dt = var{f } ≥ 0 . This result is the formal analogue of Fisher’s fundamental theorem for da asexual reproduction. =−a k + l x x + r (a − a) , i mod n and i i i+1 i 0 dt i=1 (9a) We remark that the deterministic kinetic Eqs. for (8b) and (8c) dx were extensively studied under the simplifying assumption of = a Q k + l x x − x r , n ij j j j+1 j i constant population size x (t)= c = const (Eigen 1971; i 0 dt i=1 (9b) j=1 Eigen and Schuster 1978; Eigen and McCaskill 1989). The deter- ministic solution curves formulated in relative concentrations i = 1, 2,… , n; j mod n . ξ (t)= x (t) x (t) are identical for the CSTR and for constant i i i i=1 population size (Schuster and Sigmund 1985). As we mentioned In  (9a) we made use of the conservation relation already in "Idealized environments" the stochastic system is unstable Q = 1 . No analytical solutions are available for  (9) for unregulated constant population size (Jones and Leung 1981) and ij i=1 a proper description has to account explicitly for the physical setup in general but numerical integration is straightforward as applied, here the flow reactor. 1 3 European Biophysics Journal (2018) 47:403–425 413 On the cooperation axis, the conditions (I) and (III) are assumption that chemical processes occur one at a time, all fulfilled and we obtain the equations of hypercycle dynamics jumps involve single steps and the particle numbers change by ±1 unless the elementary steps involve more than a single da molecule per class. The jumps S  → S or S → S  are =−a l x x + r (a − a) and (11a) i i i+1 0 denoted by the shorthand notation dt i=1 =(m ± 1, s , … , s ) ≡ (; m ± 1) or 1 n dx i � =(l x a − r) x , i = 1, 2, … , n, i mod n ,  =(m, s , … , s ± 1, … , s ) ≡ (; s ± 1). (11b) i i+1 i 1 k n k dt Then the master equation of mechanism (8) takes on the which were studied in great detail in a number of previous form papers (Eigen and Schuster 1978; Schuster 1978; Schuster et al. 1979, 1980; Hofbauer et al. 1991). dP The third case—with (I) and (II) being true—yields a = a r P − P + r (m + 1)P − mP 0 (;m−1)  (;m+1) dt degenerate deterministic solution: All distributions of sub- species with fixed population size x = c yield the same i=1 + r (s + 1)P − s P i (;s +1) i solutions and therefore constitute an (n − 1)-dimensional i=1 invariant manifold fulfilling + Q (k + l s ) (m + 1)(s − 1)P − ms P ii i i i+1 i (;m+1,s −1) i da i=1 =−ak c + r (a − a) and (12a) n n dt + Q (k + l s )s (m + 1)P − mP . ij j j j+1 j (;m+1,s −1) i=1 j=1,j≠i dc = ak c − cr =(ak − r) c and (14) dt (12b) =(a − ) r with  = a + c dt Each reaction step involving S changes the probability to be in state S , P , in two ways: It increases the prob- In the absence of selection and cooperation neutral evolu- ability through reactions or pseudoreactions S → S tion in the sense of Motoo Kimura (1955, 1983) is observed and decreases the probability through the reaction steps on the mutation axis. For an adequate description of the S → S where  is summed over all states from which process a stochastic treatment is required. Random selection S can be reached or vice versa. The two terms in the first of a single arbitrarily chosen subspecies is observed in the line, for example, describe the two pseudoreactions mod- no-mutation limit, lim p → 0 . For non-vanishing mutation eling inflow and outflow of the material A , and further each rates once selected subspecies may be replaced by other sub- reaction is represented by two steps. It is also worth noticing species and the mean time of replacement of one randomly that stoichiometry requires two slightly different replication −1 selected subspecies is T =  where  is the mutation rep terms depending on whether the copy is correct or incorrect. rate per generation (Kimura 1955, 1983). Translated into our Master equations are easily written down and station- model, we find for single point mutations  ≈ p. ary solutions can be derived by generally applicable techniques as it was shown for the flow equilibrium in Master equation and simulation "Deterministic and stochastic approaches" but explicit time-dependent solutions are very hard to obtain and Reaction  (8) can be cast into chemical master equations. known only in exceptionally simple cases [Schuster 2016e, The particle numbers of the molecular species, []= A(t) pp 347–568]. Here, we shall use the simulation technique and [ ]= X (t) with i = 1,… , n , are integers and in the i i of sampling trajectories introduced already in "Autocataly- absence of flows they fulfil the conservation relations, sis in the batch reactor". Expectation values and second A(t)+ X (t)= C . The variables of the master equation i=1 moments of variables can be computed through sampling are the probabilities of trajectories—an example is shown in Fig. 6—but often this approach exceeds the available computational facili- P (t)= Prob A(t)= m and m ties and it is necessary to interpret single trajectories. As examples we consider single trajectories in Fig.  4 and (13) P (t)= Prob X (t)= s ; i = 1, … , n , s i i in Fig.  7. The process for convenience starting from an empty flow reactor is split into four phases: (i) establish- and the indices are subsumed in an index vector, ment of the flow equilibrium of A , (ii) random decision =(m, s , … , s ) , which characterizes the state S of 1 n on the (quasi)stationary state towards which the trajectory the system. The chemical master equation is based on the 1 3 414 European Biophysics Journal (2018) 47:403–425 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 pro- cess is assumed to start with an empty reactor except seeds for the two autocatalysts  and  . It can be partitioned into four phases: 1 2 (I) fast raise in the concentration of A, (II) a random phase where (1) (2) the decision is made into which final state— S , S , S or S —the 0 2 1 1 trajectory progresses, (III) the approach towards the final state— here S —and (IV) fluctuations around the values of the (quasi) stationary state. Comparison with the simple autocatalytic pro- cess in Fig.  4 reveals great similarity. Choice of parameter values: −1 −1 −2 −1 k = 0.099 , k = 0.101  [M t ], l = 0.0050 , l = 0.0045  [M t ], 1 2 1 2 −1 Fig. 6 Quasispecies formation. The two plots show the forma- a = 200 , r = 4.0  [V  t ], pseudorandom number generator: Extend- tion of quasispecies from an initially empty reactor, A(0)= 0 , with edCA, Mathematica, seeds s = 631 . Initial conditions: A(0)= 0 , replicators in seeding amounts. The system in the upper plot was X (0)= X (0)= 1 . Color code: A(t) black, X (t) red, and X (t) yellow. 1 2 1 2 initiated to be a single copy of the master sequence, X (0)= 1 and The figure is adapted from Schuster (2017a) X (0)= X (0)= X (0)= 0 . At the beginning the system might be 2 3 4 extinguished by a single dilution event X → ∅ and the high prob- ability of extinction gives rise to enormously broad and overlapping fittest subspecies and at non-zero mutation rates this sce- one-standard-deviation bands. The initial values in the lower plot nario yields selection of a fittest ensemble of subspecies, were chosen such that the probability of extinction is zero for all prac- tical purposes: X (0)= 10 and X (0)= X (0)= X (0)= 0 . Choice which has been characterized as quasispecies (Eigen and 1 2 3 4 −1 −1 of parameters: N = 2000 ; k = 0.011[M   t ], k = k = 0.010 1 2 3 Schuster 1977; Domingo and Schuster 2016). Precisely, the −1 −1 −1 −1 −1 [M   t ], k = 0.009[M   t ], r = 0.5[V  t ], Color code: A black, quasispecies is the stable stationary long-time distribution red,  yellow,  green, and  blue. Expectation values are 1 2 3 4 of a population of subspecies undergoing replication and shown as full lines, deterministic solutions as broken lines. Since x (t)= x (t) is fulfilled for the parameter values used here the curve is mutation: A population that consists of several genotypes 2 3 shown as a yellow–green broken line present in time-dependent concentrations, Π(t)= x (t) ⊕ x (t) ⊕ … ⊕ x (t) (15) 1 1 2 2 n n has the stationary solution, lim Π(t)=  = x  ⊕ x  ⊕ … ⊕ x  , which is t→∞ 1 1 2 2 n n converges, (iii) the approach towards this (quasi)station- called the quasispecies . ary state and (iv) fluctuations around the (quasi)stationary Deterministic quasispecies The deterministic or continu- state. The separation into phases is made possible by the ous quasispecies represents the unique deterministic long- choice of suitable initial conditions. time solution of the replication mutation problem, which is described in the flow reactor by the ODE (Schuster and Competition, mutation and quasispecies Sigmund 1985): The bottom face of the three-dimensional parameter space— da  in Fig. 5—is dealing with selection provided by the com- =−a k x + r(a − a) (16a) i i 0 dt bination of competition and mutation. Natural selection at i=1 zero mutation rate leads to a homogeneous population of the 1 3 European Biophysics Journal (2018) 47:403–425 415 an approximate uniform distribution. An explanation is dx i straightforward: Above this critical transition, mutations = a Q k x − rx , i = 1,… , n . (16b) ij j j i dt occur too often to sustain sufficiently accurate reproduction j=1 of the template sequence and the result is random replica- tion: In the long run, every sequence is obtained with the All early works on quasispecies dynamics were per- same probability. The transition has been characterized as formed with the constraint of constant population size: error threshold (Eigen 1971; Eigen and Schuster 1977) since x (t)= c = const with a(t)= a and f = k a that gives i 0 i i 0 i=1 evolutionary dynamics does not sustain a structured long- rise to the differential equation time population at higher error rates, p > p . On typical cr dx fitness landscapes, the error threshold sharpens with increas- = Q f x − (t) x , ij j j i ing chain length l and becomes a first-order phase transition dt j=1 in the limit l → ∞ (Tarazona 1992; Park et al. 2010; Huang f x 7 j j j=1 et al. 2017). with  = ; i = 1, … , n , (16’) Discrete quasispecies In the continuous description, the j=1 quasispecies contains all species at finite positive concentra- which fulfils dc∕ dt = dx ∕ dt = 0 .   Equations (16) i=1 tions no matter how small the concentrations might be. Deal- and (16’) have identical solutions in normalized variables ing with less than one molecule per reactor volume, how- (t)= x ∕ x (t) but the stability properties of the corre- i i i i=1 ever, is unrealistic. Numbers of molecules are non-negative sponding stochastic systems are different. [For more details integers, X ∈ ℕ , subspecies distributions are truncated in see (Thompson 1974; Jones et al. 1976; Ebeling and Mahnke reality and we call them stochastic or discrete quasispecies: 1979; Jones and Leung 1981; Schuster and Sigmund 1985)]. The quasispecies as a function of the mutation rate param- ⌊x ⌋ if x ≥ 1 i i = X  ⊕ X  ⊕ … ⊕ X  with X = . 1 1 2 2 n n i eter (p) begins at as a homogeneous population con- 0 if x < 1 p = 0 taining exclusively the fittest subspecies  and becomes a (17) distribution of subspecies at nonzero p. This distribution The discreteness of the stochastic variables leads to a consists of a most frequent sequence called master sequence, modification of the scenario for the mutation dependence which is surrounded by a mutant cloud. Commonly, the most of quasispecies (p) . In the mutation-free system, p = 0 , frequent or master sequence is also the fittest one, but this survival of the fittest is observed and the quasispecies con- is not necessarily so: In case of strong stationary mutation sists of a single sequence:  = X  = N with m indi- m m m flow from the mutant cloud to the master sequence, a less cating maximal fitness, f = max{f ;i = 1, … , n} , a nd N m i fit master may outgrow a fitter sequence with less efficient being the population size. At sufficiently small values of the mutational backflow. At further increase in p the distribution mutation rate parameter p, no subspecies except the master broadens, and eventually ends up as the uniform distribu- exceeds the threshold x ≥ 1 and the discrete quasispecies tion  ∶ { x =  x = ⋯ =  x } at  p = 1 −( − 1) p∕ = 1∕ still consists of a single fittest genotype  . The scenario 1 2 n l l l where n =  and  x = 1∕ (i = 1, …  ) for sequences of in the small mutation rate regime—populations are almost chain length l over an alphabet with  letters. At the muta- always homogeneous except short non-stationary periods tion rate , correct and wrong digits are incorporated p =  p during which advantageous mutations take over the popula- with the same frequency. The uniform distribution is the tion—is tantamount to the strong selection–weak mutation dynamical answer to the absence of any preference for nucle- scenario discussed in evolutionary biology (Gillespie 1983; otide assignment: All subspecies in the stationary distribu- Joyce et al. 2008). If the p-value is so large that for one or tion occur with the same frequency. The quasispecies in the intermediate range is determined by the fitness landscape— We mention that some model landscapes like the additive land- the distribution of fitness values f in sequence space, by scape or the multiplicative landscape sustain smooth rather than sharp transitions (Wiehe 1997). For a recent update and a review of the the move set of allowed mutations as well as the mutation state of the art in the relation between fitness landscape and selection rate p. Typically, sharp transitions occur at some critical dynamics see (Schuster 2016d). mutation rates p = p : The distribution changes smoothly tr Biologists (Dean and Thronton 2007; Sniegowski and Gerrish from to , where the distribution turns abruptly p = 0 p = p tr 2010) and computer scientists commonly distinguish strong and weak into another distribution. At the transition with the larg- mutation. The weak mutation scenario assumes that adaptive muta- tions are sufficiently rare and do not interfere with the selection pro- est p-value , the quasispecies becomes p = 1∕𝜅> p > p cr tr cess but initiate replacements of currently fittest genotypes by still fitter variants. The strong mutation scenario is characterized by suf- ficiently large values of p that give rise to the quasispecies dynam- ics described here [for details see (Domingo and Schuster 2016)] and The value  = 2 is used here and applies for binary sequences. For to mutation induced cooperative dynamics ("Cooperation and muta- the natural AUGC -nucleotide alphabet we have  = 4. tion"). 1 3 416 European Biophysics Journal (2018) 47:403–425 more subspecies x ≥ 1 is fulfilled, selection of a family of Fluctuations at small particle numbers have different ori- genotypes is observed:  consist of a master genotype  , gins ("Autocatalysis in the batch reactor"): (i) conventional together with a stationary distribution of sufficiently frequent thermal fluctuations, (ii) enhanced fluctuations related to mutants  (i ≠ m ). The quasispecies becomes broader with autocatalytic self-enhancement, and (iii) anomalous fluctua- increasing mutation parameter p until a threshold value p tions in the stochastic variables arising from two or more cr is reached above which error propagation does not sustain quasistationary states. The standard deviations (t) fulfil the a stationary state and the population Π drifts randomly N-law for the resource A but are larger for the autocata- through sequence space (Huynen et al. 1996; Higgs and Der- lysts  ,  ,  , and  . The stationary states of the stochas- 1 2 3 4 rida 1991). Interestingly, the critical mutation rate p can be tic system are extinction S and quasispecies S . Fig. 6 shows cr 0 1 derived from the continuous quasispecies theory. a typical example: The anomalous fluctuations in the upper As calculations of the time dependencies of the first two plot are in full analogy to first-order autocatalysis (Fig.  4). moments of the probability distribution of subspecies in the Since the initial condition X (0)= 1 was chosen, one out- (Π) population, P (t) , reveal (Jones and Leung 1981)  (16’) is flow step may extinguish the population and the probabil- only marginally stable: ity of dying out is indeed as large as 20%. The fluctuations bands are extremely broad, and large differences between � � d⟨N⟩ the deterministic solution and the corresponding expecta- = f ⟨X ⟩ − ⟨ (t) N⟩ = 0, (18a) i i tion values are observed. In the lower plot initial conditions dt i=1 X (0)= 10 were chosen, which are sufficient to reduce the contributions of anomalous fluctuations practically to zero. � � d⟨N ⟩ Then, for a population size of N = 2000 the concentration = 2 f ⟨X N⟩ − 2⟨ (t) N ⟩ + i i dt a(t) coincides with the expectation value E A(t) for all prac- i=1 (18b) n n tical purposes and the fluctuations fall into a typical ± N � � d⟨N⟩ + 2 f ⟨X ⟩ − = 2 f ⟨X ⟩ , and -band. The autocatalysts, X ,… , X , show broader than i i i i 1 4 dt i=1 i=1 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 2 2 var(N)= ⟨N ⟩ − ⟨N⟩ = 2 f ⟨X ()⟩ d (18c) i i an open system and does not approach an equilibrium state i=1 (see also Fig. 4). The time derivative of the first moment vanishes as expected It is worth recalling what means stochasticity for qua- since the condition of a constant population size is fulfilled sispecies: (i) continuous concentrations are replaced by by the differential equation (16’). The variance, however, discrete particle numbers, (ii) fluctuations replace single increases with time since the integrand in (18b) is always line trajectories by bands within which trajectories follow a positive. After sufficiently long time, the fluctuation band probability distribution, (iii) subspecies can be diluted out becomes so large that the expectation value is irrelevant for of the flow reactor and if this happens for all subspecies the the description of the system. The instability in Eigen’s qua- population goes extinct giving rise to anomalous fluctua- sispecies equation (Eigen 1971) is well known from similar tions, and (iv) error thresholds introduce random reproduc- problems: (i) the neutral birth-and-death process with equal tion that is closely related to Motoo Kimura’s random drift. birth and death parameters,  =  , and (ii) the Wiener pro- An increase in the error rate up to the error threshold leads to cess. In both cases a constant expectation value is jeopard- broadening of the mutant spectrum surrounding the master ized by a variance that grows linearly with time. Stability sequence. Above the thresholds, the populations migrate by against fluctuation is easily introduced into (16’): one needs random drift through sequence space. only to give up the condition of strictly constant population size and to replace the denominator in (t) by a constant Θ, Competition and cooperation The kinetic equations for replication describing the template (t)= f x Θ , j j induced, uncatalyzed and catalyzed processes are obtained j=1 from Eq. () by neglect of mutation. Then the kinetic differ - ential equations for competition and cooperation of competi- where Θ is the population size towards which the population tors result from (9) by setting Q =  : converges after sufficiently long time. The approach cho- sen here, the implementation of a flow reactor as a physical da device, yields a stable system as well. =−a k + l x x + r (a − a) and (19a) i i i+1 i 0 dt i=1 1 3 European Biophysics Journal (2018) 47:403–425 417 Table 2 Asymptotically stable Symbol Stationary values Stability range stationary states of Eq. (19) with n = 3Schuster (2016) a x x x 1 2 3 S a 0 0 0 0 ≤ a ≤ 0 0 0 3 (3) 0 0 a −   ≤ a ≤  + S 3 0 3 3 0 3 32 (1) 0 a −  −    +  ≤ a ≤  +  + 3 0 3 32 32 3 32 0 3 32 31 S      +  +  ≤ a 3 3 1 2 3 32 31 0 The four stationary states are ordered with respect to increasing a -values of their asymptotically stable regime. The relations k < k < k and l > l > l between the rate parameters were assumed. The abbre- 1 2 3 1 2 3 viations  = r∕k ,  =(k − k )∕l ) , and  =(k − k )∕l ) are used for the combination of param- 3 3 31 3 1 1 32 3 2 2 eters. For the cooperative state S the stationary concentration of A is obtained as one root of a quad- ∑ ∑ 3 3 ratic equation (20) with two combinations of the rate parameters,  = k ∕l and  = 1∕l , and i i i i=1 i=1 =(r − k )∕(l ) for i = 1, 2, 3 and  = a from   (20). The existence of the non-trivial stationary state i i i 1 requires a sufficiently small flow rate: r ≤ (a + ) ∕4 a and is stable for a < r∕k , the state of selection of  , 0 0 3 3 dx i (3) S , is stable in the range r∕k < a < r∕k +(k − k )∕l , = a k + l x − r x ; i = 1… , n; i mod n . 3 0 3 3 2 2 i i i+1 i 1 dt (1) followed by ex c l u s i o n of  , S , in the range (19b) k +(k − k )∕l < a < k +(k − k )∕l +(k − k )∕l . 3 3 2 2 0 3 3 2 2 3 1 1 Both equations contain terms of different molecularity in Above this value for a cooperation of all three sub- and this has the consequence that the dynamic behav- ∑ species is observed provided the flow rate is not too ior depends strongly on the total concentration c = x , i 2 i=1 large: r < r =(a + 𝜓 ) ∕4𝜙 with  = k ∕l and cr 0 i i i=1 which in turn is determined by the amount of available = 1∕l . The free concentration of A is obtained as i=1 resources, a : At sufficiently low concentration, the first- 0 9 solution of a quadratic equation order terms, k x , dominate whereas the second-order terms, i i l x x become important at high concentrations. No direct 1 i i+1 i a = a +  ∓ (a + ) − 4r , (20) 1,2 0 0 analytical solutions are available but exhaustive qualitative analysis is possible and the concentrations of the stationary where the minus sign corresponds to the cooperative state S . solutions, a, x (i = 1, … , n) , can be computed analytically The second solution belongs to an unstable state S , which Schuster (2016a, d) (Table 2). 3 separates the basins of attraction of the states of coopera- The basis of the calculations of stationary states tion and extinction. Above the critical flow rate, r > r the cr is the factorability of the r.h.s. of (19b). The relation states S and S do not exist. For n > 3 , the situation becomes x ⋅ (k + l x )a − r = 0 with i mod n is compatible with 3 i i i i+1 more complex since solutions may oscillate. Many systems stationary states that fulfil either with n = 4 have oscillations with very weak damping fac- 1 r (i) (ii) tors, n ≥ 5 commonly leads to undamped oscillations. The x = 0 or x = − k , ∀ i = 1, … , n , i−1 i i l a i−1 properties of these systems have been discussed extensively n in previous publications to which we refer here (Schuster and the conservation condition a + x = a . In total i 0 i=1 2016e; Schuster and Sigmund 1985; Schuster 1978). there are 2 possible states out of which a single one is Stochasticity has common effects on competition–coop- asymptotically stable for a given set of parameters. In eration systems like thermal fluctuations and fluctuation Table 2 the results for n = 3 are shown. The number of sub- through autocatalytic enhancement. In addition, there are species present at the stationary state, n , is suitable for char- many more quasistationary states than the asymptotically acterizing the state: n = 0 is the state of extinction, n = 1 S S stable states of the deterministic system. For example, states is the state of selection, n = 2 is the state of exclusion, in which less efficient subspecies are selected show up as etc., and finally n = 3 = n occurs at the state of coopera- quasistationary states as well. In case of the smallest pos- tion where all three subspecies are present. The numbers of sible system with n = 2 and k > k , we have four states: 2 1 long-time subspecies n depend on the resource input a . S 0 (i) the absorbing boundary as state of extinction S , (ii) the The relative size ordering of the parameters determines the (2) state of natural selection S where the fittest variant  is identity of the selected, excluded, etc., subspecies. For the calculations shown in Table 2 the orderings k < k < k and 1 2 3 l > l > l were applied and hence  is selected and  is 1 2 3 3 1 It is worth mentioning that Eq.  (20) and the solutions for a are 1,2 excluded. Considering steady states as functions of the avail- valid for arbitrary n provided the summations in  and  are adjusted able resources, the state of extinction S comes first at small accordingly. 1 3 418 European Biophysics Journal (2018) 47:403–425 Table 3 Probabilities to reach Initial values Counted numbers of states in final outcomes quasistationary states in the cooperative regime with n = 2 X (0) X (0) N N (1) N (2) N 1 2 S S S S 0 2 1 1 and different initial conditions 1 1 385.1 ± 23.6 1481.0 ± 36.8 1719.6 ± 37.8 6414.3 ± 53.8 2 1 77.4 ± 9.1 1822.6 ± 41.6 367.6 ± 17.0 7733.3 ± 38.3 1 2 71.6 ± 8.5 280.6 ± 20.0 2075.8 ± 28.9 7572.0 ± 39.2 3 1 15.0 ± 2.9 1900.4 ± 30.9 74.69 ± 10.0 8009.0 ± 35.3 1 3 14.0 ± 3.7 53.1 ± 4.8 2180.5 ± 48.4 7752.3 ± 53.8 2 2 14.9 ± 2.6 303.7 ± 16.0 354.5 ± 23.8 9326.8 ± 44.9 3 3 0 70.2 ± 10.0 106.2 ± 10.9 9823.4 ± 15.7 4 4 0 12.1 ± 2.6 28.0 ± 5.0 9959.9 ± 6.4 5 5 0 2.5 ± 1.1 6.3 ± 2.6 9991.2 ± 3.0 The table provides probabilities of occurrence for all four possible long-term states: extinction S , selec- (1) (2) tion of  in S , selection of  in S , or cooperation S . The counted numbers of events are sample 1 2 2 1 1 means and unbiased standard deviations calculated from ten packages, each of them containing 10 000 tra- jectories 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, Mathemat- −1 −1 −1 −1 −2 −1 −2 −1 ica). Choice of parameters: k = 0.09[M t ], k = 0.11[M t ], l = 0.0011[M t ], l = 0.0009[M t ], 1 2 1 2 −1 −4 a = 200 , r = 0.5[V t ]. Initial value A(0)= 0 . Probabilities are obtained by multiplication by 10 Table 4 Long-time behavior Mutation Initial values Expectation values at t = 30 Counts in the cooperation-mutation system with n = 3 and different p X (0) X (0) X (0) E(A) E(X ) E(X ) E(X ) P(S ) 1 2 3 1 2 3 3 initial conditions 0 1 1 1 279.3 44.56 36.06 39.58 0.301 0 2 2 2 81.53 117.7 94.9 105.4 0.814 0 4 4 4 2.023 146.1 119.6 132.2 0.996 0 10 10 10 0.376 147.0 120.0 132.6 1.0 0 Deterministic 0.377 147.0 120.3 132.3 1 0.05 1 1 1 177.2 81.66 68.16 72.64 0.432 0.05 2 2 2 28.57 136.1 113.8 121.5 0.937 0.05 4 4 4 0.383 147.0 122.2 130.4 1.0 0.05 Deterministic 0.377 146.7 122.3 130.5 1 0.1 1 1 1 161.6 86.95 74.39 76.96 0.599 0.1 2 2 2 28.57 136.1 113.8 121.5 0.957 0.1 4 4 4 0.383 147.0 122.2 130.4 1.0 0.1 Deterministic 0.377 146.7 122.3 130.5 1 The table provides expectation values at the time of the end of the simulation (t = 30 ) for different end mutation rate parameters p and different initial conditions. Choice of parameters: l = 0.011 , l = 0.010 , 1 2 −2 −1 −1 l = 0.009[M t ], a = 400 , r = 0.5[V t ]. Initial value A(0)= 0 3 0 selected, (iii) the state of selection of the less fit subspecies solution. In case of n = 3 we present expectation values of , and eventually (iv) the state of cooperation S . The rela- the four stochastic variables, A(t) and X (t) ( i = 1, 2, 3 ) at 1 2 i tive stabilities of the individual states are reflected by the some predefined time of the simulation end, t for differ - end probabilities to reach these states by randomly chosen trajec- ent initial conditions and mutation rate parameters p. The tories (Table 3). Parameters for the calculations shown in the values for p = 0.0 refer to the pure competition–coopera- table were chosen such that the corresponding deterministic tion case discussed here (Table 4). As expected extinction system is situated in the cooperative domain in parameter plays a major role at small initial values, X (0)= 1, 2, 3 space, and indeed the approach towards the cooperative state ( i = 1, 2, 3 ), but X (0)= 4 is already sufficient for coming S has always the largest probability. Again, initial particle close to the expectation values obtained for large initial numbers around X (0)+ X (0)= 4 are sufficient for strong values, X (0)= 10 is enough for reaching the deterministic 1 2 i dominance of the state corresponding to the deterministic values for practical purposes. 1 3 European Biophysics Journal (2018) 47:403–425 419 In systems with n ≥ 4 deterministic and stochastic solu- tion curves oscillate. The solutions of the ODE’s are differ - ent for n = 4 where weakly damped oscillations occur and for n ≥ 5 showing undamped relaxation oscillations (Phil- lipson 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 domi- nant 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 param- eter p was chosen and where accordingly quasispecies for- mation 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 diver- Fig. 8 Times to extinction as a function of available resources in the sity through variation. In principle, mutation can reintroduce five membered cooperative system ( n = 5 ). Extinction times T of the extinguished subspecies into the population. Here, we shall populations Π are shown for different amounts of available resources focus on this aspect and, in particular, study the influence of measured as inflow concentrations a or A when expressed in num- 0 0 the mutation rate parameter p on extinction times. To study bers of molecules per unit volume. The upper diagram presents the data at a resolution of ten molecules (Δ A = 10 ; 100, 110, 120, ) for the role of mutation in low-dimensional cooperative systems 0 four different values of the mutation rate parameter: p = 0.0000 (red), ( n = 2, 3 ) expectation values of the stochastic variables A(t), p = 0.0005 (yellow), p = 0.0010 (green), and p = 0.0020 (blue). T X (t) and X (t) , or X (t) and X (t) were calculated at prede- 1 2 2 3 -values above 1000 are truncated at this value. The lower diagram fined times ( t ) and compared with the probabilities of tra- shows the two plots p = 0.0000 (red) and p = 0.0020 (blue) at the end highest possible resolution (Δ A = 1 ). Choice of other parameters: jectories to end up in the cooperative state, S or S , respec- 0 2 3 −1 −1 −2 −1 r = 0.5[V t ], l = l = l = l = l = 0.01[M t ]. Pseudorandom 1 2 3 4 5 tively. For the initial conditions X (0)= X (0) = X (0) > 4 1 2 3 number generator: Extended CA, Mathematica, seed: s = 491 . Initial the results are practically indistinguishable from the deter- conditions: A(0)= 0 , X (0)= X (0)= X (0)= X (0)= X (0)= 5 1 2 3 4 5 ministic values. An increase in the mutation rate parameter p shows the expected inu fl ence: Extinguished subspecies can be reintroduced and this increases the probability of reaching are less time consuming and provide in essence the same the cooperative state, S or S . The case n = 3 is shown as an results. The plots shown in Fig. 8 show enormous scatter but, 2 3 example in Table 4. nevertheless, allow for drawing two conclusions: (i) In the The oscillating systems are more difficult to investigate. mutation-free case the extinction time T is independent of Here, we consider the time of extinction of the entire popu- the amount of available resources up to a value A ≈ 1300 , lation as a function of the mutation rate and the available and (ii) for non-zero mutation rates a kind of noisy or sto- resources, T (p, A ) . The results are shown in Fig  8: The 0 0 chastic threshold phenomenon is observed. Considering the extinction times T show very strong scatter and their appear- noisy function T (p, A ) and taking A at the first value of 0 0 0 ance is dependent on the resolution of the calculations. By T (p, A ) ≥ 1000 we find for the parameter values applied in 0 0 (T ≥1000) (cr) resolution, we mean here the number of molecules A in A Fig. 8: A = A = 1130, 690, 360 for the mutation 0 0 between two neighboring points. The highest resolution is rates p = 0.0005, 0.0010, 0.0020 , respectively. As expected (cr) achieved when the calculations are performed with every the threshold moves to lower A -values with increasing (integer) number of A molecules, e.g., ΔA = 1 yields 100, 0 0 mutation rate. The behavior of the extinction times T (A ) is 0 0 101, 102,  . Computations at somewhat lower resolutions 1 3 420 European Biophysics Journal (2018) 47:403–425 (cr) similar for n = 4 but the critical concentrations for the different p-values lie much closer together and the analysis is more difficult. Considering survival at constant resources A 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 rever- tant 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 (Albery and 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 strat- egy 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 coop- erative reproduction, or with the introduction of cooperation into Darwinian systems. The characteristic dependence of the population dynam- ics 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 dis- cussed in "Competition and cooperation". An increase of the Fig. 9 The transitions from competition to cooperation in the sys- cooperation parameter h = lA(t) leads in steps from selec- tem with n = 5 The transitions from competition to cooperation in the system with n = 5 . The three plots show single trajectories for tion of the fittest to a cooperative state with all subspecies three different scenarios in the flow reactor: (i, topmost plot) the qua- present. Oscillating systems ( n ≥ 4 ) are more spectacular sispecies scenario with a dominant master sequence (  ), (ii, mid- since the hypercycles are unstable at p = 0 in the stochas- dle plot) an intermediate scenario with irregular dynamics and two tic approach and raising the cooperation parameter h leads dominating species (  and  ), and (iii, bottom plot) the stochastic 1 5 hypercycle scenario with irregular, undamped oscillations. Choice from selection to extinction. In the intermediate param- −1 of parameter: k = 0.150 , k = k = 0.125 , k = k = 0.100 [M 1 2 5 3 4 eter range where the deterministic system shows stepwise −1 t ], l = l = l = l = l = l , l = 0.0 (topmost plot), l = 0.002 1 2 3 4 5 −2 −1 increase in the number of coexisting subspecies (1, 2,… , n (middle plot), l = 0.01  [M t ] (bottom plot), a = 800 , −1 or, expressed phenomenologically, selection, exclusion, , r = 0.5 [V t ], p = 0.075 , pseudorandom number generator: Extend- edCA, Mathematica, seeds s = 491 . Initial conditions: A(0)= 0 , cooperation) the stochastic approach yields highly irregular X (0)= X (0)= X (0)= X (0)= X (0)= 5 . Color code: (t) black, 1 2 3 4 5 dynamics with different numbers of non-extinguished sub- (t) red,  (t) yellow,  (t) green,  (t) blue, and  (t) cyan 1 2 2 2 5 species whereby the number of species present increases 1 3 European Biophysics Journal (2018) 47:403–425 421 with increasing values of the cooperation parameter h. The of asexually reproducing genotypes in a population of con- scenario in parameter space is completed by considering stant size N. Correct replication and mutation are seen as the series of states at different mutation rate parameters p . parallel chemical reactions leading to a uniquely den fi ed sta - The case n = 5 , which for large h-values leads to undamped tionary population called quasispecies (Eigen and Schuster oscillations with a stochastic contribution to the amplitude, 1977). RNA replication catalyzed by single virus specific was chosen to facilitate the illustration (Fig.  9): At zero enzymes from RNA bacteriophages provides a bridge from mutation rate p = 0 the series with increasing h is selection chemistry to biology: The mechanism of the replication pro- → irregular dynamics with two species → irregular dynam- cess is well understood in all molecular details (Biebricher ics with three species → … → extinction. At intermediate 1983; Biebricher et al. 1984, 1985) and an appropriate rep- p-values, we find quasispecies → irregular dynamics → lication assay serves for in vitro evolution studies (Mills oscillations with highly irregular spacings. At high mutation et al. 1967; Biebricher 1983). The mutation-selection sce- rates, we are dealing with quasispecies → irregular dynamics nario was found to provide an appropriate molecular basis with increasing numbers of dominant subspecies → stochas- for understanding also virus evolution [For a recent survey tic hypercycle dynamics (Fig. 9). see the contributed volume (Domingo and Schuster 2016)]. More complex systems, for example, bacteria and popula- tions of cancer cells, were found to be describable by qua- Concluding remarks sispecies theory as well (Bertels et al. 2017; Covacci and Rappuoli 1998; Napoletani et al. 2013; Brumer et al. 2006). The model presented here has been conceived with modu- In the strong mutation scenario Darwin’s view of evo- lar structure in the sense that different mechanisms can be lution has to be modified. Not a single fittest genotype is applied for each of the three basic components. Here, it has selected but a uniquely defined distribution of genotypes, been presented in its simplest form. Each of the three mod- which is represented by the largest eigenvector of a value ules, competition, cooperation and variation, can be made matrix that represents the long time or stationary solution arbitrarily complex. Variation, for example, can be extended of Eigen’s mutation-selection equation. The mean fitness of to include more elaborate mutation mechanisms and recom- populations is not always optimized since situations can be bination as well as environmental influences. Even in case constructed in which the fitness is decreasing in the approach of viruses the reproduction mechanism is commonly much towards the stationary state. A trivial but illustrative example more elaborate and comprises a whole molecular machinery of decreasing fitness during evolution considers a homoge- instead of a single enzyme. Virus reproduction may include neous population consisting exclusively of fittest genotypes also the defense system of the host, epigenetic phenomena as initial condition: Mutations introduce mutants into the could be taken into account through simultaneous consid- population and since they have lower fitness by definition the eration of several generations, and for higher organisms the mean fitness is doomed to decrease. Such situations, how - real challenge in reproduction is to deal with the enormous ever, are rather rare and Darwinian optimization still remains complexity of development in a form that is simple enough a very powerful heuristic that applies to almost all scenarios. for modeling. Cooperation at the molecular level could also For error rate parameters exceeding a critical value p , t he cr involve reproductive autocatalytic networks whereas social largest eigenvector approaches the uniform distribution over phenomena in reproductive groups or societies represent the the entire sequence space, which is the exact solution for currently highest step in the open ended complexity increase the value p =  p leading to incorporations of correct and of biological evolution. Cooperation has been frequently incorrect nucleotides with equal probabilities—for binary modeled by game theory Maynard Smith (1982); Hofbauer sequences this happens at  p = 1 −  p = . In realistic popu- and Sigmund (1998). There is no limitation to make the lations, the uniform distribution is incompatible with a dis- model more complex, the problem evidently is to include crete quasispecies (17). Instead populations are observed the desired phenomena but to keep the model sufficiently that migrate randomly through sequence space Higgs and simple for mathematical analysis or simulation. Derrida (1991); Huynen et al. (1996). In the simple form in which the model was introduced In the second half of the twentieth century, most of the here, it has been tested experimentally by in vitro evolution molecular insights into reproduction and inheritance came experiments [For an overview of early works on this subject from viruses and bacteria and a high percentage of molecu- see (Spiegelman 1971; Biebricher 1983); as a recent review lar biologists thought that the basic regulation mechanisms we mention (Joyce 2007)]. The kinetic equations describing of gene activities are understood. Eukaryotic cells, however, replication and mutation were introduced 1971 by Manfred are not “giant bacteria”. Although the genetic code is the Eigen in his scholarly written paper on self-organization of same, the gene expression and inheritance system of higher biological macromolecules (Eigen 1971). Eigen’s mutation- organisms are different from the prokaryotic one and much selection equation describes the evolution of the distribution more complex. A true wealth of information on eukaryotic 1 3 422 European Biophysics Journal (2018) 47:403–425 cells has been discovered in the past fifty years, but gene autocatalytic systems consist almost always of complex expression in animals, plants, and fungi is still a subject of reaction networks rather than single-step reactions [as current cutting-edge research. Most of the recently revealed examples for attempts to construct simple systems of this gene expression regulation mechanisms are subsumed under kind see (McCaskill 1997; Wlotzka and McCaskill 1997)]. the notion of epigenetics for which an operational definition John Maynard Smith and Eörs Szathmáry collected a true has been proposed at the Cold Spring Harbor Meeting in the wealth of evidence for the historic occurrence of such year 2008 (Berger et al. 2009): major evolutionary transitions (Maynard Smith and Sza- An epigenetic trait is a stably heritable phenotype result- thmáry 1995) in the evolution of life. ing from changes in a chromosome without alterations in It is illustrative to think about transitions in terms of the DNA sequence. thresholds: Up to a certain value of the transition determin- The diversity of epigenetic effects on gene regulation is ing parameter the typical feature is hardly detectable and enormous. It ranges from specific methylation of DNA, in does not become evident before the parameter exceeds the particular cytosine methylation in position 5 of CpG ele- transitions value. Accordingly, such thresholds correspond ments (Zemach et al. 2010), histone methylation and acety- to sharp transitions. Nevertheless, it appears useful to be less lation (Lawrence et al. 2016) to post-transcriptional RNA- fussy and to accept the notion of threshold also for smooth methylation of adenine in position 6 (Barbosa Dogini et al. transitions. On the three faces of the coordinate system 2014; Yue et al. 2015) and small interfering RNAs (He and (Fig. 5) we observe an error threshold between selection and Hanon 2004). Epigenetics provides an extremely diverse, random replication, a cooperation threshold between selec- complex and flexible richness of regulatory actions on genes, tion and symbiosis, and a mutation threshold that separates which so far was not yet cast into a comprehensive theory the regime of independent replication of all subspecies from and precisely this is one of the greatest challenges for the mutual support through frequent mutation. future of evolutionary biology. Understanding evolution implies knowledge on the rela- There is neither a convincing theoretical model nor tion between genotypes being DNA or RNA sequences and experimental evidence that Darwinian evolution leads to phenotypes giving rise to all fitness relevant parameters. The an obligatory increase in complexity. The combination metaphor of an abstract fitness landscape has been originally of competitive selection and cooperation, however, may introduced by Sewall Wright for the purpose of illustration lead from one level of complexity to the next higher one (Wright 1922). In formal mathematical terms such a relation by integration of competitors through cooperation (Sza- can be modeled as a mapping from a genotype or sequence thmáry and Maynard Smith 1995; Schuster 1996). The space into fitness values. In molecular structural biology evolution model presented here proposes a mechanism such a mapping is split into two parts: (i) a mapping from for this integration of competitors and identifies the abun- sequences into structures or genotypes into phenotypes, and dance of resources as one driving force towards higher (ii) a second mapping assigning fitness values to structures complexity. This simple model distinguishes four steps or phenotypes (Schuster 2016). Computer models of RNA (Schuster 1996): (i) Initially the systems consists of inde- evolution with sequence–secondary structure–fitness map- pendent replicators competing for a single resource, (ii) pings have been studied in the past (Fontana and Schuster the capability of cooperative interaction allows to form 1987; Fontana et al. 1989; Fontana and Schuster 1998a, b an autocatalytic network, which couples the replicators and these studies provided the basis for a definition of evo- and suppresses competition but makes the network vulner- lutionary nearness of phenotypes in the presence of neu- able to exploitation by parasites, which consume resources trality (Stadler et al. 2001). With more and more data on without contributing a share to the common properties, sequences and fitness values of mutants becoming currently (iii) the members of the autocatalytic network are sepa- available fitness landscapes may also be determined directly rated from the environment by means if a suitable bound- by experiment (Kouyos et al. 2012) and it is not risky to ary that prevents the system from exploitation and allows predict that genotype–phenotype relations will become an for the formation of a new unit at a hierarchically higher integral part of evolution models in the near future. Then level, and (iv) the individual units at the higher level diver- evolution can be described in a self-contained manner where sify by variation, compete for common resources, Darwin- the genotype–phenotype mapping is a genuine part of the ian selection sets in and takes place now a the higher level. evolving system. The previously autonomous units at the lower level lost Acknowledgements Open access funding provided by University of their autonomy at least in part when they were integrated Vienna. This contribution presents the continuation of a twenty five into the higher unit of selection. Although modeling major years lasting, wonderful cooperation with Manfred Eigen, which has transitions as shown in "Competition and cooperation" is been initiated in 1968 when the author was Post-Doc in the Eigen labo- ratory at the Max Planck-Institute for Physical Chemistry in Göttingen, not difficult, suitable experimental molecular models are Germany. The author wants to thank for the unique opportunity of very hard to conceive, because second- and higher-order 1 3 European Biophysics Journal (2018) 47:403–425 423 participating in the development of Manfred Eigen’s visionary molecu- Dobzhansky T, Ayala FJ, Stebbins GL, Valentine JW (1977) Evolution. lar view of life. Many years of cooperation and continuous discussions W. H. Freeman & Co., San Francisco with Professors Karl Sigmund, Josef Hofbauer, Walter Fontana, Peter Domingo E, Schuster P (eds) (2016) Quasispecies: from theory to F. Stadler, Ivo L. Hofacker, Christoph Flamm and others are gratefully experimental systems, current topics in microbiology and immu- acknowledged. nology, vol 392. Springer, Berlin Domingo, E., Schuster, P.: What is a quasispecies? Historical origins and current scope. In: E. Domingo, P. Schuster (eds.) Qua- Open Access This article is distributed under the terms of the Crea- sispecies: From Theory to Experimental Systems, Current tive Commons Attribution 4.0 International License (http://creat iveco Topics in Microbiology and Immunology, vol. 392, chap. 1, mmons.or g/licenses/b y/4.0/), which permits unrestricted use, distribu- pp. 1–22. Springer-Verlag, Berlin (2016) tion, and reproduction in any medium, provided you give appropriate Doob JL (1942) Topics in the theory of Markoff chains. Trans Am credit to the original author(s) and the source, provide a link to the Math Soc 52:37–64 Creative Commons license, and indicate if changes were made. Doob JL (1945) Markoff chains—Denumerable case. Trans Am Math Soc 58:455–473 Dykhuizen DE, Hartl DL (1983) Selection in chemostats. Microbiol Rev 46:150–168 References Ebeling W, Mahnke R (1979) Kinetics of molecular replication and selection. Zagadnienia Biofizyki Współczesnej 4:119–128 Albery WJ, Knowles JR (1976) Evolution of enzyme function and the Eigen M (1971) Selforganization of matter and the evolution of bio- development of catalytic efficiency. Biochemistry 15:5631–5640 logical macromolecules. Naturwissenschaften 58:465–523 Albery WJ, Knowles JR (1977) Efficiency and evolution of enzyme Eigen M, McCaskill J, Schuster P (1989) The molecular quasispe- catalysis. Angew Chem Internat Ed Engl 16:285–293 cies. Adv Chem Phys 75:149–263 Baake E, Gabriel W (1999) Biological evolution through mutation, Eigen M, Schuster P (1977) The hypercycle. A principle of natural selection, and drift: an introductory review. In: Stauffer D (ed) self-organization. Part A: emergence of the hypercycle. Natur- Annual review of computational physics VII. World Scientific, wissenschaften 64:541–565 Singapore, pp 203–264 Eigen M, Schuster P (1978) The hypercycle. A principle of natural Berger SL, Kouzarides T, Shiekhattar R, Shilatifard A (2009) An oper- self-organization. Part B: the abstract hypercycle. Naturwis- ational definition of epigenetics. Genes Dev 23:781–783 senschaften 65:7–41 Bertels F, Gokhale CS, Traulsen A (2017) Discovering complete qua- Eigen M, Schuster P (1978) The hypercycle. A principle of natural sispecies in bacterial genomes. Genetics 206:2149–2157 self-organization. Part C: the realistic hypercycle. Naturwis- Biebricher CK (1983) Darwinian selection of self-replicating RNA senschaften 65:341–369 molecules. In: Hecht MK, Wallace B, Prance GT (eds) Evolu- Ewens WJ, Lessard S (2015) On the interpretation and relevance tionary biology, vol 16. Plenum Publishing Corporation, New of the fundamental theorem of natural selection. Theor Popul York, pp 1–52 Biol 104:59–67 Biebricher CK, Eigen M, Gardiner WC Jr (1983) Kinetics of RNA Feller W (1940) On the integro-differential equations of purely replication. Biochemistry 22:2544–2559 discontinuous Markoff processes. Trans Am Math Soc Biebricher CK, Eigen M, William C, Gardiner J (1984) Kinetics of 48:488–515 RNA replication: plus-minus asymmetry and double-strand for- Fisher RA (1930) The genetical theory of natural selection. Oxford mation. Biochemistry 23:3186–3194 University Press, Oxford Biebricher CK, Eigen M, William C, Gardiner J (1985) Kinetics of Fisher RA (1958) The genetical theory of natural selection, 2nd edn. RNA replication: competition and selection among self-repli- Dover Publications Inc, New York cating RNA species. Biochemistry 24:6550–6560 Fontana W, Schnabl W, Schuster P (1989) Physical aspects of evolu- Brumer Y, Michor F, Shakhnovich EI (2006) Genetic instability and tionary optimization and adaptation. Phys Rev A 40:3301–3321 the quasispecies model. J Theor Biol 241:216–222 Fontana W, Schuster P (1987) A computer model of evolutionary opti- Bryson V, Szybalski W (1952) Microbial selection. Science 116:45–46 mization. Biophys Chem 26:123–147 Covacci A, Rappuoli R (1998) Helicobacter pylori: molecular evolu- Fontana W, Schuster P (1998a) Continuity in evolution. On the nature tion of a bacterial quasi-species. Curr Opt Microbiol 1:96–102 of transitions. Science 280:1451–1455 Crow JF, Kimura M (1970) An introduction to population genetics Fontana W, Schuster P (1998b) Shaping space. The possible and the theory. Sinauer Associates, Sunderland attainable in RNA genotype-phenotype mapping. J Theor Biol Dogini BD, Pascoal VD, Avansini SH, Vieira AS, Campos TC, Lopes- 194:491–515 Cendes I (2014) The new world of RNAs. Genet Mol Biol 37(1 Gillespie DT (1976) A general method for numerically simulating the suppl):285–293 stochastic time evolution of coupled chemical reactions. J Comp de Pasquale F, Tartaglia P, Tombesi P (1980) Stochastic approach to Phys 22:403–434 chemical instabilities. Anomalous fluctuations transient behavior. Gillespie DT (1977) Exact stochastic simulation of coupled chemical Lettere al Nuovo Cimento 28:141–145 reactions. J Phys Chem 81:2340–2361 de Pasquale F, Tartaglia P, Tombesi P (1982) New expansion technique Gillespie DT (2007) Stochastic simulation of chemical kinetics. Annu for the decay of an unstable state. Phys Rev A 25:466–471 Rev Phys Chem 58:35–55 de Pasquale F, Tartaglia P, Tombesi, P (1982) Transient-behavior near Gillespie JH (1983) Some properties of finite populations experienc- an instability point. Vectorial stochastic representation of the ing strong selection and weak mutation. Am Nat 121:691–708 Malthus. Nuovo Cimento B 69:228–238 Goel NS, Richter-Dyn N (1974) Stochastic models in biology. Aca- Dean AM, Thronton JW (2007) Mechanistgic approaches to the study demic Press, New York of evolution: the functional synthesis. Nat Rev Genet 8:675–688 Gorini L, Beckwith JR (1966) Suppression. Annu Rev Microbiol Deuflhard P, Huisinga W, Jahnke T, Wulkow M (2008) Adaptive dis- 20:401–422 crete Galerkin methods applied to the chemical master equation. Hamming RW (1950) Error detecting and error correcting codes. Bell SIAM J Sci Comput 30:2990–3011 Syst Tech J 29:147–160 1 3 424 European Biophysics Journal (2018) 47:403–425 Hamming RW (1986) Coding and information theory, 2nd edn. Pren- Nicolis G, Prigogine I (1977) Self-organization in nonequilibrium sys- tice-Hall, Englewood Cliffs tems. Wiley, New York Hartman PE, Roth JR (1973) Mechanisms of suppression. Adv Genet Novick A, Szillard L (1950) Description of the chemostat. Science 17:1–105 112:715–716 He L, Hanon GJ (2004) Micro RNAs: small RNAs with a big role in Novick A, Szillard L (1950) Experiments with the chemostat on gene regulation. Nat Rev Genet 5:522–531 spontaneous mutations of bacteria. Proc Natl Acad Sci USA Higgs PG, Derrida B (1991) Stochastic models for species formation in 36:708–719 evolving populations. J Phys A Math Gen 24:L985–L991 Noyes RM, Field RJ, Körös E (1972) Oscillations in chemical sys- Hofbauer J, Mallet-Paret J, Smith HL (1991) Stable periodic solutions tems. I. Detailed mechanism in a system showing temporal for the hypercycle system. J Dyn Diff Equ 3:423–436 oscillations. J Am Chem Soc 94:1394–1395 Hofbauer J, Sigmund K (1998) Evolutionary games and population Okasha S (2008) Fisher’s fundamental theorem of natural selec- dynamics. Cambridge University Press, Cambridge tion—a philosophical analysis. Br J Philos Sci 59:319–351 Huang CI, Tu MF, Lin HH, Chen CC (2017) Variation approach Park JM, Muñoz E, Deem MW (2010) Quasispecies theory for finite to error threshold in generic fitness landscape. Chin J Phys populations. Phys Rev E 81:e011,902 55:606–618 Phillipson PE, Schuster P (2009) Modeling by Nonlinear Differen- Husimi Y (1989) Selection and evolution of a bacteriophages in cell- tial Equations. Dissipative and Conservative Processes, World stat. Adv Biophys 25:1–43 Scientific Series on Nonlinear Science A, 69th edn. World Sci- Husimi Y, Nishigaki K, Kinoshita Y, Tanaka T (1982) Cellstat—a con- entific, Singapore tinuous culture system of a bacteriophage for the study of the Phillipson PE, Schuster P, Kemler F (1984) Dynamical machinery of mutation rate and the selection process at the DNA level. Rev a biochemical clock. Bull Math Biol 46:339–355 Sci Instr 53:517–522 Plutynski A (2006) What was Fisher’ fundamental theorem of natural Huynen MA, Stadler PF, Fontana W (1996) Smoothness within rug- selection and what was it for? Stud Hist Phil Biol Biomed Sci gedness. The role of neutrality in adaptation. Proc Natl Acad Sci 37:50–82 USA 93:397–401 Prelich G (1999) Suppression mechanisms. themes from variations. Jahnke T, Huisinga W (2007) Solving the chemical master equation Trends Genet 15:261–266 for monomolecular reaction systems analytically. J Math Biol Price GR (1972) Fisher’s ‘fundamental theorem’ made clear. Ann 54:1–26 Hum Genet 36:129–140 Jones BL, Enns RH, Rangnekar SS (1976) On the theory of selection Sagués F, Epstein IR (2003) Nonlinear chemical dynamics. J Chem of coupled macromolecular systems. Bull Math Biol 38:15–28 Soc Dalton Trans 2003:1201–1217 Jones BL, Leung HK (1981) Stochastic analysis of a non-linear model Schmidt LD (2005) The engineering of chemical reactions, 2nd edn. for selection of biological macromolecules. Bull Math Biol Oxford University Press, New York 43:665–680 Schnabl W, Stadler PF, Forst C, Schuster P (1991) Full characteri- Joyce GF (2007) Forty years of in vitro evolution. Angew Chem Inter- zatin of a strang attractor. Chaotic dynamics on low-dimen- nat Ed 46:6420–6436 sional replicator systems. Phys D 48:65–90 Joyce P, Rokyta DR, Beisel CL, Orr HA (2008) A general extreme Schuster P (1996) How does complexity arise in evolution. Nature’s value theory model for the adaptation of DNA sequences under recipe for mastering scarcity, abundance, and unpredictability. strong selectiion and weak mutation. Genetics 180:1627–1643 Complexity 2/1(1):22–30 Kimura M (1955) Solution of a process of random genetic drift with a Schuster P (2016) Increase in complexity and information through continuous model. Proc Natl Acad Sci USA 41:144–150 molecular evolution. Entropy 18:e397 Kimura M (1955) Stochastic processes and distribution of gene fre- Schuster P (2016) Major transitions in evolution and in technology. quencies under natural selection. Cold Spring Harb Symp Quant What they have in common and where they differ. Complexity Biol 20:33–53 21/4(4):7–13 Kimura M (1983) The neutral theory of molecular evolution. Cam- Schuster P (2016) Quasispecies on fitness landscapes. In: bridge University Press, Cambridge E. Domingo, P. Schuster (eds) Quasispecies: from theory to Kolmogorov AN (1931) Über die analytischen Methoden in der Wahr- experimental systems, Current Topics in Microbiology and scheinlichkeitsrechnung. Math Ann 104:415–458 In German Immunology, vol. 392, chap. 4. Springer-Verlag, Berlin, pp Koltermann A, Kettling U (1997) Principles and methods of evolution- 61–120. ary biotechnology. Biophys Chem 66:159–177 Schuster P (2016) Some mechanistic requirements for major transi- Kouyos RD, Leventhal GE, Hinkley T, Haddad M, Whitcomb JM, tions. Phil Trans Roy Soc Lond B 371:e20150,439 Petropoulos CJ, Bonhoeffer S (2012) Exploring the complexity Schuster P (2016) Stochasticity in processes. Fundamentals and appli- of the HIV-1 fitness landscape. PLos Genet 8:e1002,551 cations in chemistry and biology. Springer series in synergetics. Lawrence M, Daujat S, Schneider R (2016) Lateral thinkiing: how Springer, Berlin histone modifications regulate gene expression. Trends Genet Schuster P (2017) A mathematical model of evolution. MATCH Com- 32:42–56 munications in Mathematical and in Computer Chemistry 78, Maynard Smith J (1982) Evolution and the theory of games. Cambridge submitted University Press, Cambridge Schuster P (2017) Molecular evolution between chemistry and biol- Maynard Smith J, Szathmáry E (1995) The major transitions in evolu- ogy. The interplay of competition, cooperation, and mutation. tion. W. H. Freeman-Spektrum, Oxford European Biophysics Journal 46, submitted McCaskill JS (1997) Spatially resolved in vitro molecular ecology. Schuster P, Sigmund K (1985) Dynamics of evolutionary optimization. Biophys Chem 66:145–158 Ber Bunsenges Phys Chem 89:668–682 Mills DR, Peterson RL, Spiegelman S (1967) An extracellular Darwin- Schuster P, Sigmund K, Wolff R (1978) Dynamical systems under ian experiment with a self-duplicating nucleic acid molecule. constant organization I. Topological analysis of a family of non- Proc Natl Acad Sci USA 58:217–224 linear differential equations - A model for catalytic hypercycles. Napoletani D, Signore M (2013) Cancer quasispecies and stem-like Bull Math Biol 40:734–769 adaptive aneuploidy. F1000Research 2:e268 1 3 European Biophysics Journal (2018) 47:403–425 425 Schuster P, Sigmund K, Wolff R (1979) Dynamical systems under con- Tarazona P (1992) Error thresholds for molecular quasispecies as phase stant organization III. Cooperative and competitive behavior of transitions: From simple landscapes to spin glasses. Phys Rev A hypercycles. J Diff Equ 32:357–368 45:6038–6050 Schuster P, Sigmund K, Wolff R (1980) Dynamical systems under con- Thompson CJ, McBride JL (1974) On Eigen’s theory of the self-organ- stant organization II. Homogenoeus growth functions of degree ization of matter and the evolution of biological macromolecules. p = 2 . SIAM J Appl Math 38:282–304 Math Biosci 21:127–142 Sniegowski PD, Gerrish PJ (2010) Beneficial mutations and the dynam- Wiehe T (1997) Model dependency of error thresholds: The role of ics of adaptation in asexual populations. Phil Trans R Soc B fitness functions and contrasts between the finite and infinite sites 356:1255–1263 models. Genet Res Camb 69:127–136 Spiegelman S (1971) An approach to the experimental analysis of pre- Wlotzka B, McCaskill JS (1997) A molecular predator and its prey: cellular evolution. Quart Rev Biophys 4:213–253 coupled isothermal amplification of nucleic acids. Chem Biol Stadler BRM, Stadler PF, Wagner GP, Fontana W (2001) The topology 4:25–33 of the possible: Formal spaces underlying patterns of evolution- Wright S (1922) Coefficients of inbreeding and relationship. Am Nat ary change. J Theor Biol 213:241–274 56:330–338 Strunk G, Ederhof T (1997) Machines for automated evolution experi- Yue Y, Liu J, He C (2015) RNA N -methyladenosine methylation ments in  vitro based on the serial-transfer concept. Biophys in post-ttranscriptional gene expression regulation. Genes Dev Chem 66:193–202 29:1343–1355 Swetina J, Schuster P (1982) Self-replication with errors - A model for Zemach A, McDaniel IE, Silva P, Zilberman D (2010) Genome-wide polynucleotide replication. Biophys Chem 16:329–345 evolutionary analysis of eukaryotic DNA methylation. Science Szathmáry E, Maynard Smith J (1995) The major evolutionary transi- 328:916–919 tions. Nature 374:227–232 1 3 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png European Biophysics Journal Springer Journals

Molecular evolution between chemistry and biology

Free
23 pages
Loading next page...
 
/lp/springer_journal/molecular-evolution-between-chemistry-and-biology-K4TAqh8mhL
Publisher
Springer International Publishing
Copyright
Copyright © 2018 by The Author(s)
Subject
Life Sciences; Biochemistry, general; Biological and Medical Physics, Biophysics; Cell Biology; Neurobiology; Membrane Biology; Nanotechnology
ISSN
0175-7571
eISSN
1432-1017
D.O.I.
10.1007/s00249-018-1281-7
Publisher site
See Article on Publisher Site

Abstract

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, ( ) competition and mutation, (  ) cooperation and competition, and (  ) cooperation and mutation, are analyzed. Changes in population dynam- ics that are induced by bifurcations and threshold phenomena are discussed. Keywords Competition · Cooperation · Darwinian optimization · Error threshold · Mutation · Neutral evolution · Selection Introduction The model focusses on three basic processes: (i) fitness- driven competition through differential reproductive success, Biology is centered around evolutionary thinking or under- (ii) reproduction-relevant cooperation between competitors, standing biology implies understanding evolution as Theo- and (iii) reproduction-induced variation. The model is con- dosius Dobzhansky pointed out clearly in his famous book ceived with a modular structure and allows for the imple- on evolution: “Nothing in biology makes sense except in the mentation of different, more or less complicated mechanisms light of evolution” (Dobzhansky et al. 1977). Pure Darwin- for the processes (i), (ii), and (iii). For example, variation ian evolution is a simple process but its embedding in nature may be implemented by mutation, by recombination or by renders it complex: Natural selection would follow uncom- both. Here, we shall apply the simplest conceivable mecha- plicated laws in a simple environment. In the light of current nisms: reproduction as single enzyme mediated replication molecular biology, there is need for a simple but comprehen- (Biebricher 1983), cooperation of competitors as hypercy- sive mathematical model of evolution to be able to account cle dynamics (Eigen and Schuster 1978a, b), and variation for modern genetics. The various epigenetic mechanisms as single point mutations based on the uniform error rate have to be part of any comprehensive model of evolution assumption (Swetina 1982). and to keep such a model amenable to analysis and handling, Usage of the notion evolution is often ambiguous and a molecular details must be reduced to a coarse-grained level. precise definition is desirable. Here evolution is understood This article deals with a flexible model of evolution under as a process based on reproduction of a genotype being a defined environmental conditions. We present a concise and DNA or an RNA sequence that carries encoded information comprehensive review of work that was published elsewhere on the formation of a phenotype, which is evaluated with (Schuster 2016a, d, 2017a) together with a few new results. respect to success in reproduction. The evolutionary pro- cess is built upon two foundations: (i) the dynamics at the population level and (ii) the environment dependent encod- Special Issue: Chemical Kinetics, Biological Mechanisms and ing of phenotypes in genotype. At the molecular level, the Molecular Evolution. latter boils down to sequence–structure–function relations * Peter Schuster (Schuster 2016). The simplest systems that are capable of pks@tbi.univie.ac.at evolution in the sense of the given definition are nothing but special autocatalytic reactions involving polynucleotides Institut für Theoretische Chemie, Universität Wien, under suitable conditions. Währingerstraße 17, 1090 Wien, Austria Vol.:(0123456789) 1 3 404 European Biophysics Journal (2018) 47:403–425 In the next "Preliminaries", some prerequisites are pre- be implemented as long as they can be cast in kinetic equa- sented for the model, which is introduced in "A minimal tions. In a separate "Basic processes", the processes that will mathematical model for the evolution of molecules". The be used to model evolution are introduced. The following model comes in a deterministic version based on kinetic dif- section describes the deterministic and stochastic methods, ferential equations or it is formulated as a stochastic process which are applied to find solutions of the kinetic equations. modeled by means of chemical master equations (Schuster Finally, we review some fundamental features of autocataly- 2016). Although almost no analytic solutions are available sis since this is the chemical counterpart of reproduction. for master equations derived from nonlinear reaction kinet- ics in two or more variables, the stochastic version of the Idealized environments model can be studied by efficient simulation methods for the calculations of trajectories (Gillespie 1977, 2007). In Environments that allow for investigations of observa- three sections, solution curves for the three two-dimensional tions as functions of one or few parameters with everything subspaces, (  ) reproduction and mutation, (  ) reproduc- else being constant require elaborate design in the form of tion & cooperation, and ( ) cooperation & mutation, are pre- sophisticated experiments since devices controlling envi- sented, analyzed and interpreted. Then follows a brief dis- ronmental conditions may be quite involved. In theoretical cussion of results obtained with the full three-dimensional approaches, often the silent assumption is made that there model, reproduction & cooperation & mutation and in the exists a hypothetical machinery, which takes care of fixing final section we discuss the simple model in the context of parameter values as needed for the mathematical analysis. the complex processes in nature. 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 evolu- Preliminaries tion, which result from reproduction, symbiotic cooperation and variation, and therefore impacts on evolution caused by Darwinian evolution is often—and incorrectly—seen as a changes in the environment are intended to be kept to a mini- synonym for optimization mainly because Ronald Fisher’s mum. To reduce the influence of the environment as much fundamental theorem of natural selection (Fisher 1930, as possible we shall assume a control device in the form of 1958; Price 1972; Ewens and Lessard 2015) focusses on a simple flow reactor (Fig.  1; see also Schmidt 2005). More non-decreasing mean fitness (t) in evolving populations. In elaborate reactors, which keep, for example, the numbers of the simplest form derived from idealized models, the time bacterial cells constant have been designed and implemented derivative of the mean fitness is the variance of the fitness (Novick and Szillard 1950a, b; Bryson 1952). The flow reac- values and therefore a non-negative quantity. In reality, mean tors called chemostat, cellstat or turbidostat and other exper- fitness is the result of many factors and as Fisher was been imental devices for monitoring and controlling evolution in certainly aware, only a few cases—like single locus genet- the laboratory may serve as examples (Husimi et al. 1982; ics—fulfil the theorem in pure form (for elaborate discus- Dykhuizen and Hartl 1983; Husimi 1989; Koltermann and sions see the more recent literature Plutynski 2006; Okasha Kettling 1997; Strunk and Ederhof 1997). 2008). Implementation of a physical device rather than appli- Evolution is intimately related to the environment in cation of idealized assumptions like constant population which it takes place and accordingly environment and envi- size is required for the stochastic description of evolution. ronmental changes are major factors shaping evolutionary An illustrative example where the deterministic approach processes. Here, we are primarily interested in the internal yields a stable solution whereas the corresponding stochastic dynamics and hence a well-defined and controllable environ- system is unstable is provided by the linear birth-and-death ment is required. For this goal, we introduce a flow reactor in process (Goel and Richter-Dyn 1974), which is described by f d "Idealized environments", which represents a simple device the kinetic equations  + ⟶ 2 and ⟶∅ . For equal that is not only useful for performing evolution experi- birth and death rate parameters, f = d , the population num- ments but also provides at the same time a suitable setup ber of the deterministic system stays constant whereas in for theoretical modeling. More complex environments can 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 The term linear reaction kinetics is used in this contribution for words, the system is unstable despite a (marginally) stable reactions and flow terms of zero and first order, which give rise to deterministic solution. The direct incorporation of constant constants or linear functions in the kinetic ODEs. First-order kinetics population size into Fisher’s selection (1930) or Eigen’s gives rise to exponential functions, which become linear in ln[]∕t -plots. equation (1971) leads to an instability of the same kind since 1 3 European Biophysics Journal (2018) 47:403–425 405 in the replication of every molecular species  —the number of catalytic terms is n 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, 1978, 1978): The catalyzed reactions 2 +  ⟵ +  +  (i = 1,… , n; i mod n) i i+1 i i+1 ⋯ ←  ←  ←  ← ⋯ ←  ←  ← … n 1 2 n−1 n 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 sim- plest 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 replica- tion kinetics of the multistep process, Q ⋅k [] ji i (1) +  + →→→ +  +  , i j i Fig. 1 The continuous-flow stirred-tank reactor (CSTR). The fig- ure sketches a device for controlling the environmental condition with E being a replicase enzyme, the reaction rate is obtained of evolution experiments. The material needed for reproduction is as the product of two parameters: Q ⋅ f where f = k [] is ji i i i subsumed by A, it flows into the reactor with a (volumetric) flow the fitness of the template that depends on the availability rate r   [V  /  t]≡[volume  /  time] in form of a solution with concentra- tion []= a or number density []= A . In the reactor molecules of resources here the concentration [ A]. The dimensionless 0 0 (i = 1, … , n) are reproduced and A is consumed. The volume V factors Q with i, j = 1,… , n are understood as the elements ji of the reactor is constant and hence reaction mixture compensating of a mutation matrix Q ={Q } and provide the probability ij the volume increase through influx of stock solution has to flow out to obtain  as a copy of the template  that can be either of the reactor. The mean residence time of a volume element in the j i CSTR is  = V∕r [t]. Inflow and outflow of materials are handled as correct (  ≡  ) or error-prone (  ,  j ≠ i ). Conservation of a ⋅r j i j virtual chemical reactions or pseudoreactions: ∗⟶ for inflow and, r r probabilities requires: Q = 1 since each copy has to be ji j=1 ⟶∅ and  ⟶∅ (i = 1, … , n) for outflow either correct, Q , or incorrect, Q = 1 − Q . For ii ji ii j=1,j≠i the sake of simplicity, binary sequences will be considered here and we distinguish only between correct and incorrect fluctuations are not self-regulating (Jones and Leung 1981). base pairs. Implementation of the system in the flow reactor provides A useful simplifying approximation is made by the uniform stability due to balance control of inflow and outflow mod- a ⋅r 0 error rate model (Swetina 1982): The error per nucleotide and eled by the pseudoreactions and ⟶∅ as well as ∗⟶ replication event, p, is assumed to be independent of the posi- ⟶∅ (i = 1, … , n) (Fig. 1). tion of the nucleotide along the sequence and the nature of the nucleotide to be complemented. Then all elements of the Basic processes mutation matrix can be expressed by a simple formula, The minimal system for modeling evolution of molecules d l−d l d ij ij ij Q = p (1 − p) = p  with  = , ji (2) is based on the three classes of processes: (i) competitive 1 − p reproduction, (ii) symbiotic cooperation, and (iii) repro- with only three parameters: (i) the sequence length of duction based variation. For the minimal model we shall the RNA molecules l, (ii) the mutation rate p, and (iii) choose the simplest possible chemical reactions. Reproduc- the Hamming distance (Hamming 1950, 1986) between tion will be modeled by simple replication,  +  ⟶ 2 i i the two sequences interrelated by the mutation process, with f being the fitness of species  and competition being i i d = d (X , X ) . Without changing important results for the ij H i j introduced through living on the same resource A. Symbi- purposes pursued here, the analysis of the model is sub- ontic cooperation is introduced as catalyzed replication, ij stantially simplified by the assumption of constant chain +  +  ⟶ 2 +   (i, j = 1, … , n) , where  repre- i j i j i lengths l, which is also consistent with the restriction to sents the template and  is the catalyst. In its most general form—every molecule  has the potential to act as catalyst 1 3 406 European Biophysics Journal (2018) 47:403–425 point mutations since point mutations do not change chain Equilibrium fluctuations in conventional chemical reac- lengths by definition. tion kinetics follow an approximate N -law and hence play As an alternative to Eq. (1), mutation can be seen as almost no role in systems with particle numbers that are a consequence of DNA change—damage and imperfect typical for chemical systems. In biology a different scenario repair—during the whole life time of an organism, which is encountered. For example, every new variant originat- is the idea in the Crow–Kimura mutation model (2009), ing from mutation has to start out from a single copy. The pp 264–266]: deterministic approach commonly uses continuous vari- ables, which can only be an approximation to reality, since k [] ji numbers of molecules or biological entities are integers by (3) +  + →→→ 2 +  and  →→→ . i i i j definition. Continuous concentrations can adopt arbitrarily Interestingly, both models (1) and (3) although different with small values whereas stochastic variables cannot pass low respect to the underlying physics give rise to identical math- values beyond unity, because then molecular species go ematical problems (see e.g. Baake and Gabriel 1999 and extinct, and deterministic solutions become unrealistic and Schuster 2016, pp 76–78). differ strongly from the stochastic results. Large population size alone is not sufficient, each variable has to be suffi- Deterministic and stochastic approaches ciently large at every instant to guarantee similarity between stochastic and deterministic solutions. Two more phenom- Reaction mechanisms are commonly analyzed by determin- ena are relevant in the stochastic approach at low particle istic and stochastic approaches. The former translate the numbers and may lead to large standard deviations in the chemical reaction equations into kinetic differential equa- variables: (i) nucleation steps in reactions involving two or tions, which can be directly solved by mathematics, studied more molecules and (ii) multiple (quasi)stationary states of by means of qualitative analysis or investigated by numeri- the stochastic system. The stochastic system may approach cal integration. The theorems of existence and uniqueness any stationary state and the distribution of the possible out- of solutions of differential equations are applicable and a comes is determined by probabilities. Then the calculations single integration provides the complete information for a of expectation values and higher distribution moments build given input set. Competitive selection with nonzero mutation averages over different final states and may give rise to enor - leads to one unique asymptotically stable stationary state mous standard deviations. (Thompson 1974; Jones et al. 1976), whereas the long-time For the purpose of illustration, we consider the equilibra- dynamics of cooperative systems is much richer and multiple tion of the flow reactor as an example of a stochastic system stationary states, oscillations or deterministic chaos may be that can be fully analyzed by analytic calculations (Schuster observed (Schuster 1978; Schnabl et al. 1991). 2016, pp 436–441). The two pseudoreactions, Stochastic analysis in general in based on searching for a a r (4a) ∗ ⟶  and stochastic process that fits the model to be studied as closely as possible. Chemical reaction kinetics prefers master equa- r (4b) ⟶ ∅ , tions although the repertoire of analytical solutions is very are converted into the easy to solve kinetic differential limited. It is not difficult to write down a multivariate master equation equation but the derivation of analytical solutions is suc- cessful only in exceptional cases, for example for networks da −r t = a ⋅ r − a ⋅ r =(a − a) r with a(t)= n + n − a e 0 0 0 0 0 of monomolecular reactions (Jahnke et al. 2007; Deuflhard dt (4c) et al. 2008. If no analytical solutions are available informa- where n = a(0) . The master equation for the probability tion on the stochastic system can be obtained by trajectory sampling. The theoretical background for trajectory har- P (t)= P A(t)= n with number density A(t) being the dis- crete pendant of concentration a(t), vesting has been laid down by Andrey Kolmogorov (1931), Willy Feller (1940), and Joe Doob (1942, 1945). With elec- tronic computers now being generally available elaborate = r a P +(n + 1)P −(a + n)P , n ∈ ℕ , 0 n−1 n+1 0 n simulations of stochastic processes became possible. The (4d) more recent conception, analysis, and implementation of a simple but highly efficient algorithm by Daniel Gillespie (1976, 1977, 2007) provides a very useful tool for investiga- A stochastic quasistationary state is a state towards which the sys- tions of stochastic effects in chemical kinetics. Distributions tem converges stochastically in the long-time limit and around which it fluctuates. It is not an absorbing boundary, and if true asymptoti- of trajectories are characterized by expectation values and cally stable stationary states exist the system converges to one of higher moments, commonly only by variances or standard them in the limit t → ∞ although the mean time of convergence may deviations. be extremely long. 1 3 European Biophysics Journal (2018) 47:403–425 407 has the analytical solution The forward reaction of (5) leads to molecular self-enhance- ment: The resource A is consumed in the synthesis of the min{n ,n} −krt −rt n +n−2k replicator X and the process is catalyzed by the presence n e (1 − e ) −rt n−k a (1−e ) P (t)= a e of further molecules X. Depending on the molecularity k (n − k)! k=0 of the autocatalytic process as expressed by the value of (4e) m autocatalysis comes in different forms that, in essence, with n, n , a ∈ ℕ for the sharp initial condition P (0)=  . 0 0 n n,n fall into two classes: (i) simple or first-order autocatalysis In the limit t → ∞ the distribution (4e) converges to a Pois- with m = 1 and (ii) second- and higher-order autocatalysis son distribution with m ≥ 2 [(Phillipson and Schuster 2009, pp 18–27)]. For consistency, we add here also the uncatalyzed reaction and lim P (t)= P = exp(−a ) . (4f) n n 0 t→∞ n! call it zeroth order autocatalytic ( m = 0). All processes in closed systems converge to an equilibrium state and show First and second moment yield expectation values and stand- mass conservation that implies a conservation relation ard deviations A(0)+ X(0)= A(t)+ X(t)= N . The reversible autocatalytic � � −rt reactions converge to the equilibrium state: E A(t) = a +(n − a ) e and 0 0 0 � � (4g) −rt −rt A(t) = (a + n e )(1 − e ) . S ∶ a = N∕(1 + K), x = N ⋅ K∕(1 + K). 0 0 Although the expressions for the equilibria are the same The standard deviation starts out from (t = 0)= 0 and for the all reactions independently of the stoichiomet- approaches the equilibrium value  = a either from ric coefficient m, K = g ∕g , there are subtle differ- → ← below for a > n or from above for a < n after having 0 0 0 0 ences in the probability distributions, which become passed a maximum at important at small concentrations: The particle numbers are discrete quantities, change by integers only, and the 2n t = ln . (4h) max stoichiometric factors are X(t) X(t)− 1 ≡ X(t) or r n − a 2 0 0 X(t) X(t)− 1 X(t)− 2 ≡ X(t) rather than X(t) and In general the maximum if it exists is very flat (see Fig.  3) X(t) , respectively . The states with X(t)= 1 for m = 1 , or the as can be easily checked from the analytical expressions states X(t)= 1 and X(t)= 2 for m = 2 cannot react because through inspection: 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 � � (n − a ) � � n − a 0 0 0 0 states with X(t)= 0 and X(t)= 1 in the first- or second-order E A(t ) = a + and  A(t ) = . max 0 max 2n 2 n autocatalytic reactions require different normalizations for the stationary probabilities of the three systems: The flow reactor is one of the few rather exceptional cases (0) where the stochastic approach can be done completely by P = ; n ∈[0, N] , (6a) (1 + K) mathematical analysis. (1) Autocatalysis in the batch reactor P = ; n ∈[0, N − 1] , (6b) N N (1 + K) − K A batch reactor is an elaborate device that allows for per- (2) forming chemical reactions under controlled conditions N K P = ; n ∈[1, N − 2] , N N N−1 (1 + K) − K − N K − 1 without inflow and outflow (Schmidt 2005). Here, the term batch reactor is used to indicate that reactions are carried out (6c) in a well-mixed closed system under constant temperature where as before the stochastic variable is n = A . Mass con- and pressure, and in the long run approach a thermodynamic servation provides X = N − n . As expected the truncation (0) equilibrium state. of P is important for small values of N only. For first- In conventional chemistry autocatalysis is a rather rare order autocatalysis with N = 10 , k = k = 1.0 we obtain → ← phenomenon but in biology it represents the most impor- tant process since multiplication is just a special form of autocatalysis that in simplest form can be expressed by the We shall use different letters for the rate parameters: For m = 0 we reversible reaction use g ≡ h , for m = 1 g ≡ k , and for m = 2 g ≡ l . It is worth recalling −1 −1 −1 −2 −1 the different dimensions: h [t ], k [M t ], and l [M t ] that are the same in both directions. −−→ + m ←−− (m + 1) . (5) g Here the expression (x) = x!∕(x − n)! denotes the falling Pochham- → n mer symbol. 1 3 408 European Biophysics Journal (2018) 47:403–425 small t the factor A(t) is large and the factor(s) contain- ing 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. Self- enhancement 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 uncata- lyzed 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{(t)} is a measure of the random scatter in the max delay in the onset of the autocatalytic reaction (Table 1). (ad iii) Multiple final states give rise to an additional sto- chastic 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 ini- tial 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, respec- tively. A total population size of N = 1000 was chosen so Fig. 2 Comparison of fluctuations in the reversible zero- and first- that the one-standard-deviation fluctuation band of order order autocatalytic reaction  + m ⇌ (m + 1) with m = 0, 1 bf in N appears small and the deterministic solutions coincide the batch reactor. The two plots show the expectation values E A(t) (black) and E X(t) (red) together with the one standard deviation with the expectation values in the reference process A ⇌ X band E(t)± (t) (gray and pink) obtained by sampling of 10,000 tra- ( m = 0 ). The transient of the autocatalytic process ( m = 1 ) jectories that were calculated by Gillespie’s simulation method for the is different: It shows substantial broadening of the fluctua- uncatalyzed ( m = 0 ; top plot) and the first-order autocatalytic reac- tion band (Table 1) because of delayed onset of the reaction tion ( m = 1 ; bottom plot). Both reactions approach almost identical thermodynamic equilibrium states. Choice of parameters and initial before it narrows down to the equilibrium value. In the inter- −1 conditions: N = 1000 ; h = h = 1.5[t ], A(0)=1000, X(0)=0 (top mediate range the expectation values differ remarkably from → ← −1 −1 plot); k = k = 0.01[M   t ], A(0)=999, X(0)=1 (bottom). Solu- → ← the deterministic solution curves. As expected an increase tions of the corresponding kinetic differential equations are shown as in X leads to a decrease in the width of the fluctuation broken lines. 0 band. Second order autocatalysis ( m = 2 ) differs from first- order autocatalysis mainly by the size of the characteristic effects: Both autocatalytic reactions show broadening of the E(A)= 4.9951 , for second-order autocatalysis E(A)= 4.9605 fluctuation bands but the band width in the second-order compared to E(A)= 5 for the uncatalyzed process. case is about four times as wide. In particular, the scatter In principle, there are three major sources of random- in the waiting times until the first reaction events occurs is ness in autocatalytic reactions: (i) thermal fluctuations, (ii) much larger because we are dealing with two small factors, delayed onset of autocatalytic reactions and (iii) multiple (X − 1) ⋅ (X − 2) , in the expression for v(0). stationary states. (ad i) In the transients towards equilib- The standard deviation in the course of the reactions, (t) , rium the thermal fluctuation bands are essentially the same is shown in Fig. 3. Because of sharp initial conditions the in autocatalytic and conventional reactions as can be seen fluctuation band starts out from zero— (0)= 0 , increases, best from the comparison of standard deviations at equi- becomes broad in the intermediate range and settles down at librium (Fig.  2). The differences between conventional the equilibrium value (6). Substantial deviations between the and autocatalytic equilibrium densities can be recognized deterministic solution and the stochastic expectation value, numerically for very small particle numbers only (). (ad a(t) and the E A(t) or x(t) and E X(t) , respectively, are ii) The reaction rate for an autocatalytic reaction of order observed in the intermediate range. Accordingly, the stand- m is v(0)= kA(0) X(0) − l X(0) (with m ≥ 1 ). At m m+1 ard deviation goes through a pronounced maximum that is 1 3 European Biophysics Journal (2018) 47:403–425 409 those shown in Fig.  2 shows similarities except a twice as wide range for the values of the stochastic variables, A(0)+ X(0) ≥ A(t) ≥ 0 and 1 ≤ X(t) ≤ A(0)+ X(0) and expectation values and standard deviations approach zero at sufficiently long times. Again, we observe a sig- moid shape of the solution curves, for small initial val- ues ( X(0) < 5 ) the standard deviation (t) becomes large in the intermediate range (Fig. 3), and the deterministic curve deviates substantially from the stochastic expecta- tion values. Autocatalysis in the flow reactor Fig. 3 Comparison of the standard deviation in the revers- Implementation of autocatalytic reactions in the flow reac - ible first-order autocatalytic reaction  +  ⇌ 2 and the revers- tor provides additional insights into the different forms of ible uncatalyzed reaction  ⇌  . The reactions are recorded for the closed system fulfilling A(t)+ X(t)= A + X = const = N where randomness. In particular we are interested in multiple 0 0 A(0)= A and X(0)= X for both reactions. The figure presents 0 0 stationary states as a source of stochasticity (item iii). The the results of statistical evaluation of 10  000 trajectories obtained reaction equations for first order autocatalysis are: by computer simulations with Gillespie’s method (Gillespie 2007) (1) (1) a r for the autocatalytic reaction  (t)=  (t) (black) and for the (0) (7a) (0) ∗ ⟶ , uncatalyzed reaction  (t)=  (t) (red). Choice of parameters: −1 −1 −1 k = k = 0.01[M   t ] and h = h = 1.5[t ]; equilibr ium → ← → ← parameters K = k ∕k = h ∕h = 1 ; initial conditions: N = 1000 , → ← → ← (7b) +  ⟶ 2, X = 1 , A = 999 and N = 1000 , X = 0 , A = 1000 yielding the 0 0 0 0 numerical equilibrium values A = X = 500 . The equilibrium value of the standard deviation is practically the same for both reactions: (7c) ⟶ ∅, and (1) (0) lim  (t)≈ lim  (t)≈ 15.8114 t→∞ t→∞ X X (7d) ⟶ ∅, with the kinetic differential equations qualitatively different from the shallow maximum observed with the standard deviation of conventional chemical pro- da dx =− k ax +(a − a) r and = x (ka − r). cesses see Eq. (4h) . (7e) dt dt The irreversible processes,  + m → (m + 1) The simple first-order autocatalytic process in the flow reac- obtained from Eq. (5) by putting g = 0 , are illustrative tor sustains two long time states: (i) The state of extinction S because they are closer to biology where replication or with lim a(t)= a = a and lim x(t)= x = 0 , and (ii) reproduction occur always under the conditions of irre- t→∞ 0 t→∞ the quasistationary state S with lim a(t)= a = r∕k and versibility. The shape of the solution curves compared to 1 t→∞ Table 1 Fluctuations in the two Autocatalysis Reaction Initial conditions Limit autocatalytic reactions A+ X order ⇌ 2 X and A+2 X ⇌ 3 X X = 1 X = 2 X = 5 X = 10 t → ∞ 0 0 0 0 Zero A ⇌ X – – – – 15.8 First A+ X ⇌ 2 X 123.2 89.8 59.6 44.4 15.8 Second A+2 X ⇌ 3 X – 245.1 226.7 189.4 15.8 The table presents maximal standard deviations  computed from 1 000 trajectories of the two autocata- max lytic reactions for different initial conditions X = X(0) . For the autocatalytic reactions the standard devia- tion (t) passes through the maximum  listed here whereas for the uncatalyzed process it increases max monotonously from (0)= 0 to the equilibrium value (see Fig. 3) The equilibrium fluctuations calculated from equations  () are practically the same for all three reactions. −1 −1 −1 Choice of parameters, N = 1000 , A ⇌ X: h = h = 1.5[t ]; A+ X ⇌ 2 X: k = k = 0.01[M  t ]; A+2 → ← → ← −2 −1 X ⇌ 3 X: l = l = 0.00001[M  t ] → ← 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) The accurate value obtained from the stationary master equation is  = 15.8114 1 3 410 European Biophysics Journal (2018) 47:403–425 determine the probabilities to end up here or there. For suf- ficiently large population sizes, the long-time expectation values of the stochastic variables can be well approximated by linear combination of the deterministic values: (0) (1) (0) (1) E A =  a +  a andE X =  x +  x 0 1 0 1 with  = N ∕(N + N ) and  = N ∕(N + N ) where N 0 0 0 1 1 1 0 1 0 and N are the counted numbers of trajectories ending up in S or S from a sufficiently large sample. Although only 0 1 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, 1982, 1982). 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 extinc- tion ( s = 919). Fig. 4 The first-order autocatalytic reaction.  +  ⇌ 2 in the flow The major difference between the two classes of auto- reactor. The two plots show a single trajectory (top plot) and statis- catalytic reactions lies in the repertoire of possible dynamic tics consisting of expectation value within the one-standard deviation band, E A(t) ±  A(t) and E X(t) ±  X(t) taken form sample of behaviors. First-order autocatalysis gives rise to exponen- 100 trajectories (bottom plot). The four phases of the approach to the tial growth in unconstrained systems and to selection and long-time solution are indicated (top plot; see text). Choice of param- optimization of mean fitness in multispecies cases with −1 −1 eters and initial conditions: N = 2000 ; k = 0.01[M   t ], r = 0.5 −1 finite resources (see ’’Competition, mutation and qua- [V  t ], A(0)= 0 , X(0)= 1 (top) and X(0)= 3 (bottom). Pseudoran- dom number generator: ExtendedCA, Mathematica, seeds s = 491 . sispecies"). Accordingly first-order autocatalysis leads to a The sample size of the bottom plot was 100 trajectories, 3 led to S Darwinian scenario of selection of the fittest. In contrast (extinction) and 97 reached the pseudo-stationary state S . Solu- to first-order autocatalytic reaction networks, the dynamics tions of the corresponding kinetic differential equations are shown as in second-order systems is very rich and includes multiple dashed lines 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 lim x(t)= x = a − r∕k . Stability of S requires that the t→∞ 0 1 models, for example the Brusselator (Nicolis and Prigogine condition r < a k is fulfilled. 1977). It provides a simple enough reaction step for studies Starting from an empty reactor containing no A and by means of rigorous mathematics. Qualitative analysis of the autocatalyst only in seeding quantities, A(0)= 0 and Brusselator dynamics is straightforward and numerical inte- X(0)= 1, 2, 3, … , the trajectory shown in the upper part of gration causes no problem provided the integration software Fig. 4 allows for the distinction of several phases: (I) the can handle stiff differential equations. In reality, however, flow reactor is filled with the resource A in phase I, (II) in single-step autocatalytic reactions are extremely rare, instead the random phase II the decision is made to which state we are commonly dealing with multistep-reaction networks the trajectory will converge, (III) the trajectory approaches (Noyes et al. 1972) (see also the review by Francesc Sagués the long-time state, and (IV) the trajectory fluctuates and Irving Epstein (2003). around the state (see Figs. 4 and 7). First-order autoca- (0) In biology, in particular in the theory of evolution, the talysis sustains the two long-time solutions S =(a ,0) (1) (1) process A+2 X → 3 X plays a special role since in the sim- and S =(a , x ) , stochastic trajectories approach either ple form of hypercycles it is the basis for suppression of of the two states, and parameters and initial conditions 1 3 European Biophysics Journal (2018) 47:403–425 411 competitive selection without destroying template-induced reproduction. It provides one fundamental mechanism for major transitions (Maynard Smith and Szathmáry 1995; Schuster 1996) and will be discussed extensively in "Compe- tition and cooperation". The enormous flexibility of second- order autocatalysis follows, for example, from Fisher’s selec- tion 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 par- ents 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". Fig. 5 A minimal model for modeling evolution. Evolution is con- sidered 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 A minimal mathematical model parameters of all three processes, (i) fitness parameters f correspond- for the evolution of molecules ing to reaction rates for competition, (ii) reaction rates h for catalyzed reproduction, and (iii) an error rate parameter p for mutation are plot- ted on the axes of a Cartesian coordinate system. On the three two- The minimal model is dealing with the time dependence dimensional faces of the coordinate system we are dealing with the of the distribution of genotypes in populations Π(t) and three fundamental evolutionary processes: ( ) competitive reproduc- hence it is rooted in chemical kinetics of (i) competitive tion and mutation are the basis of Darwinian optimization through reproduction, (ii) symbiontic cooperation, and (iii) genetic natural selection and give rise to the formation of quasispecies and eventually to the occurrence of error thresholds (Eigen 1971; Eigen variation. The quantification of these three properties yields and Schuster 1977;  ) the interplay of competition and cooperation three parameters or sets of parameters, which can be plot- allows for the description of major transitions, which are seen as the ted on the three axes of a Cartesian coordinate system (see consequences of crossing resource thresholds (Maynard Smith and Fig. 5 and the previous publications (Schuster 2016a, b, d, Szathmáry 1995; Szathmáry and Maynard Smith 1995; Schuster 1996, 2016); ( ) the combination of cooperation and mutation ena- 2017b). We consider the simplest case here, where the three bles reintroduction of extinguished species provided the error rate is quantities are a fitness parameter f , a cooperation parameter sufficiently large such that a mutation threshold for stochastic survival h, and an mutation rate parameter p. For an implementa- can be recognized (Schuster 20160 tion 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 l Q i ji of a (hypothetical) compound A, A or a , respectively, and +  +  ⟶  +  +  , i, j = 1,… , n ; i mod n ; 0 0 i i+1 i j i+1 the mean residence time  = V∕r of the reaction mixture (8c) in the reactor where V is the volume of the reactor and r is (8d) ⟶ ∅ ; and the (volumetric) flow rate. The parameter  defines the time resolution of the reactor since slow reactions, which do not progress appreciably during the time interval  cannot be ⟶ ∅ , i = 1,… , n. (8e) R i studied. The process (8a) supplies the material required for repro- In the next step, the model is implemented by means of duction. A solution with A at concentration a flows into a a suitable and simple reaction mechanism. Based on "Ideal- continuously stirred tank reactor (CSTR) with a flow rate ized environments" and "Basic processes" we consider the 2 parameter r (Schmidt 2005, p. 87ff). The reactor operates at following set of 2n + n + 2 reactions in the flow reactor: constant volume and this implies that the volume of solu- a r (8a) ∗ ⟶  ; tion flowing into the reactor per time unit [ t] is compensated exactly by an outflow, which is described by the Eqs. (8d) k Q i ji and (8e) and concerns all (n + 1) molecular species, A and (8b) +  ⟶  +  , i, j = 1,… , n ; i i j , i = 1,… , n . Inflow and outflow are often characterized 1 3 412 European Biophysics Journal (2018) 47:403–425 as pseudoreactions because they are no chemical reactions long as n is not too large. In absence of cooperation terms, in the strict sense, which are converting reactants are into l = 0∀ i = 1,… , n , Eq.  (9) can be transformed into an products. The two classes of reactions producing progeny, eigenvalue problem of a symmetric matrix, which is readily template induced replication (8b) and catalyzed template diagonalized provided n is not very large ( n < 10 ; "Com- induced replication (8c), represent the core of the evolu- petition, mutation and quasispecies"). In mutation-free tion model. In agreement with the conditions in biology, systems, p = 0 ("Competition and cooperation"), qualitative both reproduction steps are implemented irreversibly in the analysis and determination of stationary states are straight- direction of polynucleotide synthesis. A basic assumption forward, and the dynamics of the complete system can be for both reproduction steps is that correct reproduction and derived by extrapolation from the error-free results to finite mutation are parallel chemical reaction channels ("Basic mutations rates. The cooperation system with mutation, (9) processes"). In other words, there is no mutation under con- with k = 0∀ i = 1,… , n is used here to study the relevance ditions that do not sustain reproduction. of mutation in symbiontic systems. It also serves as an exam- As an alternative to the Eigen model (1) mutation can ple for the study of unconventional consequences of repli- be seen, for example, as the result of DNA damage and cation with frequent errors in the strong mutation scenario imperfect damage repair during the whole life span of an ("Cooperation and mutation"). organism, which is the idea underlying the Crow–Kimura mutation model (Crow and Kimura 2009, pp  264–266). Processes on individual coordinate axes Then reproduction and mutation are completely independ- ent processes, The processes along the individual coordinate axes are considered in order to verify the initial statement on ji the three basic processes. For this goal it is easiest to set +  ⟶ 2 and  ⟶  , (8b’) i i i j certain parameters zero: (I) no natural selection implies and in the kinetic differential equations they appear as addi- k = k =…= k = k (= 0) , (II) no cooperative coupling 1 2 n tive terms. Interestingly, the Eigen and the Crow–Kimura requires l = l =…= l = l = 0 , and (III) no mutation 1 2 n model although being fundamentally different with respect leads to Q =  where  is the unit matrix. to the assumptions about the nature of the mutation pro- A process taking place on the selection axis is given by cess give rise to the same mathematical problems [see, e.g., (II) and (III) being true leads to the ODE (Schuster 2016d, pp.76-78)]. The mutation matrix Q cor- responds formally to the mutation matrix  but there are da =−a k x + r (a − a) and (10a) i i 0 non-negligible differences Q covers correct and error-prone dt i=1 replication but the process  → 2 is handled separately i i in the Crow–Kimura model and hence all diagonal terms are dx zero μ = 0∀ i = 1,… , n. ii =(k a − r) x ; i = 1, 2, … , n . (10b) i i dt The equations that will be applied in the analysis of the dynamics of the model (8) implement three processes along Competition between the reproducing elements the coordinate axes: (i) Darwinian selection of the fittest leads to survival of the fittest subspecies  , which is on the competition axis, (ii) hypercycle dynamics on the defined by k a = f = max{f ;i = 1, … , n} . For constant m m i cooperation axis, and (iii) neutral evolution on the mutation []= a = const , the mean fitness of the population, ∑ ∑ n n axis. The kinetic differential equations of the model mecha- (t)= f x (t) c(t) with c(t)= x (t) , is a non- i i i i=1 i=1 nism (8) are of the form: decreasing function of time: d∕ dt = var{f } ≥ 0 . This result is the formal analogue of Fisher’s fundamental theorem for da asexual reproduction. =−a k + l x x + r (a − a) , i mod n and i i i+1 i 0 dt i=1 (9a) We remark that the deterministic kinetic Eqs. for (8b) and (8c) dx were extensively studied under the simplifying assumption of = a Q k + l x x − x r , n ij j j j+1 j i constant population size x (t)= c = const (Eigen 1971; i 0 dt i=1 (9b) j=1 Eigen and Schuster 1978; Eigen and McCaskill 1989). The deter- ministic solution curves formulated in relative concentrations i = 1, 2,… , n; j mod n . ξ (t)= x (t) x (t) are identical for the CSTR and for constant i i i i=1 population size (Schuster and Sigmund 1985). As we mentioned In  (9a) we made use of the conservation relation already in "Idealized environments" the stochastic system is unstable Q = 1 . No analytical solutions are available for  (9) for unregulated constant population size (Jones and Leung 1981) and ij i=1 a proper description has to account explicitly for the physical setup in general but numerical integration is straightforward as applied, here the flow reactor. 1 3 European Biophysics Journal (2018) 47:403–425 413 On the cooperation axis, the conditions (I) and (III) are assumption that chemical processes occur one at a time, all fulfilled and we obtain the equations of hypercycle dynamics jumps involve single steps and the particle numbers change by ±1 unless the elementary steps involve more than a single da molecule per class. The jumps S  → S or S → S  are =−a l x x + r (a − a) and (11a) i i i+1 0 denoted by the shorthand notation dt i=1 =(m ± 1, s , … , s ) ≡ (; m ± 1) or 1 n dx i � =(l x a − r) x , i = 1, 2, … , n, i mod n ,  =(m, s , … , s ± 1, … , s ) ≡ (; s ± 1). (11b) i i+1 i 1 k n k dt Then the master equation of mechanism (8) takes on the which were studied in great detail in a number of previous form papers (Eigen and Schuster 1978; Schuster 1978; Schuster et al. 1979, 1980; Hofbauer et al. 1991). dP The third case—with (I) and (II) being true—yields a = a r P − P + r (m + 1)P − mP 0 (;m−1)  (;m+1) dt degenerate deterministic solution: All distributions of sub- species with fixed population size x = c yield the same i=1 + r (s + 1)P − s P i (;s +1) i solutions and therefore constitute an (n − 1)-dimensional i=1 invariant manifold fulfilling + Q (k + l s ) (m + 1)(s − 1)P − ms P ii i i i+1 i (;m+1,s −1) i da i=1 =−ak c + r (a − a) and (12a) n n dt + Q (k + l s )s (m + 1)P − mP . ij j j j+1 j (;m+1,s −1) i=1 j=1,j≠i dc = ak c − cr =(ak − r) c and (14) dt (12b) =(a − ) r with  = a + c dt Each reaction step involving S changes the probability to be in state S , P , in two ways: It increases the prob- In the absence of selection and cooperation neutral evolu- ability through reactions or pseudoreactions S → S tion in the sense of Motoo Kimura (1955, 1983) is observed and decreases the probability through the reaction steps on the mutation axis. For an adequate description of the S → S where  is summed over all states from which process a stochastic treatment is required. Random selection S can be reached or vice versa. The two terms in the first of a single arbitrarily chosen subspecies is observed in the line, for example, describe the two pseudoreactions mod- no-mutation limit, lim p → 0 . For non-vanishing mutation eling inflow and outflow of the material A , and further each rates once selected subspecies may be replaced by other sub- reaction is represented by two steps. It is also worth noticing species and the mean time of replacement of one randomly that stoichiometry requires two slightly different replication −1 selected subspecies is T =  where  is the mutation rep terms depending on whether the copy is correct or incorrect. rate per generation (Kimura 1955, 1983). Translated into our Master equations are easily written down and station- model, we find for single point mutations  ≈ p. ary solutions can be derived by generally applicable techniques as it was shown for the flow equilibrium in Master equation and simulation "Deterministic and stochastic approaches" but explicit time-dependent solutions are very hard to obtain and Reaction  (8) can be cast into chemical master equations. known only in exceptionally simple cases [Schuster 2016e, The particle numbers of the molecular species, []= A(t) pp 347–568]. Here, we shall use the simulation technique and [ ]= X (t) with i = 1,… , n , are integers and in the i i of sampling trajectories introduced already in "Autocataly- absence of flows they fulfil the conservation relations, sis in the batch reactor". Expectation values and second A(t)+ X (t)= C . The variables of the master equation i=1 moments of variables can be computed through sampling are the probabilities of trajectories—an example is shown in Fig. 6—but often this approach exceeds the available computational facili- P (t)= Prob A(t)= m and m ties and it is necessary to interpret single trajectories. As examples we consider single trajectories in Fig.  4 and (13) P (t)= Prob X (t)= s ; i = 1, … , n , s i i in Fig.  7. The process for convenience starting from an empty flow reactor is split into four phases: (i) establish- and the indices are subsumed in an index vector, ment of the flow equilibrium of A , (ii) random decision =(m, s , … , s ) , which characterizes the state S of 1 n on the (quasi)stationary state towards which the trajectory the system. The chemical master equation is based on the 1 3 414 European Biophysics Journal (2018) 47:403–425 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 pro- cess is assumed to start with an empty reactor except seeds for the two autocatalysts  and  . It can be partitioned into four phases: 1 2 (I) fast raise in the concentration of A, (II) a random phase where (1) (2) the decision is made into which final state— S , S , S or S —the 0 2 1 1 trajectory progresses, (III) the approach towards the final state— here S —and (IV) fluctuations around the values of the (quasi) stationary state. Comparison with the simple autocatalytic pro- cess in Fig.  4 reveals great similarity. Choice of parameter values: −1 −1 −2 −1 k = 0.099 , k = 0.101  [M t ], l = 0.0050 , l = 0.0045  [M t ], 1 2 1 2 −1 Fig. 6 Quasispecies formation. The two plots show the forma- a = 200 , r = 4.0  [V  t ], pseudorandom number generator: Extend- tion of quasispecies from an initially empty reactor, A(0)= 0 , with edCA, Mathematica, seeds s = 631 . Initial conditions: A(0)= 0 , replicators in seeding amounts. The system in the upper plot was X (0)= X (0)= 1 . Color code: A(t) black, X (t) red, and X (t) yellow. 1 2 1 2 initiated to be a single copy of the master sequence, X (0)= 1 and The figure is adapted from Schuster (2017a) X (0)= X (0)= X (0)= 0 . At the beginning the system might be 2 3 4 extinguished by a single dilution event X → ∅ and the high prob- ability of extinction gives rise to enormously broad and overlapping fittest subspecies and at non-zero mutation rates this sce- one-standard-deviation bands. The initial values in the lower plot nario yields selection of a fittest ensemble of subspecies, were chosen such that the probability of extinction is zero for all prac- tical purposes: X (0)= 10 and X (0)= X (0)= X (0)= 0 . Choice which has been characterized as quasispecies (Eigen and 1 2 3 4 −1 −1 of parameters: N = 2000 ; k = 0.011[M   t ], k = k = 0.010 1 2 3 Schuster 1977; Domingo and Schuster 2016). Precisely, the −1 −1 −1 −1 −1 [M   t ], k = 0.009[M   t ], r = 0.5[V  t ], Color code: A black, quasispecies is the stable stationary long-time distribution red,  yellow,  green, and  blue. Expectation values are 1 2 3 4 of a population of subspecies undergoing replication and shown as full lines, deterministic solutions as broken lines. Since x (t)= x (t) is fulfilled for the parameter values used here the curve is mutation: A population that consists of several genotypes 2 3 shown as a yellow–green broken line present in time-dependent concentrations, Π(t)= x (t) ⊕ x (t) ⊕ … ⊕ x (t) (15) 1 1 2 2 n n has the stationary solution, lim Π(t)=  = x  ⊕ x  ⊕ … ⊕ x  , which is t→∞ 1 1 2 2 n n converges, (iii) the approach towards this (quasi)station- called the quasispecies . ary state and (iv) fluctuations around the (quasi)stationary Deterministic quasispecies The deterministic or continu- state. The separation into phases is made possible by the ous quasispecies represents the unique deterministic long- choice of suitable initial conditions. time solution of the replication mutation problem, which is described in the flow reactor by the ODE (Schuster and Competition, mutation and quasispecies Sigmund 1985): The bottom face of the three-dimensional parameter space— da  in Fig. 5—is dealing with selection provided by the com- =−a k x + r(a − a) (16a) i i 0 dt bination of competition and mutation. Natural selection at i=1 zero mutation rate leads to a homogeneous population of the 1 3 European Biophysics Journal (2018) 47:403–425 415 an approximate uniform distribution. An explanation is dx i straightforward: Above this critical transition, mutations = a Q k x − rx , i = 1,… , n . (16b) ij j j i dt occur too often to sustain sufficiently accurate reproduction j=1 of the template sequence and the result is random replica- tion: In the long run, every sequence is obtained with the All early works on quasispecies dynamics were per- same probability. The transition has been characterized as formed with the constraint of constant population size: error threshold (Eigen 1971; Eigen and Schuster 1977) since x (t)= c = const with a(t)= a and f = k a that gives i 0 i i 0 i=1 evolutionary dynamics does not sustain a structured long- rise to the differential equation time population at higher error rates, p > p . On typical cr dx fitness landscapes, the error threshold sharpens with increas- = Q f x − (t) x , ij j j i ing chain length l and becomes a first-order phase transition dt j=1 in the limit l → ∞ (Tarazona 1992; Park et al. 2010; Huang f x 7 j j j=1 et al. 2017). with  = ; i = 1, … , n , (16’) Discrete quasispecies In the continuous description, the j=1 quasispecies contains all species at finite positive concentra- which fulfils dc∕ dt = dx ∕ dt = 0 .   Equations (16) i=1 tions no matter how small the concentrations might be. Deal- and (16’) have identical solutions in normalized variables ing with less than one molecule per reactor volume, how- (t)= x ∕ x (t) but the stability properties of the corre- i i i i=1 ever, is unrealistic. Numbers of molecules are non-negative sponding stochastic systems are different. [For more details integers, X ∈ ℕ , subspecies distributions are truncated in see (Thompson 1974; Jones et al. 1976; Ebeling and Mahnke reality and we call them stochastic or discrete quasispecies: 1979; Jones and Leung 1981; Schuster and Sigmund 1985)]. The quasispecies as a function of the mutation rate param- ⌊x ⌋ if x ≥ 1 i i = X  ⊕ X  ⊕ … ⊕ X  with X = . 1 1 2 2 n n i eter (p) begins at as a homogeneous population con- 0 if x < 1 p = 0 taining exclusively the fittest subspecies  and becomes a (17) distribution of subspecies at nonzero p. This distribution The discreteness of the stochastic variables leads to a consists of a most frequent sequence called master sequence, modification of the scenario for the mutation dependence which is surrounded by a mutant cloud. Commonly, the most of quasispecies (p) . In the mutation-free system, p = 0 , frequent or master sequence is also the fittest one, but this survival of the fittest is observed and the quasispecies con- is not necessarily so: In case of strong stationary mutation sists of a single sequence:  = X  = N with m indi- m m m flow from the mutant cloud to the master sequence, a less cating maximal fitness, f = max{f ;i = 1, … , n} , a nd N m i fit master may outgrow a fitter sequence with less efficient being the population size. At sufficiently small values of the mutational backflow. At further increase in p the distribution mutation rate parameter p, no subspecies except the master broadens, and eventually ends up as the uniform distribu- exceeds the threshold x ≥ 1 and the discrete quasispecies tion  ∶ { x =  x = ⋯ =  x } at  p = 1 −( − 1) p∕ = 1∕ still consists of a single fittest genotype  . The scenario 1 2 n l l l where n =  and  x = 1∕ (i = 1, …  ) for sequences of in the small mutation rate regime—populations are almost chain length l over an alphabet with  letters. At the muta- always homogeneous except short non-stationary periods tion rate , correct and wrong digits are incorporated p =  p during which advantageous mutations take over the popula- with the same frequency. The uniform distribution is the tion—is tantamount to the strong selection–weak mutation dynamical answer to the absence of any preference for nucle- scenario discussed in evolutionary biology (Gillespie 1983; otide assignment: All subspecies in the stationary distribu- Joyce et al. 2008). If the p-value is so large that for one or tion occur with the same frequency. The quasispecies in the intermediate range is determined by the fitness landscape— We mention that some model landscapes like the additive land- the distribution of fitness values f in sequence space, by scape or the multiplicative landscape sustain smooth rather than sharp transitions (Wiehe 1997). For a recent update and a review of the the move set of allowed mutations as well as the mutation state of the art in the relation between fitness landscape and selection rate p. Typically, sharp transitions occur at some critical dynamics see (Schuster 2016d). mutation rates p = p : The distribution changes smoothly tr Biologists (Dean and Thronton 2007; Sniegowski and Gerrish from to , where the distribution turns abruptly p = 0 p = p tr 2010) and computer scientists commonly distinguish strong and weak into another distribution. At the transition with the larg- mutation. The weak mutation scenario assumes that adaptive muta- tions are sufficiently rare and do not interfere with the selection pro- est p-value , the quasispecies becomes p = 1∕𝜅> p > p cr tr cess but initiate replacements of currently fittest genotypes by still fitter variants. The strong mutation scenario is characterized by suf- ficiently large values of p that give rise to the quasispecies dynam- ics described here [for details see (Domingo and Schuster 2016)] and The value  = 2 is used here and applies for binary sequences. For to mutation induced cooperative dynamics ("Cooperation and muta- the natural AUGC -nucleotide alphabet we have  = 4. tion"). 1 3 416 European Biophysics Journal (2018) 47:403–425 more subspecies x ≥ 1 is fulfilled, selection of a family of Fluctuations at small particle numbers have different ori- genotypes is observed:  consist of a master genotype  , gins ("Autocatalysis in the batch reactor"): (i) conventional together with a stationary distribution of sufficiently frequent thermal fluctuations, (ii) enhanced fluctuations related to mutants  (i ≠ m ). The quasispecies becomes broader with autocatalytic self-enhancement, and (iii) anomalous fluctua- increasing mutation parameter p until a threshold value p tions in the stochastic variables arising from two or more cr is reached above which error propagation does not sustain quasistationary states. The standard deviations (t) fulfil the a stationary state and the population Π drifts randomly N-law for the resource A but are larger for the autocata- through sequence space (Huynen et al. 1996; Higgs and Der- lysts  ,  ,  , and  . The stationary states of the stochas- 1 2 3 4 rida 1991). Interestingly, the critical mutation rate p can be tic system are extinction S and quasispecies S . Fig. 6 shows cr 0 1 derived from the continuous quasispecies theory. a typical example: The anomalous fluctuations in the upper As calculations of the time dependencies of the first two plot are in full analogy to first-order autocatalysis (Fig.  4). moments of the probability distribution of subspecies in the Since the initial condition X (0)= 1 was chosen, one out- (Π) population, P (t) , reveal (Jones and Leung 1981)  (16’) is flow step may extinguish the population and the probabil- only marginally stable: ity of dying out is indeed as large as 20%. The fluctuations bands are extremely broad, and large differences between � � d⟨N⟩ the deterministic solution and the corresponding expecta- = f ⟨X ⟩ − ⟨ (t) N⟩ = 0, (18a) i i tion values are observed. In the lower plot initial conditions dt i=1 X (0)= 10 were chosen, which are sufficient to reduce the contributions of anomalous fluctuations practically to zero. � � d⟨N ⟩ Then, for a population size of N = 2000 the concentration = 2 f ⟨X N⟩ − 2⟨ (t) N ⟩ + i i dt a(t) coincides with the expectation value E A(t) for all prac- i=1 (18b) n n tical purposes and the fluctuations fall into a typical ± N � � d⟨N⟩ + 2 f ⟨X ⟩ − = 2 f ⟨X ⟩ , and -band. The autocatalysts, X ,… , X , show broader than i i i i 1 4 dt i=1 i=1 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 2 2 var(N)= ⟨N ⟩ − ⟨N⟩ = 2 f ⟨X ()⟩ d (18c) i i an open system and does not approach an equilibrium state i=1 (see also Fig. 4). The time derivative of the first moment vanishes as expected It is worth recalling what means stochasticity for qua- since the condition of a constant population size is fulfilled sispecies: (i) continuous concentrations are replaced by by the differential equation (16’). The variance, however, discrete particle numbers, (ii) fluctuations replace single increases with time since the integrand in (18b) is always line trajectories by bands within which trajectories follow a positive. After sufficiently long time, the fluctuation band probability distribution, (iii) subspecies can be diluted out becomes so large that the expectation value is irrelevant for of the flow reactor and if this happens for all subspecies the the description of the system. The instability in Eigen’s qua- population goes extinct giving rise to anomalous fluctua- sispecies equation (Eigen 1971) is well known from similar tions, and (iv) error thresholds introduce random reproduc- problems: (i) the neutral birth-and-death process with equal tion that is closely related to Motoo Kimura’s random drift. birth and death parameters,  =  , and (ii) the Wiener pro- An increase in the error rate up to the error threshold leads to cess. In both cases a constant expectation value is jeopard- broadening of the mutant spectrum surrounding the master ized by a variance that grows linearly with time. Stability sequence. Above the thresholds, the populations migrate by against fluctuation is easily introduced into (16’): one needs random drift through sequence space. only to give up the condition of strictly constant population size and to replace the denominator in (t) by a constant Θ, Competition and cooperation The kinetic equations for replication describing the template (t)= f x Θ , j j induced, uncatalyzed and catalyzed processes are obtained j=1 from Eq. () by neglect of mutation. Then the kinetic differ - ential equations for competition and cooperation of competi- where Θ is the population size towards which the population tors result from (9) by setting Q =  : converges after sufficiently long time. The approach cho- sen here, the implementation of a flow reactor as a physical da device, yields a stable system as well. =−a k + l x x + r (a − a) and (19a) i i i+1 i 0 dt i=1 1 3 European Biophysics Journal (2018) 47:403–425 417 Table 2 Asymptotically stable Symbol Stationary values Stability range stationary states of Eq. (19) with n = 3Schuster (2016) a x x x 1 2 3 S a 0 0 0 0 ≤ a ≤ 0 0 0 3 (3) 0 0 a −   ≤ a ≤  + S 3 0 3 3 0 3 32 (1) 0 a −  −    +  ≤ a ≤  +  + 3 0 3 32 32 3 32 0 3 32 31 S      +  +  ≤ a 3 3 1 2 3 32 31 0 The four stationary states are ordered with respect to increasing a -values of their asymptotically stable regime. The relations k < k < k and l > l > l between the rate parameters were assumed. The abbre- 1 2 3 1 2 3 viations  = r∕k ,  =(k − k )∕l ) , and  =(k − k )∕l ) are used for the combination of param- 3 3 31 3 1 1 32 3 2 2 eters. For the cooperative state S the stationary concentration of A is obtained as one root of a quad- ∑ ∑ 3 3 ratic equation (20) with two combinations of the rate parameters,  = k ∕l and  = 1∕l , and i i i i=1 i=1 =(r − k )∕(l ) for i = 1, 2, 3 and  = a from   (20). The existence of the non-trivial stationary state i i i 1 requires a sufficiently small flow rate: r ≤ (a + ) ∕4 a and is stable for a < r∕k , the state of selection of  , 0 0 3 3 dx i (3) S , is stable in the range r∕k < a < r∕k +(k − k )∕l , = a k + l x − r x ; i = 1… , n; i mod n . 3 0 3 3 2 2 i i i+1 i 1 dt (1) followed by ex c l u s i o n of  , S , in the range (19b) k +(k − k )∕l < a < k +(k − k )∕l +(k − k )∕l . 3 3 2 2 0 3 3 2 2 3 1 1 Both equations contain terms of different molecularity in Above this value for a cooperation of all three sub- and this has the consequence that the dynamic behav- ∑ species is observed provided the flow rate is not too ior depends strongly on the total concentration c = x , i 2 i=1 large: r < r =(a + 𝜓 ) ∕4𝜙 with  = k ∕l and cr 0 i i i=1 which in turn is determined by the amount of available = 1∕l . The free concentration of A is obtained as i=1 resources, a : At sufficiently low concentration, the first- 0 9 solution of a quadratic equation order terms, k x , dominate whereas the second-order terms, i i l x x become important at high concentrations. No direct 1 i i+1 i a = a +  ∓ (a + ) − 4r , (20) 1,2 0 0 analytical solutions are available but exhaustive qualitative analysis is possible and the concentrations of the stationary where the minus sign corresponds to the cooperative state S . solutions, a, x (i = 1, … , n) , can be computed analytically The second solution belongs to an unstable state S , which Schuster (2016a, d) (Table 2). 3 separates the basins of attraction of the states of coopera- The basis of the calculations of stationary states tion and extinction. Above the critical flow rate, r > r the cr is the factorability of the r.h.s. of (19b). The relation states S and S do not exist. For n > 3 , the situation becomes x ⋅ (k + l x )a − r = 0 with i mod n is compatible with 3 i i i i+1 more complex since solutions may oscillate. Many systems stationary states that fulfil either with n = 4 have oscillations with very weak damping fac- 1 r (i) (ii) tors, n ≥ 5 commonly leads to undamped oscillations. The x = 0 or x = − k , ∀ i = 1, … , n , i−1 i i l a i−1 properties of these systems have been discussed extensively n in previous publications to which we refer here (Schuster and the conservation condition a + x = a . In total i 0 i=1 2016e; Schuster and Sigmund 1985; Schuster 1978). there are 2 possible states out of which a single one is Stochasticity has common effects on competition–coop- asymptotically stable for a given set of parameters. In eration systems like thermal fluctuations and fluctuation Table 2 the results for n = 3 are shown. The number of sub- through autocatalytic enhancement. In addition, there are species present at the stationary state, n , is suitable for char- many more quasistationary states than the asymptotically acterizing the state: n = 0 is the state of extinction, n = 1 S S stable states of the deterministic system. For example, states is the state of selection, n = 2 is the state of exclusion, in which less efficient subspecies are selected show up as etc., and finally n = 3 = n occurs at the state of coopera- quasistationary states as well. In case of the smallest pos- tion where all three subspecies are present. The numbers of sible system with n = 2 and k > k , we have four states: 2 1 long-time subspecies n depend on the resource input a . S 0 (i) the absorbing boundary as state of extinction S , (ii) the The relative size ordering of the parameters determines the (2) state of natural selection S where the fittest variant  is identity of the selected, excluded, etc., subspecies. For the calculations shown in Table 2 the orderings k < k < k and 1 2 3 l > l > l were applied and hence  is selected and  is 1 2 3 3 1 It is worth mentioning that Eq.  (20) and the solutions for a are 1,2 excluded. Considering steady states as functions of the avail- valid for arbitrary n provided the summations in  and  are adjusted able resources, the state of extinction S comes first at small accordingly. 1 3 418 European Biophysics Journal (2018) 47:403–425 Table 3 Probabilities to reach Initial values Counted numbers of states in final outcomes quasistationary states in the cooperative regime with n = 2 X (0) X (0) N N (1) N (2) N 1 2 S S S S 0 2 1 1 and different initial conditions 1 1 385.1 ± 23.6 1481.0 ± 36.8 1719.6 ± 37.8 6414.3 ± 53.8 2 1 77.4 ± 9.1 1822.6 ± 41.6 367.6 ± 17.0 7733.3 ± 38.3 1 2 71.6 ± 8.5 280.6 ± 20.0 2075.8 ± 28.9 7572.0 ± 39.2 3 1 15.0 ± 2.9 1900.4 ± 30.9 74.69 ± 10.0 8009.0 ± 35.3 1 3 14.0 ± 3.7 53.1 ± 4.8 2180.5 ± 48.4 7752.3 ± 53.8 2 2 14.9 ± 2.6 303.7 ± 16.0 354.5 ± 23.8 9326.8 ± 44.9 3 3 0 70.2 ± 10.0 106.2 ± 10.9 9823.4 ± 15.7 4 4 0 12.1 ± 2.6 28.0 ± 5.0 9959.9 ± 6.4 5 5 0 2.5 ± 1.1 6.3 ± 2.6 9991.2 ± 3.0 The table provides probabilities of occurrence for all four possible long-term states: extinction S , selec- (1) (2) tion of  in S , selection of  in S , or cooperation S . The counted numbers of events are sample 1 2 2 1 1 means and unbiased standard deviations calculated from ten packages, each of them containing 10 000 tra- jectories 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, Mathemat- −1 −1 −1 −1 −2 −1 −2 −1 ica). Choice of parameters: k = 0.09[M t ], k = 0.11[M t ], l = 0.0011[M t ], l = 0.0009[M t ], 1 2 1 2 −1 −4 a = 200 , r = 0.5[V t ]. Initial value A(0)= 0 . Probabilities are obtained by multiplication by 10 Table 4 Long-time behavior Mutation Initial values Expectation values at t = 30 Counts in the cooperation-mutation system with n = 3 and different p X (0) X (0) X (0) E(A) E(X ) E(X ) E(X ) P(S ) 1 2 3 1 2 3 3 initial conditions 0 1 1 1 279.3 44.56 36.06 39.58 0.301 0 2 2 2 81.53 117.7 94.9 105.4 0.814 0 4 4 4 2.023 146.1 119.6 132.2 0.996 0 10 10 10 0.376 147.0 120.0 132.6 1.0 0 Deterministic 0.377 147.0 120.3 132.3 1 0.05 1 1 1 177.2 81.66 68.16 72.64 0.432 0.05 2 2 2 28.57 136.1 113.8 121.5 0.937 0.05 4 4 4 0.383 147.0 122.2 130.4 1.0 0.05 Deterministic 0.377 146.7 122.3 130.5 1 0.1 1 1 1 161.6 86.95 74.39 76.96 0.599 0.1 2 2 2 28.57 136.1 113.8 121.5 0.957 0.1 4 4 4 0.383 147.0 122.2 130.4 1.0 0.1 Deterministic 0.377 146.7 122.3 130.5 1 The table provides expectation values at the time of the end of the simulation (t = 30 ) for different end mutation rate parameters p and different initial conditions. Choice of parameters: l = 0.011 , l = 0.010 , 1 2 −2 −1 −1 l = 0.009[M t ], a = 400 , r = 0.5[V t ]. Initial value A(0)= 0 3 0 selected, (iii) the state of selection of the less fit subspecies solution. In case of n = 3 we present expectation values of , and eventually (iv) the state of cooperation S . The rela- the four stochastic variables, A(t) and X (t) ( i = 1, 2, 3 ) at 1 2 i tive stabilities of the individual states are reflected by the some predefined time of the simulation end, t for differ - end probabilities to reach these states by randomly chosen trajec- ent initial conditions and mutation rate parameters p. The tories (Table 3). Parameters for the calculations shown in the values for p = 0.0 refer to the pure competition–coopera- table were chosen such that the corresponding deterministic tion case discussed here (Table 4). As expected extinction system is situated in the cooperative domain in parameter plays a major role at small initial values, X (0)= 1, 2, 3 space, and indeed the approach towards the cooperative state ( i = 1, 2, 3 ), but X (0)= 4 is already sufficient for coming S has always the largest probability. Again, initial particle close to the expectation values obtained for large initial numbers around X (0)+ X (0)= 4 are sufficient for strong values, X (0)= 10 is enough for reaching the deterministic 1 2 i dominance of the state corresponding to the deterministic values for practical purposes. 1 3 European Biophysics Journal (2018) 47:403–425 419 In systems with n ≥ 4 deterministic and stochastic solu- tion curves oscillate. The solutions of the ODE’s are differ - ent for n = 4 where weakly damped oscillations occur and for n ≥ 5 showing undamped relaxation oscillations (Phil- lipson 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 domi- nant 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 param- eter p was chosen and where accordingly quasispecies for- mation 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 diver- Fig. 8 Times to extinction as a function of available resources in the sity through variation. In principle, mutation can reintroduce five membered cooperative system ( n = 5 ). Extinction times T of the extinguished subspecies into the population. Here, we shall populations Π are shown for different amounts of available resources focus on this aspect and, in particular, study the influence of measured as inflow concentrations a or A when expressed in num- 0 0 the mutation rate parameter p on extinction times. To study bers of molecules per unit volume. The upper diagram presents the data at a resolution of ten molecules (Δ A = 10 ; 100, 110, 120, ) for the role of mutation in low-dimensional cooperative systems 0 four different values of the mutation rate parameter: p = 0.0000 (red), ( n = 2, 3 ) expectation values of the stochastic variables A(t), p = 0.0005 (yellow), p = 0.0010 (green), and p = 0.0020 (blue). T X (t) and X (t) , or X (t) and X (t) were calculated at prede- 1 2 2 3 -values above 1000 are truncated at this value. The lower diagram fined times ( t ) and compared with the probabilities of tra- shows the two plots p = 0.0000 (red) and p = 0.0020 (blue) at the end highest possible resolution (Δ A = 1 ). Choice of other parameters: jectories to end up in the cooperative state, S or S , respec- 0 2 3 −1 −1 −2 −1 r = 0.5[V t ], l = l = l = l = l = 0.01[M t ]. Pseudorandom 1 2 3 4 5 tively. For the initial conditions X (0)= X (0) = X (0) > 4 1 2 3 number generator: Extended CA, Mathematica, seed: s = 491 . Initial the results are practically indistinguishable from the deter- conditions: A(0)= 0 , X (0)= X (0)= X (0)= X (0)= X (0)= 5 1 2 3 4 5 ministic values. An increase in the mutation rate parameter p shows the expected inu fl ence: Extinguished subspecies can be reintroduced and this increases the probability of reaching are less time consuming and provide in essence the same the cooperative state, S or S . The case n = 3 is shown as an results. The plots shown in Fig. 8 show enormous scatter but, 2 3 example in Table 4. nevertheless, allow for drawing two conclusions: (i) In the The oscillating systems are more difficult to investigate. mutation-free case the extinction time T is independent of Here, we consider the time of extinction of the entire popu- the amount of available resources up to a value A ≈ 1300 , lation as a function of the mutation rate and the available and (ii) for non-zero mutation rates a kind of noisy or sto- resources, T (p, A ) . The results are shown in Fig  8: The 0 0 chastic threshold phenomenon is observed. Considering the extinction times T show very strong scatter and their appear- noisy function T (p, A ) and taking A at the first value of 0 0 0 ance is dependent on the resolution of the calculations. By T (p, A ) ≥ 1000 we find for the parameter values applied in 0 0 (T ≥1000) (cr) resolution, we mean here the number of molecules A in A Fig. 8: A = A = 1130, 690, 360 for the mutation 0 0 between two neighboring points. The highest resolution is rates p = 0.0005, 0.0010, 0.0020 , respectively. As expected (cr) achieved when the calculations are performed with every the threshold moves to lower A -values with increasing (integer) number of A molecules, e.g., ΔA = 1 yields 100, 0 0 mutation rate. The behavior of the extinction times T (A ) is 0 0 101, 102,  . Computations at somewhat lower resolutions 1 3 420 European Biophysics Journal (2018) 47:403–425 (cr) similar for n = 4 but the critical concentrations for the different p-values lie much closer together and the analysis is more difficult. Considering survival at constant resources A 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 rever- tant 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 (Albery and 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 strat- egy 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 coop- erative reproduction, or with the introduction of cooperation into Darwinian systems. The characteristic dependence of the population dynam- ics 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 dis- cussed in "Competition and cooperation". An increase of the Fig. 9 The transitions from competition to cooperation in the sys- cooperation parameter h = lA(t) leads in steps from selec- tem with n = 5 The transitions from competition to cooperation in the system with n = 5 . The three plots show single trajectories for tion of the fittest to a cooperative state with all subspecies three different scenarios in the flow reactor: (i, topmost plot) the qua- present. Oscillating systems ( n ≥ 4 ) are more spectacular sispecies scenario with a dominant master sequence (  ), (ii, mid- since the hypercycles are unstable at p = 0 in the stochas- dle plot) an intermediate scenario with irregular dynamics and two tic approach and raising the cooperation parameter h leads dominating species (  and  ), and (iii, bottom plot) the stochastic 1 5 hypercycle scenario with irregular, undamped oscillations. Choice from selection to extinction. In the intermediate param- −1 of parameter: k = 0.150 , k = k = 0.125 , k = k = 0.100 [M 1 2 5 3 4 eter range where the deterministic system shows stepwise −1 t ], l = l = l = l = l = l , l = 0.0 (topmost plot), l = 0.002 1 2 3 4 5 −2 −1 increase in the number of coexisting subspecies (1, 2,… , n (middle plot), l = 0.01  [M t ] (bottom plot), a = 800 , −1 or, expressed phenomenologically, selection, exclusion, , r = 0.5 [V t ], p = 0.075 , pseudorandom number generator: Extend- edCA, Mathematica, seeds s = 491 . Initial conditions: A(0)= 0 , cooperation) the stochastic approach yields highly irregular X (0)= X (0)= X (0)= X (0)= X (0)= 5 . Color code: (t) black, 1 2 3 4 5 dynamics with different numbers of non-extinguished sub- (t) red,  (t) yellow,  (t) green,  (t) blue, and  (t) cyan 1 2 2 2 5 species whereby the number of species present increases 1 3 European Biophysics Journal (2018) 47:403–425 421 with increasing values of the cooperation parameter h. The of asexually reproducing genotypes in a population of con- scenario in parameter space is completed by considering stant size N. Correct replication and mutation are seen as the series of states at different mutation rate parameters p . parallel chemical reactions leading to a uniquely den fi ed sta - The case n = 5 , which for large h-values leads to undamped tionary population called quasispecies (Eigen and Schuster oscillations with a stochastic contribution to the amplitude, 1977). RNA replication catalyzed by single virus specific was chosen to facilitate the illustration (Fig.  9): At zero enzymes from RNA bacteriophages provides a bridge from mutation rate p = 0 the series with increasing h is selection chemistry to biology: The mechanism of the replication pro- → irregular dynamics with two species → irregular dynam- cess is well understood in all molecular details (Biebricher ics with three species → … → extinction. At intermediate 1983; Biebricher et al. 1984, 1985) and an appropriate rep- p-values, we find quasispecies → irregular dynamics → lication assay serves for in vitro evolution studies (Mills oscillations with highly irregular spacings. At high mutation et al. 1967; Biebricher 1983). The mutation-selection sce- rates, we are dealing with quasispecies → irregular dynamics nario was found to provide an appropriate molecular basis with increasing numbers of dominant subspecies → stochas- for understanding also virus evolution [For a recent survey tic hypercycle dynamics (Fig. 9). see the contributed volume (Domingo and Schuster 2016)]. More complex systems, for example, bacteria and popula- tions of cancer cells, were found to be describable by qua- Concluding remarks sispecies theory as well (Bertels et al. 2017; Covacci and Rappuoli 1998; Napoletani et al. 2013; Brumer et al. 2006). The model presented here has been conceived with modu- In the strong mutation scenario Darwin’s view of evo- lar structure in the sense that different mechanisms can be lution has to be modified. Not a single fittest genotype is applied for each of the three basic components. Here, it has selected but a uniquely defined distribution of genotypes, been presented in its simplest form. Each of the three mod- which is represented by the largest eigenvector of a value ules, competition, cooperation and variation, can be made matrix that represents the long time or stationary solution arbitrarily complex. Variation, for example, can be extended of Eigen’s mutation-selection equation. The mean fitness of to include more elaborate mutation mechanisms and recom- populations is not always optimized since situations can be bination as well as environmental influences. Even in case constructed in which the fitness is decreasing in the approach of viruses the reproduction mechanism is commonly much towards the stationary state. A trivial but illustrative example more elaborate and comprises a whole molecular machinery of decreasing fitness during evolution considers a homoge- instead of a single enzyme. Virus reproduction may include neous population consisting exclusively of fittest genotypes also the defense system of the host, epigenetic phenomena as initial condition: Mutations introduce mutants into the could be taken into account through simultaneous consid- population and since they have lower fitness by definition the eration of several generations, and for higher organisms the mean fitness is doomed to decrease. Such situations, how - real challenge in reproduction is to deal with the enormous ever, are rather rare and Darwinian optimization still remains complexity of development in a form that is simple enough a very powerful heuristic that applies to almost all scenarios. for modeling. Cooperation at the molecular level could also For error rate parameters exceeding a critical value p , t he cr involve reproductive autocatalytic networks whereas social largest eigenvector approaches the uniform distribution over phenomena in reproductive groups or societies represent the the entire sequence space, which is the exact solution for currently highest step in the open ended complexity increase the value p =  p leading to incorporations of correct and of biological evolution. Cooperation has been frequently incorrect nucleotides with equal probabilities—for binary modeled by game theory Maynard Smith (1982); Hofbauer sequences this happens at  p = 1 −  p = . In realistic popu- and Sigmund (1998). There is no limitation to make the lations, the uniform distribution is incompatible with a dis- model more complex, the problem evidently is to include crete quasispecies (17). Instead populations are observed the desired phenomena but to keep the model sufficiently that migrate randomly through sequence space Higgs and simple for mathematical analysis or simulation. Derrida (1991); Huynen et al. (1996). In the simple form in which the model was introduced In the second half of the twentieth century, most of the here, it has been tested experimentally by in vitro evolution molecular insights into reproduction and inheritance came experiments [For an overview of early works on this subject from viruses and bacteria and a high percentage of molecu- see (Spiegelman 1971; Biebricher 1983); as a recent review lar biologists thought that the basic regulation mechanisms we mention (Joyce 2007)]. The kinetic equations describing of gene activities are understood. Eukaryotic cells, however, replication and mutation were introduced 1971 by Manfred are not “giant bacteria”. Although the genetic code is the Eigen in his scholarly written paper on self-organization of same, the gene expression and inheritance system of higher biological macromolecules (Eigen 1971). Eigen’s mutation- organisms are different from the prokaryotic one and much selection equation describes the evolution of the distribution more complex. A true wealth of information on eukaryotic 1 3 422 European Biophysics Journal (2018) 47:403–425 cells has been discovered in the past fifty years, but gene autocatalytic systems consist almost always of complex expression in animals, plants, and fungi is still a subject of reaction networks rather than single-step reactions [as current cutting-edge research. Most of the recently revealed examples for attempts to construct simple systems of this gene expression regulation mechanisms are subsumed under kind see (McCaskill 1997; Wlotzka and McCaskill 1997)]. the notion of epigenetics for which an operational definition John Maynard Smith and Eörs Szathmáry collected a true has been proposed at the Cold Spring Harbor Meeting in the wealth of evidence for the historic occurrence of such year 2008 (Berger et al. 2009): major evolutionary transitions (Maynard Smith and Sza- An epigenetic trait is a stably heritable phenotype result- thmáry 1995) in the evolution of life. ing from changes in a chromosome without alterations in It is illustrative to think about transitions in terms of the DNA sequence. thresholds: Up to a certain value of the transition determin- The diversity of epigenetic effects on gene regulation is ing parameter the typical feature is hardly detectable and enormous. It ranges from specific methylation of DNA, in does not become evident before the parameter exceeds the particular cytosine methylation in position 5 of CpG ele- transitions value. Accordingly, such thresholds correspond ments (Zemach et al. 2010), histone methylation and acety- to sharp transitions. Nevertheless, it appears useful to be less lation (Lawrence et al. 2016) to post-transcriptional RNA- fussy and to accept the notion of threshold also for smooth methylation of adenine in position 6 (Barbosa Dogini et al. transitions. On the three faces of the coordinate system 2014; Yue et al. 2015) and small interfering RNAs (He and (Fig. 5) we observe an error threshold between selection and Hanon 2004). Epigenetics provides an extremely diverse, random replication, a cooperation threshold between selec- complex and flexible richness of regulatory actions on genes, tion and symbiosis, and a mutation threshold that separates which so far was not yet cast into a comprehensive theory the regime of independent replication of all subspecies from and precisely this is one of the greatest challenges for the mutual support through frequent mutation. future of evolutionary biology. Understanding evolution implies knowledge on the rela- There is neither a convincing theoretical model nor tion between genotypes being DNA or RNA sequences and experimental evidence that Darwinian evolution leads to phenotypes giving rise to all fitness relevant parameters. The an obligatory increase in complexity. The combination metaphor of an abstract fitness landscape has been originally of competitive selection and cooperation, however, may introduced by Sewall Wright for the purpose of illustration lead from one level of complexity to the next higher one (Wright 1922). In formal mathematical terms such a relation by integration of competitors through cooperation (Sza- can be modeled as a mapping from a genotype or sequence thmáry and Maynard Smith 1995; Schuster 1996). The space into fitness values. In molecular structural biology evolution model presented here proposes a mechanism such a mapping is split into two parts: (i) a mapping from for this integration of competitors and identifies the abun- sequences into structures or genotypes into phenotypes, and dance of resources as one driving force towards higher (ii) a second mapping assigning fitness values to structures complexity. This simple model distinguishes four steps or phenotypes (Schuster 2016). Computer models of RNA (Schuster 1996): (i) Initially the systems consists of inde- evolution with sequence–secondary structure–fitness map- pendent replicators competing for a single resource, (ii) pings have been studied in the past (Fontana and Schuster the capability of cooperative interaction allows to form 1987; Fontana et al. 1989; Fontana and Schuster 1998a, b an autocatalytic network, which couples the replicators and these studies provided the basis for a definition of evo- and suppresses competition but makes the network vulner- lutionary nearness of phenotypes in the presence of neu- able to exploitation by parasites, which consume resources trality (Stadler et al. 2001). With more and more data on without contributing a share to the common properties, sequences and fitness values of mutants becoming currently (iii) the members of the autocatalytic network are sepa- available fitness landscapes may also be determined directly rated from the environment by means if a suitable bound- by experiment (Kouyos et al. 2012) and it is not risky to ary that prevents the system from exploitation and allows predict that genotype–phenotype relations will become an for the formation of a new unit at a hierarchically higher integral part of evolution models in the near future. Then level, and (iv) the individual units at the higher level diver- evolution can be described in a self-contained manner where sify by variation, compete for common resources, Darwin- the genotype–phenotype mapping is a genuine part of the ian selection sets in and takes place now a the higher level. evolving system. The previously autonomous units at the lower level lost Acknowledgements Open access funding provided by University of their autonomy at least in part when they were integrated Vienna. This contribution presents the continuation of a twenty five into the higher unit of selection. Although modeling major years lasting, wonderful cooperation with Manfred Eigen, which has transitions as shown in "Competition and cooperation" is been initiated in 1968 when the author was Post-Doc in the Eigen labo- ratory at the Max Planck-Institute for Physical Chemistry in Göttingen, not difficult, suitable experimental molecular models are Germany. The author wants to thank for the unique opportunity of very hard to conceive, because second- and higher-order 1 3 European Biophysics Journal (2018) 47:403–425 423 participating in the development of Manfred Eigen’s visionary molecu- Dobzhansky T, Ayala FJ, Stebbins GL, Valentine JW (1977) Evolution. lar view of life. Many years of cooperation and continuous discussions W. H. Freeman & Co., San Francisco with Professors Karl Sigmund, Josef Hofbauer, Walter Fontana, Peter Domingo E, Schuster P (eds) (2016) Quasispecies: from theory to F. Stadler, Ivo L. Hofacker, Christoph Flamm and others are gratefully experimental systems, current topics in microbiology and immu- acknowledged. nology, vol 392. Springer, Berlin Domingo, E., Schuster, P.: What is a quasispecies? Historical origins and current scope. In: E. Domingo, P. Schuster (eds.) Qua- Open Access This article is distributed under the terms of the Crea- sispecies: From Theory to Experimental Systems, Current tive Commons Attribution 4.0 International License (http://creat iveco Topics in Microbiology and Immunology, vol. 392, chap. 1, mmons.or g/licenses/b y/4.0/), which permits unrestricted use, distribu- pp. 1–22. Springer-Verlag, Berlin (2016) tion, and reproduction in any medium, provided you give appropriate Doob JL (1942) Topics in the theory of Markoff chains. Trans Am credit to the original author(s) and the source, provide a link to the Math Soc 52:37–64 Creative Commons license, and indicate if changes were made. Doob JL (1945) Markoff chains—Denumerable case. Trans Am Math Soc 58:455–473 Dykhuizen DE, Hartl DL (1983) Selection in chemostats. Microbiol Rev 46:150–168 References Ebeling W, Mahnke R (1979) Kinetics of molecular replication and selection. Zagadnienia Biofizyki Współczesnej 4:119–128 Albery WJ, Knowles JR (1976) Evolution of enzyme function and the Eigen M (1971) Selforganization of matter and the evolution of bio- development of catalytic efficiency. Biochemistry 15:5631–5640 logical macromolecules. Naturwissenschaften 58:465–523 Albery WJ, Knowles JR (1977) Efficiency and evolution of enzyme Eigen M, McCaskill J, Schuster P (1989) The molecular quasispe- catalysis. Angew Chem Internat Ed Engl 16:285–293 cies. Adv Chem Phys 75:149–263 Baake E, Gabriel W (1999) Biological evolution through mutation, Eigen M, Schuster P (1977) The hypercycle. A principle of natural selection, and drift: an introductory review. In: Stauffer D (ed) self-organization. Part A: emergence of the hypercycle. Natur- Annual review of computational physics VII. World Scientific, wissenschaften 64:541–565 Singapore, pp 203–264 Eigen M, Schuster P (1978) The hypercycle. A principle of natural Berger SL, Kouzarides T, Shiekhattar R, Shilatifard A (2009) An oper- self-organization. Part B: the abstract hypercycle. Naturwis- ational definition of epigenetics. Genes Dev 23:781–783 senschaften 65:7–41 Bertels F, Gokhale CS, Traulsen A (2017) Discovering complete qua- Eigen M, Schuster P (1978) The hypercycle. A principle of natural sispecies in bacterial genomes. Genetics 206:2149–2157 self-organization. Part C: the realistic hypercycle. Naturwis- Biebricher CK (1983) Darwinian selection of self-replicating RNA senschaften 65:341–369 molecules. In: Hecht MK, Wallace B, Prance GT (eds) Evolu- Ewens WJ, Lessard S (2015) On the interpretation and relevance tionary biology, vol 16. Plenum Publishing Corporation, New of the fundamental theorem of natural selection. Theor Popul York, pp 1–52 Biol 104:59–67 Biebricher CK, Eigen M, Gardiner WC Jr (1983) Kinetics of RNA Feller W (1940) On the integro-differential equations of purely replication. Biochemistry 22:2544–2559 discontinuous Markoff processes. Trans Am Math Soc Biebricher CK, Eigen M, William C, Gardiner J (1984) Kinetics of 48:488–515 RNA replication: plus-minus asymmetry and double-strand for- Fisher RA (1930) The genetical theory of natural selection. Oxford mation. Biochemistry 23:3186–3194 University Press, Oxford Biebricher CK, Eigen M, William C, Gardiner J (1985) Kinetics of Fisher RA (1958) The genetical theory of natural selection, 2nd edn. RNA replication: competition and selection among self-repli- Dover Publications Inc, New York cating RNA species. Biochemistry 24:6550–6560 Fontana W, Schnabl W, Schuster P (1989) Physical aspects of evolu- Brumer Y, Michor F, Shakhnovich EI (2006) Genetic instability and tionary optimization and adaptation. Phys Rev A 40:3301–3321 the quasispecies model. J Theor Biol 241:216–222 Fontana W, Schuster P (1987) A computer model of evolutionary opti- Bryson V, Szybalski W (1952) Microbial selection. Science 116:45–46 mization. Biophys Chem 26:123–147 Covacci A, Rappuoli R (1998) Helicobacter pylori: molecular evolu- Fontana W, Schuster P (1998a) Continuity in evolution. On the nature tion of a bacterial quasi-species. Curr Opt Microbiol 1:96–102 of transitions. Science 280:1451–1455 Crow JF, Kimura M (1970) An introduction to population genetics Fontana W, Schuster P (1998b) Shaping space. The possible and the theory. Sinauer Associates, Sunderland attainable in RNA genotype-phenotype mapping. J Theor Biol Dogini BD, Pascoal VD, Avansini SH, Vieira AS, Campos TC, Lopes- 194:491–515 Cendes I (2014) The new world of RNAs. Genet Mol Biol 37(1 Gillespie DT (1976) A general method for numerically simulating the suppl):285–293 stochastic time evolution of coupled chemical reactions. J Comp de Pasquale F, Tartaglia P, Tombesi P (1980) Stochastic approach to Phys 22:403–434 chemical instabilities. Anomalous fluctuations transient behavior. Gillespie DT (1977) Exact stochastic simulation of coupled chemical Lettere al Nuovo Cimento 28:141–145 reactions. J Phys Chem 81:2340–2361 de Pasquale F, Tartaglia P, Tombesi P (1982) New expansion technique Gillespie DT (2007) Stochastic simulation of chemical kinetics. Annu for the decay of an unstable state. Phys Rev A 25:466–471 Rev Phys Chem 58:35–55 de Pasquale F, Tartaglia P, Tombesi, P (1982) Transient-behavior near Gillespie JH (1983) Some properties of finite populations experienc- an instability point. Vectorial stochastic representation of the ing strong selection and weak mutation. Am Nat 121:691–708 Malthus. Nuovo Cimento B 69:228–238 Goel NS, Richter-Dyn N (1974) Stochastic models in biology. Aca- Dean AM, Thronton JW (2007) Mechanistgic approaches to the study demic Press, New York of evolution: the functional synthesis. Nat Rev Genet 8:675–688 Gorini L, Beckwith JR (1966) Suppression. Annu Rev Microbiol Deuflhard P, Huisinga W, Jahnke T, Wulkow M (2008) Adaptive dis- 20:401–422 crete Galerkin methods applied to the chemical master equation. Hamming RW (1950) Error detecting and error correcting codes. Bell SIAM J Sci Comput 30:2990–3011 Syst Tech J 29:147–160 1 3 424 European Biophysics Journal (2018) 47:403–425 Hamming RW (1986) Coding and information theory, 2nd edn. Pren- Nicolis G, Prigogine I (1977) Self-organization in nonequilibrium sys- tice-Hall, Englewood Cliffs tems. Wiley, New York Hartman PE, Roth JR (1973) Mechanisms of suppression. Adv Genet Novick A, Szillard L (1950) Description of the chemostat. Science 17:1–105 112:715–716 He L, Hanon GJ (2004) Micro RNAs: small RNAs with a big role in Novick A, Szillard L (1950) Experiments with the chemostat on gene regulation. Nat Rev Genet 5:522–531 spontaneous mutations of bacteria. Proc Natl Acad Sci USA Higgs PG, Derrida B (1991) Stochastic models for species formation in 36:708–719 evolving populations. J Phys A Math Gen 24:L985–L991 Noyes RM, Field RJ, Körös E (1972) Oscillations in chemical sys- Hofbauer J, Mallet-Paret J, Smith HL (1991) Stable periodic solutions tems. I. Detailed mechanism in a system showing temporal for the hypercycle system. J Dyn Diff Equ 3:423–436 oscillations. J Am Chem Soc 94:1394–1395 Hofbauer J, Sigmund K (1998) Evolutionary games and population Okasha S (2008) Fisher’s fundamental theorem of natural selec- dynamics. Cambridge University Press, Cambridge tion—a philosophical analysis. Br J Philos Sci 59:319–351 Huang CI, Tu MF, Lin HH, Chen CC (2017) Variation approach Park JM, Muñoz E, Deem MW (2010) Quasispecies theory for finite to error threshold in generic fitness landscape. Chin J Phys populations. Phys Rev E 81:e011,902 55:606–618 Phillipson PE, Schuster P (2009) Modeling by Nonlinear Differen- Husimi Y (1989) Selection and evolution of a bacteriophages in cell- tial Equations. Dissipative and Conservative Processes, World stat. Adv Biophys 25:1–43 Scientific Series on Nonlinear Science A, 69th edn. World Sci- Husimi Y, Nishigaki K, Kinoshita Y, Tanaka T (1982) Cellstat—a con- entific, Singapore tinuous culture system of a bacteriophage for the study of the Phillipson PE, Schuster P, Kemler F (1984) Dynamical machinery of mutation rate and the selection process at the DNA level. Rev a biochemical clock. Bull Math Biol 46:339–355 Sci Instr 53:517–522 Plutynski A (2006) What was Fisher’ fundamental theorem of natural Huynen MA, Stadler PF, Fontana W (1996) Smoothness within rug- selection and what was it for? Stud Hist Phil Biol Biomed Sci gedness. The role of neutrality in adaptation. Proc Natl Acad Sci 37:50–82 USA 93:397–401 Prelich G (1999) Suppression mechanisms. themes from variations. Jahnke T, Huisinga W (2007) Solving the chemical master equation Trends Genet 15:261–266 for monomolecular reaction systems analytically. J Math Biol Price GR (1972) Fisher’s ‘fundamental theorem’ made clear. Ann 54:1–26 Hum Genet 36:129–140 Jones BL, Enns RH, Rangnekar SS (1976) On the theory of selection Sagués F, Epstein IR (2003) Nonlinear chemical dynamics. J Chem of coupled macromolecular systems. Bull Math Biol 38:15–28 Soc Dalton Trans 2003:1201–1217 Jones BL, Leung HK (1981) Stochastic analysis of a non-linear model Schmidt LD (2005) The engineering of chemical reactions, 2nd edn. for selection of biological macromolecules. Bull Math Biol Oxford University Press, New York 43:665–680 Schnabl W, Stadler PF, Forst C, Schuster P (1991) Full characteri- Joyce GF (2007) Forty years of in vitro evolution. Angew Chem Inter- zatin of a strang attractor. Chaotic dynamics on low-dimen- nat Ed 46:6420–6436 sional replicator systems. Phys D 48:65–90 Joyce P, Rokyta DR, Beisel CL, Orr HA (2008) A general extreme Schuster P (1996) How does complexity arise in evolution. Nature’s value theory model for the adaptation of DNA sequences under recipe for mastering scarcity, abundance, and unpredictability. strong selectiion and weak mutation. Genetics 180:1627–1643 Complexity 2/1(1):22–30 Kimura M (1955) Solution of a process of random genetic drift with a Schuster P (2016) Increase in complexity and information through continuous model. Proc Natl Acad Sci USA 41:144–150 molecular evolution. Entropy 18:e397 Kimura M (1955) Stochastic processes and distribution of gene fre- Schuster P (2016) Major transitions in evolution and in technology. quencies under natural selection. Cold Spring Harb Symp Quant What they have in common and where they differ. Complexity Biol 20:33–53 21/4(4):7–13 Kimura M (1983) The neutral theory of molecular evolution. Cam- Schuster P (2016) Quasispecies on fitness landscapes. In: bridge University Press, Cambridge E. Domingo, P. Schuster (eds) Quasispecies: from theory to Kolmogorov AN (1931) Über die analytischen Methoden in der Wahr- experimental systems, Current Topics in Microbiology and scheinlichkeitsrechnung. Math Ann 104:415–458 In German Immunology, vol. 392, chap. 4. Springer-Verlag, Berlin, pp Koltermann A, Kettling U (1997) Principles and methods of evolution- 61–120. ary biotechnology. Biophys Chem 66:159–177 Schuster P (2016) Some mechanistic requirements for major transi- Kouyos RD, Leventhal GE, Hinkley T, Haddad M, Whitcomb JM, tions. Phil Trans Roy Soc Lond B 371:e20150,439 Petropoulos CJ, Bonhoeffer S (2012) Exploring the complexity Schuster P (2016) Stochasticity in processes. Fundamentals and appli- of the HIV-1 fitness landscape. PLos Genet 8:e1002,551 cations in chemistry and biology. Springer series in synergetics. Lawrence M, Daujat S, Schneider R (2016) Lateral thinkiing: how Springer, Berlin histone modifications regulate gene expression. Trends Genet Schuster P (2017) A mathematical model of evolution. MATCH Com- 32:42–56 munications in Mathematical and in Computer Chemistry 78, Maynard Smith J (1982) Evolution and the theory of games. Cambridge submitted University Press, Cambridge Schuster P (2017) Molecular evolution between chemistry and biol- Maynard Smith J, Szathmáry E (1995) The major transitions in evolu- ogy. The interplay of competition, cooperation, and mutation. tion. W. H. Freeman-Spektrum, Oxford European Biophysics Journal 46, submitted McCaskill JS (1997) Spatially resolved in vitro molecular ecology. Schuster P, Sigmund K (1985) Dynamics of evolutionary optimization. Biophys Chem 66:145–158 Ber Bunsenges Phys Chem 89:668–682 Mills DR, Peterson RL, Spiegelman S (1967) An extracellular Darwin- Schuster P, Sigmund K, Wolff R (1978) Dynamical systems under ian experiment with a self-duplicating nucleic acid molecule. constant organization I. Topological analysis of a family of non- Proc Natl Acad Sci USA 58:217–224 linear differential equations - A model for catalytic hypercycles. Napoletani D, Signore M (2013) Cancer quasispecies and stem-like Bull Math Biol 40:734–769 adaptive aneuploidy. F1000Research 2:e268 1 3 European Biophysics Journal (2018) 47:403–425 425 Schuster P, Sigmund K, Wolff R (1979) Dynamical systems under con- Tarazona P (1992) Error thresholds for molecular quasispecies as phase stant organization III. Cooperative and competitive behavior of transitions: From simple landscapes to spin glasses. Phys Rev A hypercycles. J Diff Equ 32:357–368 45:6038–6050 Schuster P, Sigmund K, Wolff R (1980) Dynamical systems under con- Thompson CJ, McBride JL (1974) On Eigen’s theory of the self-organ- stant organization II. Homogenoeus growth functions of degree ization of matter and the evolution of biological macromolecules. p = 2 . SIAM J Appl Math 38:282–304 Math Biosci 21:127–142 Sniegowski PD, Gerrish PJ (2010) Beneficial mutations and the dynam- Wiehe T (1997) Model dependency of error thresholds: The role of ics of adaptation in asexual populations. Phil Trans R Soc B fitness functions and contrasts between the finite and infinite sites 356:1255–1263 models. Genet Res Camb 69:127–136 Spiegelman S (1971) An approach to the experimental analysis of pre- Wlotzka B, McCaskill JS (1997) A molecular predator and its prey: cellular evolution. Quart Rev Biophys 4:213–253 coupled isothermal amplification of nucleic acids. Chem Biol Stadler BRM, Stadler PF, Wagner GP, Fontana W (2001) The topology 4:25–33 of the possible: Formal spaces underlying patterns of evolution- Wright S (1922) Coefficients of inbreeding and relationship. Am Nat ary change. J Theor Biol 213:241–274 56:330–338 Strunk G, Ederhof T (1997) Machines for automated evolution experi- Yue Y, Liu J, He C (2015) RNA N -methyladenosine methylation ments in  vitro based on the serial-transfer concept. Biophys in post-ttranscriptional gene expression regulation. Genes Dev Chem 66:193–202 29:1343–1355 Swetina J, Schuster P (1982) Self-replication with errors - A model for Zemach A, McDaniel IE, Silva P, Zilberman D (2010) Genome-wide polynucleotide replication. Biophys Chem 16:329–345 evolutionary analysis of eukaryotic DNA methylation. Science Szathmáry E, Maynard Smith J (1995) The major evolutionary transi- 328:916–919 tions. Nature 374:227–232 1 3

Journal

European Biophysics JournalSpringer Journals

Published: Mar 2, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off