## Abstract

Major transitions in nature and human society are accompanied by a substantial change towards higher complexity in the core of the evolving system. New features are established, novel hierarchies emerge, new regulatory mechanisms are required and so on. An obvious way to achieve higher complexity is integration of autonomous elements into new organized systems whereby the previously independent units give up their autonomy at least in part. In this contribution, we reconsider the more than 40 years old hypercycle model and analyse it by the tools of stochastic chemical kinetics. An open system is implemented in the form of a flow reactor. The formation of new dynamically organized units through integration of competitors is identified with transcritical bifurcations. In the stochastic model, the fully organized state is quasi-stationary whereas the unorganized state corresponds to a population with natural selection. The stability of the organized state depends strongly on the number of individual subspecies, *n*, that have to be integrated: two and three classes of individuals, and , readily form quasi-stationary states. The four-membered deterministic dynamical system, , is stable but in the stochastic approach self-enhancing fluctuations drive it into extinction. In systems with five and more classes of individuals, , the state of cooperation is unstable and the solutions of the deterministic ODEs exhibit large amplitude oscillations. In the stochastic system self-enhancing fluctuations lead to extinction as observed with . Interestingly, cooperative systems in nature are commonly two-membered as shown by numerous examples of binary symbiosis. A few cases of symbiosis of three partners, called *three-way symbiosis*, have been found and were analysed within the past decade. Four-way symbiosis is rather rare but was reported to occur in fungus-growing ants. The model reported here can be used to illustrate the interplay between competition and cooperation whereby we obtain a hint on the role that resources play in major transitions. Abundance of resources seems to be an indispensable prerequisite of radical innovation that apparently needs substantial investments. Economists often claim that scarcity is driving innovation. Our model sheds some light on this apparent contradiction. In a nutshell, the answer is: scarcity drives optimization and increase in efficiency but abundance is required for radical novelty and the development of new features.

This article is part of the themed issue ‘The major synthetic evolutionary transitions’.

## 1. Properties of major transitions

Biological evolution is not a smooth process and changes are almost always abrupt. The large steps in this development have been called *major transitions* by John Maynard Smith & Eörs Szathmáry [1]. Major transitions in nature and human society represent the formation of new, higher *units of selection* and accordingly require a substantial change towards higher complexity in the core of the evolving system. New features are established, novel hierarchies emerge, new regulatory mechanisms are required and so on. In the transition, the system makes use of functions that are already available but in a different context and builds the novel organization upon them. The development of new capabilities requires sufficient resources even on the molecular level [2]. The high costs are, for example, clearly seen in the case of the recruitment of new enzymatic function through gene duplication [3,4]: a full duplication of the entire yeast genome resulted in only two additional genes and proteins. Two typical examples of major transitions are mentioned here: (i) the transition from prokaryotic to eukaryotic life and (ii) the industrial revolution. In both cases, cheap energy from new sources became available, in particular, the presence of a sufficiently high concentration of free oxygen in the atmosphere allowed for the development of oxidative phosphorylation, which increases the energy output by ATP production per molecule of glucose from two in fermentation to 30 in mitochondrial oxidative phosphorylation, that is by a factor of 15. In order to make oxidative phosphorylation possible, however, a large number of *molecular investments* leading to the complex eukaryotic cell were necessary [5,6]. Formally, the industrial revolution [7] is characterized by a similar situation: an enormous resource of energy has been available in the form of fossil fuels, especially coal, but the exploitation of this resource on an industrial scale required huge investments. Abundance of cheap resources appears to be an indispensable prerequisite of radical innovations because new functions or novel technologies need investments. Economists often raise the issue encapsulated in the familiar quotation *‘necessity is the mother of invention’* and claim that scarcity is driving innovation. Natural selection when understood as an optimization process is not bound to major investments and can be based indeed on a multitude of small improvements that lead to an increase in fitness. A related question that can be asked from the point of view of population dynamics concerns the dominant terms in the kinetic equations: *‘Are competition or cooperation arising with increased resources allowing for higher population density?’*

The most natural way to achieve higher complexity is to integrate competitors into a new functional unit. What is needed, however, is a mechanism that allows for suppression of competition. In other words, we have to look for a dynamical interaction that results in cooperation of competitors. Other than in economics, multiplication in the sense of copying is a *conditio sine qua non* for survival in biology. The notion of *replicator* was created for individuals in species that have to reproduce in order to avoid going extinct [8]. Hence, all biological agents are replicators and replication implies autocatalysts in the language of chemical kinetics. About 40 years ago a simple class of mechanisms, which provides dynamical coupling of replicating units, was proposed and the simplest chemical reaction networks of this class were called *hypercycles* [9–11]. The principle of hypercyclic coupling is straightforward: competitors are forced to cooperate, because their reproductive success is bound to the presence of agents from another class. A stable functional unit is formed, for example, through closing a cycle of mutual dependence, which yields the hypercycle. Needless to say, the transition from competition to cooperation is only the first step of major transitions as we shall explain in §6.

Empirical evidence for dynamical coupling of reproducing species is abundant for the smallest possible systems consisting of two cooperating species in the form of symbiosis. The presumably most common form is endosymbiosis [6,12] in the eukaryotic cells of animals and fungi where the cellular nucleus and the mitochondria reproduce autonomously but strong mutual dependence is caused by the majority of mitochondrial genes being stored in the nuclear genome and strong metabolic interaction because oxidative phosphorylation is performed only in mitochondria. The extension to three cooperating partners has happened in the cells of plants and algae where the chloroplasts represent a second class of endosymbionts [13,14]. Several other examples of three-way symbiosis are known, for example [15,16] and the systematic studies on ants–fungi–bacteria systems [17,18]. Examples of four-way symbiosis seem to be rare [19, p. 71].

Mathematical analysis of dynamical systems derived from hypercycles and more general replicator equations [8] has been performed and was reported 40 years ago [20–22]. Few attempts of modelling cooperative dynamics as a stochastic process were made (examples are [23–28]) and no systematic study on the interplay of competition and cooperation is known. The computational facilities for studying dynamical systems and simulating stochastic processes have increased enormously since the 1980s and therefore it appears worth revisiting the dynamical models of cooperation in order to make predictions based upon assumptions which are now somewhat closer to reality. In this contribution, we shall address two open questions: (i) is it possible to conceive a sufficiently simple model that shows the transition from competition and selection to cooperation, and (ii) which role is played by stochasticity in cooperative dynamical system based on replicators? Being able to find such a model and to answer question (ii), we hope to find at the same time an explanation why multipartite symbiosis becomes unstable with increasing numbers of players.

## 2. Deterministic model for transitions

In order to be able to do rigorous mathematical analysis, we have to conceive a model for transitions from competition to cooperation. Because of the enormous complexity of real world biology, we have to look for a system that fulfils the rule of Albert Einstein's version of Occam's razor: *‘Everything should be made as simple as possible, but not simpler’*. We assume a population of *n* classes of replicators or *subspecies*, . The subspecies are reproducing in an open system consisting of a flow reactor, which is fed from a large reservoir filled with stock solution containing the material A in concentration where A stands for the whole set of materials that are required for reproduction. The mechanisms of reproduction and catalysed reproduction in the flow reactor are encapsulated in chemical reaction equations
2.1*a*
2.1*b*
2.1*c*
2.1*d*
2.1*e*

Reaction (2.1*a*) supplies the material required for reproduction. The solution with A at concentration flows into a continuously stirred tank reactor (CSTR) with a flow rate *r* [29, p. 87ff.]. The reactor operates at constant volume and this implies that the same volume of solution per unit time [*t*] flows into the reactor and out from it. The outflow is described by the equations (2.1*d*) and (2.1*e*) and concerns all molecular species, A and in the same manner. Reaction (2.1*b*) describes the process of correct copying of replicators and reaction (2.1*c*) finally introduces cooperation between otherwise competing types: supports the process of copying . On the molecular level we might consider as a catalyst for the replication of .

The choice of a flow reactor as an appropriate open system is not accidental. The two parameters and *r* are easily controlled in real model experiments and they provide a tool for characterizing the open system: the resource in the form of is a measure for the amount of material available for reproduction and thus it is related to the carrying capacity of an ecosystem; the flow rate *r* provides a time constraint as is the mean residence time of a volume element in the reactor and subspecies which need more time for replication and cannot survive in the reactor. There are two limits of interest: (i) leading to the closed system and (ii) , the condition of no reaction because the mean residence time in the reactor goes to zero.

The kinetic differential equations according to mechanisms (2.1*a–e*) are formulated in continuous variables and
2.2*a*and
2.2*b*The rate parameter representing the analogue to the Althusian fitness value of and is the rate parameter for replication of catalysed by . We shall refer to the *f*_{i} values as fitness parameters and to the *k*_{ij} values as cooperation parameters despite the fact that certain combinations of *k*_{ij} values may lead to strong selection, for example the multidimensional Schlögl model [30,31].

The reaction equation (2.1*c*) takes into account catalysis for all binary combinations of replicators in the system. The occurrence of efficient catalysis, however, is a rare event and accordingly the number of catalytic interactions in mechanism (2.1*a–e*) is unrealistically large. In the hypercycle model, catalysis is restricted to a single cooperation partner. The system shows special properties when mutual catalysis forms a closed loop
2.1*c*where we write instead of as by definition. In other words, catalysis is confined to a one-dimensional cyclic array
2.3and this simplifies the ordinary differential equations (ODEs) to
2.4*a*
2.4*b*The ODEs (2.4*a,b*) combine competition and cooperation kinetics and are candidates for modelling transitions from Darwinian selection to stable symbiosis. At the same time, it seems that nothing in the model can be omitted without losing the capability to describe the transition and therefore we conjecture that it is also the simplest model possible.

Steady state analysis reveals the long-time behaviour of the system. Properly, we start from equation (2.2*b*) and find
2.4*c*Accordingly, we have two solutions for the stationary concentration of every subspecies: (i) extinction of given by , and (ii) survival of expressed by , requiring . Provided all combinations of solutions are compatible, we can take either (i) or (ii) for each subspecies and this yields stationary states in total. Three states are of particular importance here: the state of extinction, , the state of Darwinian selection of the fittest subspecies , , and the state of cooperation, .^{1} For the two smallest systems with and , we shall analyse all states, but for arbitrary *n* only extinction, selection and cooperation will be considered.

From equation (2.4*c*) follow expressions for the stationary values of the variables () which are functions of . For the calculation of the stationary values we make use of the total concentration , which fulfils
2.5

At the state of cooperation, , no subspecies vanishes and the stationary concentrations fulfil . For arbitrary *n*, the concentration of A is obtained from the quadratic equation
2.6The root
2.7belongs to . The second root belongs to an unstable state , which lies often at negative *x*_{i} values. It will be mentioned in the case of pure hypercycles where implying .

In summary, the three steady states of particular interest are

(i) , the state of extinction: for all and .

(ii) , the state of selection: , for , and with , the following constraint on the flow rate , and for we have extinction.

(iii) , the state of cooperation: , , where is one of the two solutions of the quadratic equation (2.6) for flow rates in the range ; the state of cooperation does not exist for and the system converges to the state of extinction.

The population converges to one of the three states, , and , or to one of the intermediate stable states with depending uniquely on the values of rate parameters and initial conditions. If is not present, any other subspecies with highest fitness among the replicators actually present can be selected. It is worth noting that the rate parameters for cooperation, , do not show up in the expression for the selection states, whereas the fitness values have an influence on the position of the cooperation equilibrium.

In figure 1, we illustrate the competition–cooperation system by means of two numerical examples. We choose and tune the flow reactor such that the concentration of the subspecies is gradually increasing from small initial values, , to the stationary concentration. Apart from the molecules for the start of the reaction the reactor is empty, , at first the concentration increases rapidly because the influx of A is a *zeroth* order process and dominates at low concentrations. As the concentrations of the subspecies begin to grow, A is consumed by the reaction (2.1*b*) and the concentration of A goes down to the stationary value . A different behaviour is observed for low and high population sizes: selection dominates at small population size and as the population becomes larger cooperative behaviour takes over and the system eventually settles down in the state of cooperation, . This interesting and appealing feature is also a result of the molecularity of the reaction. The reproduction reaction in the selection scenario is bimolecular, , whereas it is trimolecular in the cooperation case, , and for this reason low concentrations in the sense of scarce resources favour selection, but abundant resources favour cooperation. The trimolecular terms give rise to a strong increase of the reproduction rates, , with population density, a phenomenon that is called the *Allee* effect in macroscopic biology [32].

Two kinds of series are illustrative: (i) increasing flow rate *r* at constant resources , and (ii) increasing at constant *r*. In order to keep the system as simple as possible only two and three subspecies, , were chosen. We shall see later on that these two cases represent the most important ones. In case of we are dealing with four states, , , and . Out of the two selection states only one is asymptotically stable, the one with the higher fitness value . At sufficiently small flow rates the cooperative state is asymptotically stable (figure 2). It loses stability through a transcritical bifurcation and it is straightforward to calculate the critical flow rate by setting equal the two stationary concentrations of A
which leads to
2.8*a*

Further increase in the flow rate *r* crosses the range of selection of the fittest species , which ends at a second transcritical bifurcation above which the state of extinction, is stable:
2.8*b*

In figure 2*a*, the two consecutive bifurcations are shown for a numerical example with . The second scenario, increasing resources at constant flow rate *r* (figure 2*b*) is relevant for the question we are raising here. Very low resources lead to extinction, implying that is globally, asymptotically stable. Gradual increase in leads to the transcritical bifurcation (2.8*b*) and further increase of results in the cooperative state via a second transcritical bifurcation (2.8*a*). Accordingly, cooperation requires more resources than selection.

The bifurcation pattern in the system with three subspecies, , and , and stationary states are very similar. Without losing generality we order the fitness values, with . We are dealing with a state of extinction, , a cooperative state , and three states of selection , and . In addition, there are three steady states where two of the three subspecies are present. As a single subspecies is excluded we have chosen the name *state of exclusion* and use the notation , and where the absent subspecies is given in the superscript. Again, we can derive analytical formulae for the three bifurcation points by equating concentration values
2.9*a*
2.9*b*
2.9*c*

The sequence of states with increasing -values is extinction, , selection of the fittest, , exclusion of the least fit, and cooperation, . The sequence of bifurcations is shown schematically in figure 3.

The deterministic approach to the description of transitions in biology suffers from a general problem. Mutants, being the ultimate source of change in properties and functions, always start off from a single copy and stochasticity plays an important role in the early phases of biological innovation. In §§3–5 we shall discuss further aspects of competition, cooperation and their interplay, which become apparent only beyond the deterministic approach and consider the stochastic phenomena at small and very small particle numbers that are relevant for transitions.

## 3. Stochastic competition

The chemical master equation is the appropriate tool for studying kinetic processes at small particle numbers. In order to avoid confusion, we shall write the discrete random variables with upper case letters: , or . It is not difficult to write down a multivariate master equation for A and the *n* subspecies , but analytical solutions are very hard to derive particularly when the system has more than one independent variable. Numerical studies, on the other hand, are quite efficient for small particle numbers. We make use of them here in order to work out the difference between the probabilistic systems and the deterministic approach. Simulation algorithms are available which yield results that were shown to converge towards the solutions of the corresponding master equations in the limit of infinite numbers of stochastic trajectories [33–35]. For the simulations reported here, we used a slightly modified version of the *xSSA* stochastic simulation algorithm in the implementation for *Mathematica* within the *xCellerator* project [36]. Bundles of stochastic trajectories were evaluated by two different strategies: (i) final state counting in the case of systems with more than one attractor whereby the final state is defined as the family of values within a band of thermal fluctuations around the stationary expectation value and (ii) statistics of trajectories calculated in the form of sample mean values and standard deviations.

In the deterministic reference competition between subspecies with different fitness values, is described by equations (2.2*a,b*) with for all . The case of selective neutrality can be modelled by Motoo Kimura's theory of neutral evolution [37] and is not considered here. The deterministic long-time solution is unique and, depending on parameters, we observe either selection of the fittest subspecies for low flow rates, or extinction for high flow rates . It is straightforward to see that only one subspecies can be present in the long run: let us assume that two species, and , are present at the stationary state, then has to fulfil the two equations and simultaneously, which is impossible unless and this is the case which has been excluded.

In the stochastic system, we have long-time states: one state is , the state of extinction with and for all , and the other states are the *n* selection states with at which A and one subspecies, , are present at the values and . The state of extinction is an absorbing boundary and it is the only absorbing boundary of the dynamical system. Accordingly, all stochastic trajectories have to end in in the limit . Precisely speaking, is absorbing only for with because A is replenished by the influx should it have gone to extinction by accident. All selection states are *quasi-stationary*. Quasi-stationarity for a state means that it is indistinguishable in its properties from a true stationary state with the exception that it cannot last to infinity.

For very small initial particle numbers, , the selection process has a strong random component. Some numerical values obtained by simulation of the smallest possible system with are shown in table 1. If the process starts with only a single copy of a subspecies the probability that this subspecies is lost before it had a chance to reproduce is not negligible. As no mutation is taken into account, a subspecies cannot come back once it has been extinguished. In reality the chance that one particular genotype is reintroduced into the system by a revertant mutation after it had been lost is extremely small and zero for most practical purposes. Some special cases with viroids and viruses may be counterexamples ([38], p.49). Interestingly, very few copies are already sufficient to make the deterministic result the only one with appreciable probability. As the data in table 1 show, is already sufficient to achieve a probability of more than 90% for *selection of the fittest*. Figure 4 shows two characteristic trajectories leading to extinction and selection, respectively.

Mathematical analysis [39, pp. 72–75] and the computer simulations fully confirm the competitive exclusion principle or Gause's law [40–42], which states that a single resource can only support a single species. More ecological niches require specialization. In other words, provided the resources are sufficient, , the range of independent reproduction is characterized by survival of the fittest. The ultimate difference between the deterministic and the stochastic system is to be seen in the fact that in the latter a subspecies can die out in the early phase of the process when particle numbers are small. If the subspecies with highest fitness dies out first another subspecies, the one with highest but one fitness, will be selected. The expectation values in the early phase of the selection process are very close to the deterministic results but then the stochastic trajectories are driven by random events into one of several basins of attraction of the steady states , or any other one of the quasi-stationary states , whereas the deterministic system approaches the unique attractor in the basin where the initial values were lying. In this context, we mention that the competitive exclusion principle does not apply to communities of species that exhibit oscillations or deterministic chaos where dozens of species can coexist on a handful of resources as in the case of phytoplankton [43].

## 4. Stochastic cooperation

Cooperative networks of replicators with hypercyclic organization in the flow reactor were considered previously in some detail [31]. The mechanism is given by equations (2.1*a–e*) where according to (2.1*c*) we use for and otherwise. The corresponding ODE's are presented in (2.4*a,b*) with for all . It is straightforward to calculate the stationary concentrations from the equations resulting from evaluating the ODEs to zero: and with . According to the catalytic sequence (2.3) all replicators go extinct sequentially after one has become extinguished as seen from: if , and vanishes followed by , , and so on. Accordingly we are dealing with only two long-time states, and . The stationary concentration of A is in and in it is a root of the quadratic equation
4.1*a*which sustains two solutions and in the range and no real solution exists for . The smaller of the two values, , belongs to the state of cooperation . There is a third stationary state, , which we conjecture as being always unstable. It corresponds to a saddle point and it separates the basins of attraction of the two states and : depending on initial conditions the deterministic solution converges either to the state of extinction, , or the state of cooperation, . The transition from extinction to cooperation is mediated by a saddle-node bifurcation (figure 5). It is worth considering the stationary concentrations of subspecies in detail. The value implies that the concentration of the subspecies is inversely proportional to the catalytic rate parameter of its precursor in the cycle
or in other words the *catalytic flow* is constant along the cycle. The phase space topology of the hypercycle in the flow reactor analysed above is quite general. The qualitative saddle-node bifurcation pattern involving the states , and does not depend on the dimension *n* of the cycle and is valid for a variety of other embeddings in open systems [28,44,45].

The dynamical behaviour of the system is reflected completely by the symmetric case, , and we shall analyse it here because symmetry simplifies the expressions for the properties of the central stationary point, . The symmetrization can be done in full generality by means of a transformation of the *x* variables: , . Then, the position of the steady state is in the centre of the simplex and has the coordinates with
4.1*b*Stability analysis is performed by means of the Jacobian matrix, which has a simple structure with only *n* entries apart from the flow rate *r*
or where is a circulant matrix containing only zeros and one *r* in every row. Calculation of the eigenvalues of the Jacobian yields [31]
4.1*c*

Accordingly, is negative for and positive for in the whole range where the two states exist and is always unstable as conjectured above. Apart from the factor *r* the eigenvalues are the complex roots of unity except the root . Hence, all real parts of the eigenvalues are negative for and is asymptotically stable; for the state is unstable because there are roots with positive real parts, and it is marginally stable for because we have the complex conjugate pair . As a matter of fact, the linear term in the stability analysis for with vanishes but there is a negative higher-order term that causes slow oscillatory convergence to the stationary state [20]. For , trajectories of hypercycles were computed by numerical integration and show *relaxation oscillations* becoming harder with increasing *n* [46]. The existence of a stable limit cycle for systems with was proven many years later in an elegant piece of analytical work [47]. Bifurcation analysis of oscillating hypercycles () has been reported recently [48]. A series of studies on small hypercycles revealed interesting dynamical phenomena such as *delayed transitions* and *ghosts* [49,50]. *Survival* of spatially extended hypercycles with follows similar rules and has been studied by means of a cellular automaton in [44].

The stochastic cooperative system in the flow reactor exhibits one absorbing boundary, the state of extinction , a quasi-stationary central cooperative state . The unstable state cannot be recognized unambiguously in the stochastic approach because of the superposition of fluctuations. In the deterministic system, the initial conditions determine uniquely the particular trajectory. The outcome of a stochastic trajectory, however, is dependent on all random events in the early phase of the process, and we can only give probabilities to end up either in state or in state , which are determined by the initial conditions, , for , as well as the sequence of random events (table 2). A detailed numerical study was performed for and the initial values , . The individual trajectories (an example is shown in figure 6) show either an approach to the state of extinction with fluctuating around the stationary value or a transient ending in a very sharp transition to the cooperative state with all three random variables fluctuating around the deterministic stationary values: , and . For constant and , and a given sequence of random events as controlled by the seeds for the pseudorandom number generator, the outcome of the stochastic process, or , depends critically on : for , the system converges to the cooperative state whereas extinction is observed for . The critical value is extremely variable, for and different seeds corresponding to different sequences of random events we recorded values from to . Apart from the probability to go extinct the initial conditions have remarkably little influence on the shape of the trajectories after the transient (figure 6).

An important issue is the dependence of the cooperative system on the number of subspecies *n*. For and 3 the trajectories converge towards the inner stationary state and then fluctuate around it provided the system did not go extinct in the initial phase. For the system starts to oscillate around , the oscillations increase in amplitude until one subspecies dies out, then extinction of subspecies proceeds along the chain (2.3), and the state of extinction is reached. The case is intermediate with a weakly stable central stationary state in the deterministic equations. For small systems (e.g. ), fluctuations destabilize the four-dimensional system, it is only a matter of time until a fluctuation wipes out one subspecies and then the whole system goes extinct (figure 7). In the case of a larger system (e.g. ), the quasi-stationary state lives longer and as expected approaches the stability behaviour of the deterministic system in the limit . In nature, cooperative systems are formed by a mutational event producing a single copy and accordingly the small size scenario is applicable. Eventually, we are left with the result that the kinetic model predicts that simple cooperative units of the hypercycle type are stable only when they have no more than three members.

## 5. Stochastic competition and cooperation

Finally, we consider the full mechanism (2.1*a–e*) in stochastic simulations. In total, stationary states are possible and, as said before, states are of particular importance: extinction, , selection with *n* different possible outcomes, with , and cooperation, . Here, we perform a complete analysis for the two smallest systems with and where we are dealing with four and eight stationary states, respectively.

First, we consider individual trajectories in the system with two subspecies (): four types of typical trajectories lead after an initial phase to one of the four states, , , and (figure 8). All four types of trajectories are very similar to those that were seen before in the pure competition or in the pure cooperation scenario. Because we chose the initial condition , all trajectories begin with a fast increase of . In the second part of the initial phase, random fluctuations determine the outcome of the process through canalizing the system into one of the four possible long-time states. In the long-time phase, the trajectories fluctuate around the deterministic stationary concentrations as illustrated in figure 8.

Again, two different evaluations of sampled trajectories are reported: (i) statistics of trajectories and (ii) final state counting. The expectation values , and for the competition-cooperation mechanism (2.1*a–e*) together with the one-standard deviation bands reflect the role of the initial phase in the determination of the long-term solution (figure 9). We were choosing two cases with different initial conditions: and . For , the random scatter originating from convergence to different quasi-stationary states is enormous, the one-error bands are extremely broad and those for and overlap almost completely. For , practically all trajectories go the state of cooperation , the one-error bands are much narrower and fulfil the expected relationship. The bands for and still overlap to some extent. The comparison of the deterministic curves with the sample mean values shows appreciable differences in the single copy case, but for 10 copies from each subspecies we find agreement within the line width. The results from sampling trajectories are put on a quantitative basis through final state counting (table 3). For very small numbers of the initially present members of the subspecies, for all , we find indeed a broad distribution over all states. The difference between one and two copies of the subspecies in the initial conditions is very impressive: in the example shown in table 3, the choice is already sufficient to end up in the state in more than 93% of all cases.

The three-membered system exhibits eight long-time or quasi-stationary states: the state of extinction , the three selection states with , the three exclusion states with and the cooperative state for . Different from the deterministic solutions, which always show *selection of the fittest* subspecies , every subspecies can be selected in the stochastic system with the fitness values determining the probabilities. The same is true for the exclusion states, where only , corresponding to *exclusion of the least fit* is stable for a range of values in the deterministic system, but all three exclusion states can appear as quasi-stationary long-time solutions of the stochastic process. The variables and fluctuate strongly at low particle numbers, and commonly the system passes over into one of the two selection states, or , respectively. At higher particle numbers the exclusion state is quasi-stationary once it has been reached in the initial random period (figure 10). For sufficiently high resources, , the stochastic system after it has successfully passed the random period approaches the cooperative state where the three particle numbers with fluctuate around the deterministic stationary values . Most trajectories progressing from competition towards cooperation show dominance of the fittest in the low concentration domain and cooperation in the high concentration range. This is precisely what we have also observed in the deterministic case (figure 1).

## 6. Discussion

The complexity of biology is commonly prohibitive for comprehensive quantitative modelling. The model presented here is to be understood as explorative [51]. It involves, of course, drastic simplifications of the real world, but the choice of the essential features is done in the same spirit as with the common models in population genetics and theoretical evolution. The model is to be understood as a kind of minimal system and there is ample room for making the system more complex in many aspects. Here, the focus is on the dynamics of reproducing populations and the open system is represented by a flow reactor of the CSTR type [29, p. 87ff.], which is easy to control in theory and accessible to experimental implementation. An apparent question concerns the appropriateness of the flow reactor for modelling external driving forces such as daily changes, seasons or slow climate changes. Tuning of the flow rate is easily introduced into the theoretical model and can also be achieved in experimental set-ups. In this sense, the flow reactor is quite flexible despite its simplicity. Reproduction is introduced as simple multiplication, , and mutation is neglected. Its consideration would not alter the general picture unless one were to assume that a subspecies that had died out is replaced readily through mutating another one, which is rather unrealistic. The meaning of catalysed reproduction, , needs some explanation. It may occur directly in the sense of a subspecies which is involved as a catalyst in the reproduction process [11]. The stoichiometry of the catalysed process may also result indirectly as the relevant overall reaction equation of complex mechanisms. Typical examples are known from chemistry: the essential reaction step of the Brusselator model, , can be seen as an abstraction from the Oregonator model, which in turn is a simplification of the Belousov–Zhabotinsky reaction with more than 20 individual reaction steps [52]. In reproduction, the presence of increases the fitness of in a such way that results in proportionality to the amount of available. Polynucleotide replication is a highly complex process involving a large number of individual steps even in the simplest cases [53]. Under suitable conditions, however, it shows the same overall kinetics as simple multiplication does and supports the usage of one-step overall kinetics in models.

It is important to comment on the notion of *cooperation* in the context of this contribution. The most common usage in biology, sociology and economics refers to the well-known *Prisoner's dilemma* game where cooperation implies a beneficial coordinated action of two prisoners with the risk that a contestant who chooses the non-cooperative strategy has a larger pay-off and the cooperative agent has a big loss. Game theory, in general, has been applied successfully to the evolution of animal behaviour that focuses on *evolutionarily stable strategies* (ESS) [54,55]. Two implicit assumptions are important for the relationship of game theory and biological evolution: (i) pay-off is related to fitness and reproductive success, and (ii) playing a given strategy has a genetic component and is inherited (for various biological and sociological applications of game theory see, for example, [56–58]). Behavioural game theory and replicator dynamics described by equation (2.4*b*) with , , and normalized concentrations ,
6.1are closely related indeed: it can be shown that every ESS corresponds to an asymptotically stable stationary point whereas the converse is not true, and there are stable attractors in the dynamical system (6.1) to which no ESS can be found [59–61]. The ecological May–Leonard model [62] may serve as an example (for the relationship between dynamical systems and ESSs, see also [63–65]). Hypercycles correspond to ESSs and in this way cooperation in game theory in related to cooperation in the sense it is used here: competition is described by independent reproduction and additional terms are added to reproduction kinetics that can suppress competition and introduce cooperation. On the molecular level, catalysis is a proper tool for this goal. Coupling between reproducing elements creates a new functional unit, species become subspecies and Gause's law does not apply any more on the lower hierarchical level. It applies, however, to the new functional units.

An attractive feature of the model is the description of selection and cooperation in a single dynamical system. There is a transition range that can be detected particularly easily in case the deterministic cooperative system is oscillating (; figure 7). One way to progress from the selection regime to cooperation is to simply increase the concentration of A in the stock solution and indeed there is a critical value of above which the long-time solution starts to oscillate. Below the critical value, we observe selection, which may be biased by differences in the values. Increasing the amount of A available for reproduction is tantamount to progressing from scarcity to abundance of resources. If we agree to identify the state of cooperation as the beginning of higher-level organization, the toy model apparently supports the initially made conjecture that plenty of resources are required for a transition to higher-level organization.

The formation of an organized system of functionally coupled replicators is only the beginning of a major transition. The next logical step concerns the creation of a *unit of selection* at the new and higher hierarchical level [1,2,66]. This is done through the formation of boundaries between the new functional unit and the environment. Commonly, spatial boundaries are formed in biology; examples are cell membranes, cell walls and various kinds of skins. The members of cooperative units in animal and human societies use signals and language to distinguish members from non-members. The formation of boundaries serves two goals: (i) it saves the new unit from exploitation by parasites, which take advantage of the new achievements but do not contribute to the commons, and (ii) it opens the route to Darwinian evolution at the new hierarchical level. In addition, the evolution of cooperative groups requires *mutual policing* and *repression of competition* [67]. A large number of investigations were and are modeling cooperation by means of game theory [56].

An illustrative example of a transition is provided by a model for the origin of life [9,11,66]: RNA molecules when reproducing independently compete for the resources required for replication and are subjected to Darwinian natural selection. Uncatalysed replication, however, is inefficient and a successful system depends upon the presence of catalytic RNA molecules called *ribozymes*. Because of the low replication fidelity, an error propagation problem known as *error threshold* [68] arises and sets a limit to the lengths of molecules that can be copied faithfully. A more accurately operating replicase ribozyme cannot be made in one piece and several RNA molecules have to *cooperate* in the formation of the catalyst [69–71]. The transition creates a new functional unit in which the RNA molecules are no longer independent. Together they cooperate in efficient reproduction and the simplest coupling mechanism for such a composite replicase is a hypercycle. From the point of genetic information the integration of several information carriers into a single functional unit increases the content of information that can be stored, a more accurate replication machinery is formed that allows for the replication of longer sequences and the error threshold experienced by the original short molecules is circumvented.

Stochastic effects make an important contribution in mainly two features. Depending on the random events in the initial phase the stochastic trajectories approach either the state of extinction , one of the selection states or the cooperative state . As expected, we find the broadest scatter with a single copy from each subspecies, because then the probability of extinction in the early phase is largest. For a population size and a flow rate of , two to three copies of each subspecies, , are already sufficient for obtaining more than 90% outcome in agreement with a solution corresponding to the deterministic result. In the calculation of expectation values and standard deviations, we are dealing with two kinds of fluctuations: (i) the natural fluctuations caused by thermal noise, which are around and (ii) the scatter due to the distribution over the different attractors, which as said becomes already small when two or more copies of each replicator are present. We must keep in mind that any more efficient variant is created by mutation; this is a singular event and hence every new subspecies has to begin with .

The second major impact of stochasticity concerns cooperative catalytic cycles with more than three members, . For , the systems oscillate and the deterministic solutions converge to a stable limit cycle. The stochastic process goes into self-enhancing oscillations until one subspecies disappears and then the entire system goes extinct. The case is more subtle: The deterministic solution shows marginal stability of the central stationary state in the linear approach and is stabilized by second-order terms. The stochastic process with develops self-enhancing oscillations even when the stationary values of the deterministic system are taken as initial conditions and becomes extinct by the same mechanism as discussed for . For very large systems, the time to extinction might be long. In the stochastic model, only two- and three-membered cooperative units can reach quasi-stationary states, four-membered systems are unlikely to last long, and units with five members and more go extinct.

## 7. Conclusion

Coming back to the questions we have raised at the beginning of this contribution, we can give answers within the frame of our simple model. The first problem was the possibility of conceiving a model that is able to describe the transition from Darwinian type selection to cooperation between former competitors. The answer is simply *‘yes’*, because the simple model containing first- and second-order terms of the form and progressed from a phase of selection towards a dynamically organized cooperative state. Considering the concentration of the nutrition compound A as a measure for the availability of resources, we found indeed that small resources leading to small values of result in selection, whereas plenty of resources leading to large values of favours the cooperative state.

Stochasticity introduces substantial modifications into the deterministic picture of evolution from selection to cooperation. There is a phase at very low particle numbers where random events dominate and the system can be driven to practically every possible outcome. Accordingly, we can only derive a probability distribution for the various long-time results. As every new variant created by mutation starts from a single copy, the initial random phase is crucial. The good news, however, is: provided enough resources are available—the reservoir of A is sufficiently rich corresponding to a sufficiently large value of —and the time constraint is not too strong—the flow rate *r* is not too high—a few copies of every subspecies are enough for a reasonably large probability of nucleation of the cooperative unit.

The number of cooperating subspecies forming a stable unit is seriously limited at least in case of the simple hypercycle mechanism. Units of four or more subspecies, , are unstable because of self-enhancing oscillations that wipe out the system sooner or later. The self-enhancement is a result of autocatalysis: in order to become effective the increase has to complete one cycle, which is tantamount to a delay becoming longer for larger cycles and delay equations have a strong tendency to oscillate. In the deterministic system, concentrations can be arbitrarily small but the stochastic process is limited by discrete numbers and less than one particle means extinction.

It is suggestive to use both results as hints for finding answers to the initially raised questions but one must not forget that we have considered reproduction kinetics in a largely simplified manner here. Needless to repeat, reproduction in nature is an enormously complex process. In major transitions, cooperative reproduction is only a first step and many other complex processes are in operation too [1]. Nevertheless, we can summarize: abundance of resources is a requirement for transitions. There is a strong preference for the formation of small symbiotic units with only two or three elements and this is supported by dynamical features from reproduction kinetics.

## Competing interests

I declare I have no competing interests.

## Funding

The work reported here has been supported by the University of Vienna and the Santa Fe Institute.

## Footnotes

One contribution of 13 to a theme issue ‘The major synthetic evolutionary transitions’.

↵1 The subscript of the states represents the number of subspecies, which are present at the stationary state.

- Accepted March 3, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.