Elias Zintzaras

Mauro Santos
email: mauro.santos@uab.es

Eörs Szathmáry
email: szathmary@colbud.hu

Selfishness versus functional cooperation in a stochastic protocell model

Keywords: Artificial cells, Functional complexity, Monte Carlo methods, Qβ replicase, Ribozymes, Stochastic corrector model

How to design an "evolvable" artificial system capable to increase in complexity? Although Darwin's theory of evolution by natural selection obviously offers a firm foundation, little hope of success seems to be expected from the explanatory adequacy of modern evolutionary theory, which does a good job at explaining what has already happened but remains practically helpless at predicting what will occur. However, the study of the major transitions in evolution clearly suggests that increases in complexity have occurred on those occasions when the conflicting interests between competing individuals were partly subjugated. This immediately raises the issue about "levels of selection" in evolutionary biology, and the idea that multi-level selection scenarios are required for complexity to emerge. After analyzing the dynamical behaviour of competing replicators within compartments, we show here that a proliferation of differentiated catalysts and/or improvement of catalytic efficiency of ribozymes can potentially evolve in properly designed artificial cells where the strong internal competition between the different species of replicators is somewhat prevented (i.e., by choosing them with equal probability). Experimental evolution in these systems will likely stand as beautiful examples of artificial adaptive systems, and will provide new insights to understand possible evolutionary paths to the evolution of metabolic complexity.

Introduction

At the turn of the 21 st century we are witnessing an ambitious scientific program to synthesize artificial cells or protocells that capture the essentials to be considered "alive" [START_REF] Szostak | Synthesizing life[END_REF]Pohorille and Deamer, 2002;[START_REF] Deamer | A giant step towards artificial life?[END_REF][START_REF] Luisi | The Emergence of Life: From Chemical Origins to Synthetic Biology[END_REF][START_REF] Luisi | Approaches to semi-synthetic minimal cells: a review[END_REF][START_REF] Mansy | Reconstructing the emergence of cellular life through the synthesis of model protocells[END_REF][START_REF] Rasmussen | Protocells. Bridging Nonliving and Living Matter[END_REF]. This move immediately raises the question: what are the criteria for minimal life? From a bottom-up point of view a minimal living entity can be conceptualized as a chemical system comprising three subsystems: a metabolic network, template replication, and a boundary structure [START_REF] Gánti | Az Élet Princípuma (The Principle of Life)[END_REF][START_REF] Gánti | The Principles of Life[END_REF][START_REF] Szathmáry | Evolutionary potential and requirements for minimal protocells[END_REF]. From a top-down perspective ("minimal cell" project;

see [START_REF] Luisi | Chemical aspects of synthetic biology[END_REF] the problem is more related to the question of what is the minimum level of complexity we can attain in actual living systems without losing crucial properties. In any case, we concur with [START_REF] Luisi | Approaches to semi-synthetic minimal cells: a review[END_REF] that such entities have to satisfy three basic requirements: self-maintenance, self-reproduction, and evolvability (i.e., the ability of a population of entities to generate diversity and experience Darwinian evolution). Implicit here is the notion that the three subsystems in [START_REF] Gánti | Az Élet Princípuma (The Principle of Life)[END_REF] chemoton concept can be combined to yield three different doublet infrabiological systems [START_REF] Szathmáry | Evolutionary potential and requirements for minimal protocells[END_REF].

One of the simplest constructs that (apparently) meets the criteria for minimal life is the so-called "RNA-cell", a purely imaginary object containing in a vesicle two ribozymes:

one with replicase activity, and the other catalyzing the synthesis of membrane components [START_REF] Szostak | Synthesizing life[END_REF][START_REF] Luisi | Approaches to semi-synthetic minimal cells: a review[END_REF]. Although this system is based on hypothetical ribozymes and, therefore, still poses formidable challenges associated with its eventual assembly, current technology allows some alternatives in the realm of infrabiological systems. For instance, artificial cells that allow for RNA replication -using the RNAdirected RNA replication enzyme Qβ replicase -and vesicle division have already been assembled [START_REF] Oberholzer | Enzymatic RNA replication in self-reproducing vesicles: an approach to a minimal cell[END_REF]. These infrabiological systems are completely devoid of metabolism, but they are however capable of replication and mutation as required for Darwinian evolution, with a caveat: the enzyme and RNA molecules in [START_REF] Oberholzer | Enzymatic RNA replication in self-reproducing vesicles: an approach to a minimal cell[END_REF] construct are not reproduced from inside and, therefore, will be eventually diluted and cause the "death" of the system (see [START_REF] Luisi | The Emergence of Life: From Chemical Origins to Synthetic Biology[END_REF][START_REF] Luisi | Chemical aspects of synthetic biology[END_REF]. Assuming that this drawback can be overcome in some infrabiological systems, evolution experiments with populations of vesicles can in principle be done, but the question is: what could be achieved with these

experiments?

The answer is that nothing very exciting unless the functional complexity of RNAs entrapped in the vesicles is allowed to increase, and vesicle growth and reproduction are somewhat linked to that complexity (i.e., multi-level selection has to be imposed). The reason is simply because natural selection acting exclusively on individual replicators fosters the evolution of molecular parasites that are better targets for replication but do not contribute to the "common good" (see [START_REF] Scheuring | Avoiding catch-22 of early evolution by stepwise increase in copying fidelity[END_REF][START_REF] Szabó | In silico simulations reveal that replicators with limited dispersal evolve towards higher efficiency and fidelity[END_REF]Takeuchi andHogeweg, 2007, 2009;[START_REF] Branciamore | The origin of life: chemical evolution of a metabolic system in a mineral honeycomb?[END_REF], as already illustrated by the classical studies of in vitro evolution of RNA molecules carried out by Spiegelman and his colleagues [START_REF] Mills | An extracellular Darwinian experiment with a self-duplicating nucleic acid molecule[END_REF][START_REF] Spiegelman | An approach to the experimental analysis of precellular evolution[END_REF]. These authors isolated and purified the more than 4,000 nucleotides-long single-stranded Qβ RNA that encodes a number of proteins, including Qβ replicase. When this RNA was added to a solution containing Qβ replicase and energy-rich nucleotide triphosphates, new infectious RNA strands were synthesized. However, when they kept the RNA population in perpetual growth by using the technique of serial transfer, competition between RNA strains was just for resources and replication: the RNA molecules replicated, with errors, faster and faster. As a correlated response these RNAs evolved shorter sequences and lost their infectiousness.

The two lessons to be learned from these experiments are [START_REF] Bell | The Basics of Selection[END_REF]: (i) that the rate of self-replication is the only attribute that can be selected directly; and (ii) that in contrast to the in vitro experiments bacteriophage Qβ cannot get rid of most of its genome in its natural state because it would be unable to infect bacteria.

An additional warning is that any successful technological construct that meets at least one of the first two requirements for minimal life has also to deal with the associated evolutionary hurdles. In other words, it is not enough for the hypothetical RNA-cell above to be capable of self-maintenance and self-reproduction; it also has to persist in the long run. For instance, real cells have an organized cell division and duplicate themselves with the same genetic content, but this is difficult to implement in artificial cells where duplication occurs through purely physico-chemical forces and the parental material is randomly transmitted to progeny. In this context, it might be rewarding to also consider the lessons learned from theoretical models that impose a spatial structure to a population of replicators by encapsulating them into vesicles (compartments, protocells). [START_REF] Niesert | Origin of life: Between Scylla and Charybdis[END_REF] were the first to propose a vesicle model, which was further elaborated by [START_REF] Szathmáry | Group selection of early replicators and the origin of life[END_REF] who described the "stochastic corrector model" (SCM; see also [START_REF] Grey | A re-examination of the stochastic corrector model[END_REF][START_REF] Zintzaras | Living' under the challenge of information decay: the stochastic corrector model vs. hypercycles[END_REF]. The initial stimulus behind vesicle models was to solve the conundrum of the evolutionarily dynamic coexistence of unlinked genes.

The dynamical behaviour of the system depends on two types of stochasticity: (i) replication of templates within protocells, and (ii) random assortment of templates into offspring protocells. Even though templates compete within compartments, selection on stochastically produced offspring variants (between-protocell selection) can rescue the population from extinction, which reaches equilibrium with a constant frequency of the optimal protocell. Gene redundancy is necessary to avoid an unsupportive assortment load;

that is, the drop in average fitness due to the random loss of any essential template after stochastic fission of templates in the two daughter protocells.

So far theoretical considerations have mainly stressed the sloppiness of vesicle models due to the complex balance among redundancy, assortment load, and mutation load [START_REF] Santos | Recombination in primeval genomes: a step forward but still a long leap from maintaining a sizeable genome[END_REF][START_REF] Fontanari | Coexistence and error propagation in package models: a group selection approach[END_REF], which causes protocells' survival to dwell somewhere between two Homer's monsters: Scylla and Charybdis [START_REF] Niesert | Origin of life: Between Scylla and Charybdis[END_REF]; see also [START_REF] Niesert | How many genes to start with? A computer simulation about the origin of life[END_REF]. Under this scenario it is difficult to envisage that vesicle models could have any potential for evolutionary novelties besides the proper balance between the two levels of selection for just long-term persistence of the population. However, we show here that under some conditions novel evolutionary directions can emerge in artificial cells through the generation of functional diversification. In other words, these systems can be evolvable because there is the potential for acquiring novel functions through genetic changes.

The remainder of the paper is organized as follows. Firstly, we somewhat depart from previous static analyses that have neglected a continuous evolution of both replication rate and functional activity of the templates. Since any group selection in the form of vesicle growth and reproduction is likely antagonistic with the RNA replication rate, as illustrated by the Qβ situation, a trade-off between replication rate and functional activity likely arises. Secondly, we show that the dynamical analysis of the interaction between the two levels of selection in the standard theoretical approach offers some insights on how functional diversification of templates could be eventually achieved in vesicle models.

Thirdly, in order to leave open the possibility for artificial cells to increase in functional complexity we claim that a deviation from the standard kinetics of the SCM has to be realized somehow, and suggest ways to do that. These considerations allow us to conclude that a proliferation of differentiated catalysts and/or improvement of catalytic efficiency of ribozymes will be likely observed in evolution experiments with artificial cells, and will stand as a beautiful example of artificial adaptive systems.

A simple trade-off model of two-level selection

We use here the basis behind the Monte Carlo implementation of the SCM prompted by [START_REF] Zintzaras | Living' under the challenge of information decay: the stochastic corrector model vs. hypercycles[END_REF] and [START_REF] Santos | Recombination in primeval genomes: a step forward but still a long leap from maintaining a sizeable genome[END_REF]. Each vesicle is assumed to consist of two types of RNAs whose joint functions are essential for growth and eventual vesicle splitting.

The enzyme Qβ replicase is assumed to catalyze the replication of RNAs in a manner similar to Qβ phage with tRNA-like 3' genome tags (i.e., a recognition site for the replicase at the end of the template; see [START_REF] Schaffner | Nanovariant RNAs: nucleotide sequence and interaction with bacteriophage Qβ replicase[END_REF]. Therefore, the RNA templates ( 1 T and 2 Τ , in multiple copies each) are organized as having a target region that defines an average affinity toward the replicase (i.e., whether they are good substrates for the replicase) plus a sequence of nucleotides involved in their function (enzyme efficiency). The problem now is to devise a set of rules that capture the essence of a twolevel process of selection, where selection acts directly on rates of replication and indirectly on enzyme efficiency and protocell grow, after imposing a trade-off between target affinity ( τ ) and enzyme efficiency ( ε ); see e.g., [START_REF] Könnyű | Prebiotic replicase evolution in a surfacebound metabolic system: parasites as a source of adaptive evolution[END_REF]. The trade-off function we chose was a typical concave relationship ( ) ( )

1 1 τ ε ε = - + ( ) ( ) 1 1 ε τ τ = - + ⎡ ⎤ ⎣ ⎦ ; with 0 01 0 98 . ε, τ . ≤ ≤
. But we emphasize that the qualitative results presented here are robust for different trade-off functions ( ) ( )

1 1 p τ ε ε = - + ( ) ( ) 1 1 p ε τ τ ⎡ ⎤ = - + ⎣ ⎦ , with "realistic" values of p ( 1 5 p -≤ ≤ ).
Let us now consider a finite population of K vesicles enclosing T and 2 Τ at 0 t . In each time step, a vesicle is randomly chosen according to its relative fitness for template replication. At the vesicle level, the fitness function we used was:

1 1 1 , 1 i d i T d i i n w ε = = = ∑ ∑ ( ) 1
where d is the number of different templates (ribozymes), and

1 T i n i i j j ε ε = = ∑ is the total
contribution of all its j copies to vesicle fitness:

T ij 2 2 max max 2 max 1 1 exp( ) exp( ) , 1 1 exp( ) ij ε ε ε ε ε ⎛ ⎞ ⎜ ⎟ - ⎜ ⎟ - ⎝ ⎠ = ⎛ ⎞ - ⎜ ⎟ ⎝ ⎠ ( ) 2 with T ij
ε being the actual enzymatic efficiency of the jth copy of type i ribozyme, and max 1 ε = the maximum enzymatic efficiency. The rationale behind the construction of the fitness function, which was first used in essentially this form by [START_REF] Szathmáry | What determines the size of the genetic alphabet?[END_REF] see also [START_REF] Zintzaras | Living' under the challenge of information decay: the stochastic corrector model vs. hypercycles[END_REF][START_REF] Santos | Recombination in primeval genomes: a step forward but still a long leap from maintaining a sizeable genome[END_REF] for the fitness of a ribo-organism, relies on the assumption that fitness isas usual for microbesessentially determined by the flux F of a pathway unsaturated by enzymes. Metabolic control theory shows that for such a case [START_REF] Kacser | The control of flux[END_REF]

: 1 , 1 n i i C F ε = = ∑ (3)
where C is a constant (set to 1) and the enzymatic efficiency of ribozyme i involves enzyme concentrations and catalytic constants. Catalytic efficiency is an exponential function of the binding strength between the enzyme and the substrate of the catalyzed reaction (Kacser and Beebe, 1984).

However, as it stands the numerator in Eq. ( 1) "imposes" to write down in the algorithm that protocell fitness is set to zero if any 0 i T n = . Although this follows our rationale that the death of vesicles happens whenever they lack an essential template, it would seem more appropriate to avoid the ad hoc assumption by making this numerator more similar to C in Eq. ( 3). Below (section 3) we discuss our numerical results assuming different fitness functions at the vesicle level.

Once a vesicle has been chosen, a random template is replicated according to its replication probability, which obviously depends on its target affinity towards the ( Qβ ) replicase (see below). If the number of templates in the vesicle is below twice the initial number at 0 t , the step ends; that is, vesicle's fitness is updated following Eq. ( 1), it is turned back to the population, and the next step starts by randomly chosen a vesicle according to its relative fitness. Otherwise, the vesicle splits and templates are randomly assorted into two daughter vesicles. One offspring replaces the parental vesicle and the other a randomly chosen one from the population. Our protocol is, therefore, based on the classical Moran process [START_REF] Moran | Random processes in genetics[END_REF][START_REF] Ewens | Mathematical Population Genetics[END_REF].

A critical feature in our Monte Carlo method is template replication. Remember, we have two conflicting levels of selection. Within-vesicle selection always selects for high target affinities and, hence, poor enzyme efficiency because of the imposed trade-off.

However, between-vesicle competition imposes a compromise between these two traits because the realized vesicle growth depends on both enzyme efficiency (increases vesicle fitness) and target affinity (increases the number of templates and, hence, growth rate as defined in the simulation process). Replication of the parental template would introduce mutations affecting the enzyme efficiency of the copy largely independently of its target affinity. In the Qβ RNA in vitro experiments infectiousness was lost mainly by the shortening of RNA strands to a sixth or so of the original size, but we can ignore here these specific details and simply assume an average reduction in the copy's enzyme efficiency.

Suppose the efficiency of the parental template (e.g., P 1 T ) is P ε , from which we want to obtain the efficiency of the copy C ε ι . The following function was conveniently used:

( ) ( ) ( ) (
) sum up, the average mutation slightly decreases enzyme activity but neutral or compensatory mutations are allowed, thus somewhat mimicking the effects of mutations on the enzymatic activity of real ribozymes [START_REF] Kun | Real ribozymes suggest a relaxed error threshold[END_REF].

max P P C P P 1 2 Sgn δ ε ε Sgn δ δε ε ε , k δε ⎡ ⎤ + - ⎣ ⎦ = + + ι ( ) 4
Assume now that the target affinity towards the replicase for the same parental template P 1 T is P τ , which can also mutate to C τ ι in the copy. A similar function to Eq. ( 4)

was used: high input of deleterious mutations and, therefore, a substantial drop in average fitness that could threaten the eventual survival of the vesicle population. However, recall that here we are not interested in the issue of whether or not vesicle models can overcome the information crisis of prebiotic evolution because of the bottleneck imposed by the error-threshold (see [START_REF] Zintzaras | Living' under the challenge of information decay: the stochastic corrector model vs. hypercycles[END_REF][START_REF] Stadler | Molecular replicator dynamics[END_REF][START_REF] Santos | Recombination in primeval genomes: a step forward but still a long leap from maintaining a sizeable genome[END_REF][START_REF] Silvestre | Package models and the information crisis of prebiotic evolution[END_REF], but in the evolutionary dynamics of artificial cells where templates are replicated by highly evolved enzymes as (e.g.) the Qβ replicase. Since with this enzyme the mutation rate per base per replication would be 3 10 μ -≈ [START_REF] Drake | Rates of spontaneous mutation among RNA viruses[END_REF], it could in principle be possible to encapsulate up to 100 different ribozymes of sequence length 70 nucleotides each [START_REF] Kun | Real ribozymes suggest a relaxed error threshold[END_REF]. . This is, however, not necessarily true and strongly depends on the redundancy levels of templates 1 T and 2 T (see Appendix). For instance, Fig. 2 plots two different scenarios: in the first case (Fig. 2a)

( ) ( ) ( ) ( ) max P P C P P 1 2 Sgn γ τ τ Sgn γ γτ τ τ , k γτ ⎡ ⎤ + - ⎣ ⎦ = + + ι ( )
20 i T n =
for both 1 T and 2 T , whereas in the second case (Fig. 2b)

1 38 T n = and 2 2 T n = .
With equal levels of redundancy vesicle growth is maximized, and mutational load [START_REF] Haldane | The effect of variation on fitness[END_REF][START_REF] Crow | Genetic loads and the cost of natural selection[END_REF] is minimized when

T1 T2 , , τ ε τ ε = .
With unequal levels of redundancy, vesicle growth increases with higher target affinities of the most abundant template, but this causes a high mutational load which is, conversely, minimized with increasing levels of enzyme efficiency for the same template. Note also that the minimum mutational load with unequal levels of redundancy

1 1 p τ ε ε = - + ( ) ( ) 1 1 p ε τ τ ⎡ ⎤ = - + ⎣ ⎦ with
"realistic" values of p (see above). To summarize, Fig. 2 and the Appendix show that the actual net growth of the vesicle, which is obviously a function of the combined dynamics of metabolism and total template replication, is not maximized when T1 T2

, , τ ε τ ε = .

However, deviation from an even concentration of both templates results in a higher assortment load and, hence, selection to maximize vesicle growth is opposed by selection at the vesicle level because it results in a higher probability of lineage extinction.

We now digress slightly about our fitness function in Eq. ( 1) and the "written" assumption that protocell fitness is set to zero if any 0 i T n = , that is, the death of vesicles happens whenever they lack an essential template. Simulations were also run by setting the numerator in Eq. ( 1) equal to 1 (see Eq. ( 3)), and by assuming a multiplicative fitness function

T T 1 1 ,1 ,2 1 1 n n i j i j w ε ε = = = * ∑ ∑
. In both cases vesicle's fitness is 0 if enzymatic efficiency is 0. Numerical results were qualitatively similar to those already reported, which clearly suggests that the way we choose vesicles (notice that Eq. ( 1) gives more weight to the rarest template type) seems to be less important than template competition for replication and eventual growth and division of vesicles. The reason for this can perhaps be better appreciated in Fig. 2b (left panel) where cell growth positively increases as a function of the target affinity of the more abundant template type.

Although the preceding analyses suggest that there is some room for divergence in templates' target affinities and enzyme efficiencies, the question naturally arises: could it be possible for templates 1 T and 2 Τ to set apart into two different clusters according to target affinities ( τ ) and enzymatic efficiencies ( ε ) if the inherent strong internal competition in vesicle models is somewhat alleviated to reduce the assortment load? If so, how could it be implemented when assembling artificial cells?

In principle, assortment load could be alleviated by assuming high levels of redundancy at 0 t , and Fig. 3 ). But the problem with high levels of redundancy is that evolvability could be prevented because of the increased risk that Darwinian selection would be stopped because of dilution of favourable mutations in an "orgy of redundancy" [START_REF] Koch | Evolution vs the number of gene copies per primitive cell[END_REF]. We show next that evolution experiments with artificial cells can potentially combine the best of both worlds: template differentiation and low assortment low while keeping redundancy at reasonably low levels for positive selection to happen.

Within-template competition

Let us imagine a way to prevent the strong internal competition in vesicle models. For instance, assume that during the replication process the Qβ replicase first "chooses" at random between templates 1 T and 2 Τ with equal probability, and then replicates one of the i T n copies according to its replication probability; that is, assume that strong internal competition occurs only within the same template type. The outcome now is a qualitative change in the internal dynamics of standard vesicle models that results in the differentiation of templates 1 T and 2 Τ into two clusters in accordance with target affinities and enzymatic efficiencies, and this happens largely independent of initial levels of redundancy (as long as redundancy is not too low; say 3 i T n ≤ at 0 t ). For instance, Fig. 4 plots the combination of values where the final population of vesicles with initially

6 20 i T n , =
copies for each template stabilizes (c.f. with Fig. 1). In the particular case when = with 32.2 copies. Note, however, that functional divergence seems less pronounced with high levels of redundancy (c.f. Figs. 1 and4), which is probably due to the "cooperative" dynamics of template growth imposed by the "withintemplate competition". In other words, there seems to be a complex balance between cell growth, which selects for unequal templates' concentration, and "cooperative" dynamics, which selects in the opposite direction and pushes T1 , τ ε to become equal to T2 , τ ε because this maximizes cell growth if template concentrations are equal (Fig. 2 and Appendix).

6 i T n = ,
However, for all situations investigated, at some time in the evolution of the vesicle population we observed that templates 1 T and 2 Τ set apart into two different clusters according to target affinities and enzymatic efficiencies. This happens only when we impose a two-level selection (within-template within-vesicle selection, and betweenvesicle selection); otherwise both kinds of templates settle down at similar values (results not shown). To sum up, assuming that there is a trade-off between target affinity and enzyme efficiency as suggested by the in vitro experiments with bacteriophage Qβ , substantial functional divergence between templates could be expected if the strong internal competition is somewhat alleviated (assuming initial redundancy is not too high).

But the obvious question is: how is this achieved in practical terms? There are two potential answers to this question. The first is related to the basic growth dynamics of replicators. The second to the fact that the preceding analysis equates vesicle growth to copy number and ignores that RNA replication produces ( ) + and ( )

-strands, where one strand, say ( ) + , is the functional ribozyme and the ( ) strand is the template ("gene").

In the field of prebiotic evolution non-conventional growth laws, such as hyperbolic and parabolic, have been widely discussed [START_REF] Von Kiedrowski | Minimal replicator theory i: parabolic versus exponential growth[END_REF][START_REF] Szathmáry | From replicators to reproducers: the first major transitions leading to life[END_REF]. Both represent departures from simple Malthusian growth: they are faster and slower than it, respectively. Parabolic growth was experimentally demonstrated to happen with small synthetic replicators [START_REF] Von Kiedrowski | A self-replicating hexadeoxynucleotide[END_REF], and its consequences for selection in a competitive setting are remarkable: survival of everybody [START_REF] Szathmáry | Sub-exponential growth and coexistence of nonenzymatically replicating templates[END_REF]. Inhibition of RNA synthesis by Qβ replicase can occur, probably due to competition between ( ) + and ( ) strands for the same enzyme molecules [START_REF] Kondo | Inhibition of Qβ replicase by excess template[END_REF]. Thus, it is known that a fraction of ( ) + and ( ) strands anneal to double-stranded form even in natural RNA replication today [START_REF] Biebricher | Kinetics of RNA replication: plusminus strand asymmetry and double-strand formation[END_REF]. This self-inhibition could lead to some coexistence, but in a very wasteful manner. A more interesting way of dynamical coexistence might be achieved by niche differentiation of the templates: obviously, if the two types of templates consisted of A:U and G:C pairs, respectively, they would not compete for the same nucleotides. Although ribozymes composed of one of these two pairs of nucleotides may be possible (see [START_REF] Reader | A ribozyme composed of only two different nucleotides[END_REF] this extreme case cannot hold generally. Realistic situations could nevertheless alleviate competition through niche differentiation within protocells. This suggestion warrants further analysis.

There is some experimental support for the assumption that vesicles split when their template concentration reaches a given threshold (whatever that may be in real vesicles), as assumed in our model. Thus, Chen et al. ( 2004) have shown that membrane growth and eventual protocell division can be driven by the osmotic pressure exerted by replicating RNA molecules encapsulated in fatty acid vesicles, and suggested that faster replication would lead to faster vesicle growth and fitter protocells. However, from the standard kinetics of the SCM we already know that faster vesicle growth cannot always be easily equated to vesicle fitness (see Appendix). Experimental evolution studies using artificial cells would mostly rely on the enzymatic function of the ( ) + strands (assumed here to be the functional ribozyme) to assess vesicle fitness. Therefore, if cell division is determined by some metabolic product other than the total number of templates, then selection for better enzymes, as well as for a minimum concentration of coexisting enzymes whose collective catalysis reaches a given threshold criterion for selection, will be guaranteed and functional divergence between templates will be a likely outcome. The major evolutionary increases in complexity have occurred on those occasions when the conflicting interests between competing individuals were partly subjugated (Maynard [START_REF] Smith | The Major Transitions in Evolution[END_REF], and it is probably not far-fetched to think of the within-template competition scenario above as a suitable form of "cooperative" dynamics in actual vesicles. We therefore believe that some fascinating findings in evolution experiments with artificial cells are waiting for us in the coming years.

5.

Concluding remarks [START_REF] Takeuchi | Multilevel selection in models of prebiotic evolution II: A direct comparison of compartmentalization and spatial self-organization[END_REF] have recently studied the evolution of a RNA-replication system (consisting of a parasite and a replicase) where multi-level selection was incorporated by simulating either (i) spatial self-organization in two-dimensional aggregates that constraint diffusion; or (ii) explicit compartmentalization as assumed here.

Parasites could switch between two conformational states: one increasing replication rate and the other decreasing replication rate and facilitating the vesicle growth. Interestingly, their model discovered an emergent trade-off for those conformational states due to the complex interaction between the two-levels of selection (within-and betweencompartments). Our present model is somewhat complementary in that we already assumed a trade-off between target affinity and enzyme efficiency (as illustrated by the in vitro experiments with bacteriophage Qβ ), and interactions between those properties arise.

We therefore concur with Takeuchi and Hogeweg's conclusion: the complex dynamics of vesicle models can indeed generate novel evolutionary directions and increase complexity in RNA replicator systems.

We think our model can also be taken as a complementary solution for the evolution of functional diversification to that proposed by [START_REF] Kacser | Evolution of catalytic proteins or on the origin of enzyme species by means of natural selection[END_REF] see also [START_REF] Beeby | Metabolic constraints in evolution[END_REF]. These authors focused on the functional properties of catalytic proteins to advance a kinetic-based mechanism where an initial enzyme with broad substrate specificity leads to a collection of enzymes with differentiated specificities.

However, their simulation experiments with protocell populations [START_REF] Beeby | Metabolic constraints in evolution[END_REF] obviously ignored the potential of genetic conflicts simply because they did not deal with replicable information carriers, and the kinetics of vesicle growth was the only selection criterion. Nevertheless, complexification of protocell genomes, based on duplication and divergence (cf. Maynard [START_REF] Smith | The Major Transitions in Evolution[END_REF] and parasite taming [START_REF] Könnyű | Prebiotic replicase evolution in a surfacebound metabolic system: parasites as a source of adaptive evolution[END_REF] should be investigated in detail in the future.

The criterion that vesicle division happens when a critical level of total template number has been reached apparently comes somewhat close to "soft group selection", also analyzed by [START_REF] Traulsen | Stochastic evolutionary dynamics on two levels[END_REF], where individual and group selection favour the same variants, since faster-growing templates yield earlier vesicle division (see our analysis in the Appendix). However, soft group selection does not prevail since protocells lacking one of the templates cannot grow at all. Also, the metabolic function renders template growth frequency-dependent, not considered by [START_REF] Traulsen | Stochastic evolutionary dynamics on two levels[END_REF].

Inspired by the natural model of cellular compartmentalisation, [START_REF] Agresti | Selection of ribozymes that catalyze multiple-turnover Diels-Alder cycloadditions by using in vitro compartmentalization[END_REF] developed a system termed in vitro compartmentalization (IVC). IVC uses compartmentalization to link genotype (a nucleic acid that can be replicated) and phenotype (a functional trait such as a binding or catalytic activity). Instead of compartmentalizing genes in cells, as in nature, in IVC the genes are compartmentalised in aqueous microdroplets dispersed in a water-in-oil emulsion. ICV has certain limitations that can be overcome by making and manipulating droplets in digital microfluidic systems [START_REF] Whitesides | The origins and the future of microfluidics[END_REF]Herold and Rasooly, 2009). We believe that droplet-based microfluidic "evolution machines" will constitute a powerful tool for quantitative studies on ribozyme evolution, and will enable us to test to what extent it is possible the proliferation of differentiated catalysts and/or to improve the catalytic efficiency of ribozymes. This will, in turn, help our understanding of how biochemical complexity could have arisen during the early steps in the origin of life.

growth, in turn favouring high replication rate and low metabolic efficiency for the abundant template. Thus one should expect spontaneous symmetry breaking, depending on which protocell lineage one follows (note the symmetry of Fig. A2a). Of course these plots should be combined multiplicatively with the internal metabolism of the protocell. Fig. A2b shows the combined plot for 1 2 T n = . The total template growth is maximized for high and low template affinity (or for low and high enzymatic efficiency) for templates 2 Τ and 1 T , respectively. The conclusion that the "richer template gets richer" (Matthew effect) is maintained, i.e. selection would favour increased competitiveness of the template that is already abundant. As discussed in the main text, these plots do not take into account the fact that cells with only one template type are completely unviable. This is the only effect that prevents the spread of selfish templates in the end. At high redundancy this strong punishment is "not felt" by a large fraction of the compartment population, hence one expects some stochastic symmetry breaking in which, across the whole protocell population, one template type would grow better than the other. Indeed this is what we find (main text).

Legends for figures

  where δ is the raw change in enzyme activity sampled from a normal distribution e., random mutations are expected to be mostly detrimental),( )Sgn δ is the sign of δ ( 1 ± ), and k is a smoothness parameter. When 0 δ > , for any of P ε value the activity of the copy C P

5

  where γ is the raw change in target affinity now sampled from ensure that the final target affinity and enzyme efficiency of the copy C C , τ ε remain in the permissible region imposed by the trade-off, and allowed a suitable shift θ a random angle in radians. The simulation programs were implemented in the MATLAB algebra program environment (V7; MathWorks, 2005), and in Compaq Visual Fortran90 (2000) using the IMS library.3. Between-template competition (standard kinetics of the SCM)Before addressing the net consequences of the trade-off model, let us first discuss about the choice of raw changes introduced for enzyme activity (δ) and target affinity (γ) in Eqs.

Fig. 1

 1 Fig.1shows some representative runs for different number of copies of each

  true after correcting for the fitness of the parental cells. The same qualitative arguments apply for different trade-off functions ( ) ( )

Figure 1 .

 1 Figure 1. Sample simulations (upper panel 6 i T n = , lower panel 20 i T n = , at 0 t ) showing

Figure 2 .

 2 Figure 2. Cell growth and mutation load (dotted lines, right y-axis scale) as a function of

Figure 3 .

 3 Figure 3. Same as Fig. 1 with 100 i T n = copies for both 1 T and 2 T . The main point here is

Figure 4 .

 4 Figure 4. Sample simulations (upper panel 6 i T n = , lower panel 20 i T n = , at 0 t ) showing

Figure A1 .

 A1 Figure A1. Metabolic vesicle efficiencies (according to Eq. (2)) with (a) equal, and (b)

Figure

  Figure A2. Total growth of templates 1 T and 2 T . (a) At various template type

  Fig 2

Acknowledgments

We thank two anonymous reviewers for very helpful constructive criticisms on an earlier draft. E.S. is supported by the National Office for Research and Technology (NAP 2005/KCKHA005) and by the National Scientific Research Fund (OTKA 73047). M.S. is supported by grants CGL2009-12912-C03-01 and CGL2010-15395 from the Ministerio de Ciencia e Innovación, and 2009SGR 636 from Generalitat de Catalunya to the GBE. E.S. and M.S. are supported by funding of the INCORE project under the Sixth Research Framework Programme of the European Union. The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement n° 225167 (project e-Flux). Support by the COST D27 action (prebiotic chemistry and early evolution) and COST CM0703 (systems chemistry) are also gratefully acknowledged. The views expressed in this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission.

Appendix

Here we analyze the growth of protocells understood as the growth of total template concentration. We keep the volume of the vesicles fixed; this is an obvious simplification since volume would grow as a result of membrane growth in an autonomous protocell (cf. [START_REF] Gánti | The Principles of Life[END_REF]. Keeping cell volume fixed helps us to focus on concentrations directly.

where the growth rate of total template number is highest. It is clear that this total template growth favours some nasty parameter combinations For example, when there are few copies of template 1 T ( 1 T n is low), then total template growth is maximized when the template affinity of the other template is high (whereas that of the other template does not matter). This means that if a template has already grown to large redundancies, its growth


Vlastimil Křivan

Evolutionary stability of optimal foraging: partial preferences in the diet and patch models

Keywords: evolutionarily stable strategy, game theory, ideal free distribution, population growth, predator-prey models

In this article the patch and diet choice models of the optimal foraging theory are re-analyzed with respect to evolutionary stability of the optimal foraging strategies. In their original setting these fundamental models consider a single consumer only and the resulting fitness functions are both frequency and density independent. Such fitness function do not allow us to apply the classical game theoretical methods to study evolutionary stability of optimal foraging strategies for competing animals. In this article frequency and density dependent fitness functions of optimal foraging are derived by separation of time scales in an underlying population dynamical model and corresponding evolutionarily stable strategies are calculated. Contrary to the classical foraging models the results of the present article predict that partial preferences occur in optimal foraging strategies as a consequence of the ecological feedback of consumer preferences on consumer fitness. In the case of the patch occupation model these partial preferences correspond to the Ideal Free Distribution concept while in the case of the diet choice model they correspond to the partial inclusion of the less profitable prey type in predators diet.

Introduction

Optimal foraging theory (OFT; e. g., [START_REF] Macarthur | On optimal use of a patchy environment[END_REF][START_REF] Charnov | Optimal foraging: attack strategy of a mantid[END_REF][START_REF] Stephens | Foraging theory[END_REF] assumes that organisms forage in such a way as to maximize their fitness measured as energy intake rate. These models assume a homogeneous (fine-grained) environment with several resource types that a consumer encounters sequentially, and predict the optimal consumer diet. This line of research led to the prey model (also called the "diet choice"; [START_REF] Charnov | Optimal foraging: attack strategy of a mantid[END_REF]. This model assumes that consumer fitness is measured as the average energy intake rate. The classical example Email address: vlastimil.krivan@gmail.com (Vlastimil Křivan) Preprint submitted to Elsevier September 8, 2010

of such a situation is the experiment with great tits where a single animal feeds on two food types delivered on a conveyor belt (e. g., [START_REF] Krebs | Optimal prey selection in the great tit (Parus major )[END_REF][START_REF] Berec | Are great tits (Parus major) really optimal foragers?[END_REF] which assures that prey are not depleted by predation. Thus, no interference or competition for food between consumers is considered. A similar reasoning can be applied to the case where a consumer forages in a heterogeneous environment consisting of discrete foraging patches. In this case the consumer should move to the patch with the highest payoff. These predictions are based on fitness functions that are both density and frequency independent. They consider a single consumer and resource densities that are not influenced by predation. Thus, they miss completely competition between consumers that can be either due to interference, or due to exploitation of resources. In fact, such models can be applied only to laboratory experiments with a single consumer that does not influence resource density (e. g., a single great tit feeding on two prey types delivered on the conveyor belt).

More realistic models must consider consumer competition (either interference competition, or exploitative competition). A density dependent approach in patchy environments led to the concept of the Ideal Free Distribution (IFD; [START_REF] Fretwell | On territorial behavior and other factors influencing habitat distribution in birds[END_REF]. These authors assumed that patch payoff decreases with increasing number of individuals in that patch (e. g., this captures interference among birds) and predicted that under the IFD no individual can increase its fitness by changing its strategy measured by residence times in different patches. Because under the IFD all patches provide animals with the same fitness the question arises whether such a situation is stable or not. Indeed, assuming a single mutant with an infinitely large resident population that is distributed according to the IFD, fitness of that mutant is independent and the same as is the fitness of resident animals. One approach to study stability in such situations is to use the concept of evolutionarily stable strategies [START_REF] Smith | The logic of animal conflict[END_REF][START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF]. In fact it has been shown that the IFD is an evolutionarily stable strategy of the habitat selection game [START_REF] Křivan | The Ideal Free Distribution: A review and synthesis of the game theoretic perspective[END_REF].

However, this concept of stability assumes that the fitness function is frequency dependent, i. e., it depends both on the mutant as well as on the resident strategy. Frequency dependence allows us to predict whether a mutant strategy will increase in the resident population or not. Here the idea is that the resident strategy sets the environment (e. g., the resident consumer strategy sets the resource densities) and if a mutant strategy with a higher fitness exists then this strategy spreads, replaces residents, sets a new environment and so on until an evolutionarily stable strategy (ESS) is reached.

Another approach that considers explicitly resource-consumer population dynamics together with optimal foraging was also used in theoretical ecology (reviewed in [START_REF] Abrams | Implications of flexible foraging for interspecific interactions: lessons from simple models[END_REF]. These models show how resident individuals that follow optimal foraging strategy influence resource levels which, in turn, influence consumer numbers and strategy. Although dynamical models capture the ecological feedback, they do not provide us with information about the evolutionary stability of optimal foraging strategies. To study this question it is necessary to study conditions under which the optimal foraging strategy is resistant to mutant invasions. Studying evolutionary stability in dynamic models of optimal foraging is more complicated than in static models of behavioral ecology (such as the Hawk-Dove model, or the Prisoners dilemma;

Maynard [START_REF] Smith | The logic of animal conflict[END_REF][START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF] because it requires to analyze population dynamics for both resources, resident and mutant consumers (e. g., [START_REF] Křivan | On evolutionary stability in prey-predator models with fast behavioral dynamics[END_REF][START_REF] Cressman | The ideal free distribution as an evolutionarily stable state in density-dependent population games[END_REF].

To sum up, the problem is that some models of optimal foraging are described by fitness functions that are both frequency and density independent, and they do not allow us to predict optimal foraging strategies when interference or resource depletion occur.

On the other hand, in the full dynamic setting that describes explicitly the effect of resident strategies on resource densities, the resulting models are multidimensional and difficult to analyze with respect to evolutionary stability of optimal strategies. My aim is to show how frequency dependent fitness functions can be obtained from frequency independent fitness functions, using a time scale argument. I start with a frequency independent fitness function defined as the per capita population growth rate. Then I consider resource-consumer population dynamics and I assume that resource dynamics run on a faster time scale when compared with consumer population dynamics. This is the case of many resource-consumer systems where resource generation time is shorter than consumer generation time (e. g., many plants have short generation times relative to mammalian herbivores). This time scale separation allows me to assume that at each consumer density resources are at the corresponding population equilibrium.

Substituting this resource equilibrium to the consumer fitness function leads to a frequency dependent fitness function that can be analyzed from evolutionary perspective.

I apply this approach to diet and patch choice paradigms of optimal foraging theory.

In both cases this approach predicts emergence of partial preferences for alternative prey/patch types. In the case of the patch model partial preferences describe the consumer Ideal Free Distribution [START_REF] Fretwell | On territorial behavior and other factors influencing habitat distribution in birds[END_REF]. Although partial preferences are not predicted by the classical frequency independent models of optimal foraging, they were observed in many foraging studies (for a review see [START_REF] Stephens | Foraging theory[END_REF]. Explanations for partial preferences range from incorrect classification of resource and sampling by consumers [START_REF] Krebs | Optimal prey selection in the great tit (Parus major )[END_REF][START_REF] Rechten | Optimal prey selection: why do great tits show partial preferences?[END_REF], resource crypsis [START_REF] Erichsen | Optimal foraging and cryptic prey[END_REF], incorrect estimation of encounter rates with resource by consumers [START_REF] Mcnamara | Partial preferences and foraging[END_REF][START_REF] Hirvonen | Significance of memory properties in prey choice decisions[END_REF], limited memory capacity of consumers [START_REF] Mangel | Dynamic information and host acceptance by a tephritid fruit fly[END_REF][START_REF] Bélisle | The effects of a limited memory capacity on foraging behavior[END_REF], and limited knowledge of the environment [START_REF] Berec | A mechanistic model for partial preferences[END_REF]. In this article I will show another mechanism that explains emergence of partial preferences by the ecological feedback of consumer preferences on their fitness.

The patch occupation model

I will consider an environment consisting of two foraging patches containing resources with abundance x i in patch i (i = 1, 2). If e i denotes the efficiency with which resources are converted to new consumers, λ i is the consumer search rate for resources in patch i, and m i is the consumer mortality rate in patch i then consumer fitness expressed as the instantaneous per capita population growth rate (cf. consumer per capita population growth rate in model (3) below) is

W = (e 1 λ 1 x 1 -m 1 )v 1 + (e 2 λ 2 x 2 -m 2 )v 2 (1)
where v i (v 1 + v 2 = 1) is the proportion of the time a consumer stays in patch i. Thus, v i measures the consumer preference for patch i. Assuming that consumer preferences are adaptive, the optimal consumer strategy is to stay in patch 1 when patch payoff is higher there than in patch 2 and vice versa. Thus, the strategy that maximizes consumer fitness W is

v 1 = 1 if e 1 λ 1 x 1 -m 1 > e 2 λ 2 x 2 -m 2 0 if e 1 λ 1 x 1 -m 1 < e 2 λ 2 x 2 -m 2 . ( 2
)
This simple prediction states that consumers should move to the patch that provides them with the highest fitness. It is clear that mutants using a different strategy obtain a lower fitness and therefore, the above strategy is evolutionarily stable. However, due to interference or resource depletion it is also obvious that as more and more consumers move to a patch, this patch payoff must decrease, which is not captured by maximization of fitness function W. To model effects of consumers on resources, I will consider the following population dynamics (e. g., Fryxell andLundberg, 1994, 1997;[START_REF] Křivan | Dynamic ideal free distribution: Effects of optimal patch choice on predator-prey dynamics[END_REF][START_REF] Křivan | Adaptive foraging and flexible food web topology[END_REF])

dx 1 dt = a 1 x 1 1 - x 1 L 1 -v 1 λ 1 x 1 y dx 2 dt = a 2 x 2 1 - x 2 L 2 -v 2 λ 2 x 2 y (3) dy dt = (e 1 λ 1 x 1 -m 1 )v 1 y + (e 2 λ 2 x 2 -m 2 )v 2 y
where y is consumer density, L i is the resource i carrying capacity, and a i is the resource per capita population growth rate. I stress here that consumer fitness (1) is indeed the per capita consumer population growth rate in model (3).

Using model (3) I derive a new frequency dependent fitness function. I consider a mutant with a strategy ṽ = (ṽ 1 , ṽ2 ) in a population of residents with a strategy v = (v 1 , v 2 ) and I re-write frequency independent fitness function (1) as a function of both mutant and resident strategies

W (ṽ, v) = (e 1 λ 1 x 1 -m 1 )ṽ 1 + (e 2 λ 2 x 2 -m 2 )ṽ 2 .
The dependence of this function on the resident strategy

(v = (v 1 , v 2 )) is indirect
through resource abundance that is set by the resident strategy in model (3). To make this dependence explicit I will assume that resource population dynamics run on a fast time scale when compared with consumer population dynamics so that for each consumer density resources relatively quickly reach their corresponding equilibrium densities. When residents use strategy v = (v 1 , v 2 ), the resource equilibrium at consumer density y, calculated from (3), is

x i = L i (1 -λ i v i y/a i ), (i = 1, 2). Substituting
this density in W leads to a frequency dependent fitness function

W (ṽ, v) = r 1 ṽ1 1 - v 1 y K 1 + r 2 ṽ2 1 - v 2 y K 2 , (4)
where

r i = e i λ i L i -m i and K i = a i (e i λ i L i -m i ) e i L i λ 2 i (5)
are the per capita consumer population growth rate and the environmental carrying capacity in patch i. It is an interesting observation that the new frequency dependent fitness function is described by logistic consumer population growth in each patch. In other words, if we assume that consumer population growth in each patch is logistic, we would arrive immediately to fitness function (4). Evolutionarily stable strategy for fitness function W was derived by [START_REF] Cressman | Ideal free distributions, evolutionary games, and population dynamics in multiple-species environments[END_REF] and I briefly recall main results

here. Assuming that r 1 > r 2 , there are two possibilities. Either consumer density is low (y < (r 1r 2 )K 1 /r 1 ) in which case only the better patch 1 will be occupied, or when consumer population is above the critical threshold then both patches will be occupied by consumers, patch payoffs will be the same and the corresponding ESS is

v * 1 = K 1 r 2 K 1 r 2 + K 2 r 1 + K 1 K 2 (r 1 -r 2 ) (K 1 r 2 + K 2 r 1 )y . (6)
Under this strategy no individual can increase its fitness by changing its strategy and this strategy cannot be invaded by individuals using a different strategy. This shows that when patch payoffs are the same the ESS depends on the consumer population density. This ESS corresponds to the IFD of [START_REF] Fretwell | On territorial behavior and other factors influencing habitat distribution in birds[END_REF]. When consumers re-distribute instantaneously according to the IFD at each consumer density, consumer population growth is described by a piecewise logistic equation

dy dt = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ r 1 y 1 - y K 1 if y ≤ (r 1 -r 2 )K 1 r 1 r 1 r 2 (K 1 + K 2 ) K 2 r 1 + K 1 r 2 y 1 - y K 1 + K 2 if y > (r 1 -r 2 )K 1 r 1 , (7)
( [START_REF] Křivan | Habitat selection by two competing species in a two-habitat environment[END_REF]. The equilibrium of this model is K 1 + K 2 and at this equilibrium consumer preferences satisfy

v * 1 v * 2 = K 1 K 2 . ( 8
)
These preferences correspond to the so called balanced dispersal under which no mutants can invade [START_REF] Mcpeek | The evolution of dispersal in spatially and temporally varying environments[END_REF][START_REF] Padrón | Environmentally induced dispersal under heterogeneous logistic growth[END_REF][START_REF] Deangelis | Feedback effects between the food chain and induced defense strategies[END_REF]. Now I compare solutions of the resource-consumer model (3) (Figure 1, solid lines)

with optimal consumer strategy given by ( 2), with solutions of the single-species logistic model (7) (Figure 1, dashed lines). [START_REF] Křivan | Adaptive foraging and flexible food web topology[END_REF] showed that for L i > m i /(e i λ i ), consumer population dynamics described by (3) converge on the equilibrium density K 1 + K 2 at which the corresponding IFD is given by formula ( 8) exactly as for model ( 7). This means that model ( 7) with optimal strategy (6) converges on the same equilibrium as the original resource-consumer model (3) (Figure 1). In addition, if resource population dynamics are much faster than consumer population dynamics (Figure 1, left panel), population trajectories of the reduced model ( 7) are close to trajectories of the original model (3). If both resource and consumer population dynamics operate on comparable time scales (Figure 1, right panel) the differences between trajectories get larger.

The diet choice model

Now I will consider the diet choice model. Assuming two prey types, consumer fitness in a fine-grained environment is [START_REF] Charnov | Optimal foraging: attack strategy of a mantid[END_REF]. Here x i is density of resource i in the environment, λ i is the consumer search rate for prey type i, u i is the consumer preference for prey type i (0 ≤ u i ≤ 1), h i is the handling time, and e i is the net energy gain obtained from prey type i. I stress here that on the contrary to the patch model, the diet choice model does not assume that the sum of strategies equals to one. Optimal foraging theory assumes that resource densities are fixed, and maximization of (9), with respect to consumer preferences for either resource, leads to the zero-one rule [START_REF] Charnov | Optimal foraging: attack strategy of a mantid[END_REF]. This rule predicts that consumers will always forage on the most profitable resource type, i. e., the resource type with the highest ratio e i /h i . In what follows I will assume that the resource type 1 is more profitable than the alternative resource type 2 (i. e., e 1 /h 1 > e 2 /h 2 ).

e 1 λ 1 u 1 x 1 + e 2 λ 2 u 2 x 2 1 + h 1 λ 1 u 1 x 1 + h 2 λ 2 u 2 x 2 , (9)
Therefore, resource 1 will be always included in the consumer diet (u 1 = 1) while the second resource type will be included only if density of the first resource decreases below the threshold density

x * 1 = e 2 λ 1 (e 1 h 2 -e 2 h 1 ) , (10)
(e. g., [START_REF] Charnov | Optimal foraging: attack strategy of a mantid[END_REF][START_REF] Stephens | Foraging theory[END_REF]). Consumer strategy is not uniquely predicted by the diet choice model when the more profitable prey density equals x * 1 , because under this condition any strategy will lead to the same fitness. Now I will derive a frequency dependent fitness function and find the corresponding ESS. For this I consider ecological feedback between consumers and resources that is described by the following Rosenzweig-MacArthur predator-prey dynamics [START_REF] Křivan | Adaptive foraging and flexible food web topology[END_REF])

dx 1 dt = a 1 x 1 1 - x 1 L 1 - λ 1 x 1 y 1 + h 1 λ 1 x 1 + h 2 λ 2 u 2 x 2 dy dt = e 1 λ 1 x 1 + e 2 λ 2 u 2 x 2 1 + h 1 λ 1 x 1 + h 2 λ 2 u 2 x 2 -m y. (11)
This model considers population dynamics of the preferred resource (x 1 ) and consumers (y). It assumes that the alternative resource type density x 2 stays constant. This can be the case where there is an influx of alternative food resources to the system from outside sources (allochthonous resource; [START_REF] Huxel | Food web stability: The influence of trophic flows across habitats[END_REF]. Parameter u 2 describes consumer preferences for the alternative resource.

Fitness of a mutant consumer with strategy ũ2 in a population of residents with

strategy u 2 is W (ũ 2 , u 2 ) = e 1 λ 1 x 1 + e 2 λ 2 ũ2 x 2 1 + h 1 λ 1 x 1 + h 2 λ 2 ũ2 x 2 (12)
and I derive a frequency dependent fitness function as in the case of the patch model.

For each fixed consumer strategy the resource 1 isocline (given by dx 1 /dt = 0) is an upside down parabola (Figure 2). Assuming that resource dynamics run on a fast time scale when compared with consumer demography, trajectories move in the direction of the arrows in Figure 2 and they quickly reach the stable part of the resource isocline (solid line)

x 1+ = -1 + h 1 L 1 λ 1 -h 2 u 2 x 2 λ 2 + H(u 2 ) 1 -c(u 2 )y 2h 1 λ 1 ( 13
)
where

H(u 2 ) = 1 + h 1 L 1 λ 1 + h 2 x 2 u 2 λ 2 , c(u 2 ) = 4h 1 L 1 λ 2 1 a 1 H 2 (u 2 ) .
I remark that the resource isocline is defined only for consumer densities that satisfy y < 1/c(u 2 ) which I will assume below. Substituting the equilibrium resource density

x 1+ in (12) leads to the frequency dependent fitness function

W (ũ 2 , u 2 ) = 2h 1 (e 1 L 1 λ 1 + e 2 x 2 ũ2 λ 2 ) -e 1 H(u 2 )(1 -1 -c(u 2 )y) 2h 1 H(ũ 2 ) -h 1 H(u 2 )(1 -1 -c(u 2 )y) . ( 14
)
The selection gradient ∂W ∂ ũ2 (ũ 2 , u 2 ) predicts the direction of selective pressure [START_REF] Vincent | Evolutionary game theory, natural selection and darwinian dynamics[END_REF][START_REF] Dercole | Analysis of evolutionary processes[END_REF]. When the gradient is positive, preference for the alternative prey type will increase, while if it is negative, preference will decrease.

The gradient is zero at the singular strategy (Appendix A)

u * 2 = L 1 λ 2 1 (e 1 h 2 -e 2 h 1 ) a 1 h 2 λ 2 x 2 (L 1 λ 1 (e 1 h 2 -e 2 h 1 ) -e 2 )
y -e 1 λ 2 x 2 (e 1 h 2e 2 h 1 ) (15) provided this value is between zero and one, i. e., consumer density satisfies y * 1 ≤ y ≤ y * 2 where

y * 1 = a 1 e 1 h 2 (L 1 λ 1 (e 1 h 2 -e 2 h 1 ) -e 2 ) (e 2 h 1 -e 1 h 2 ) 2 L 1 λ 2 1 and y * 2 = a 1 h 2 (L 1 λ 1 (e 1 h 2 -e 2 h 1 ) -e 2 )(e 1 + (e 1 h 2 -e 2 h 1 )x 2 λ 2 ) (e 2 h 1 -e 1 h 2 ) 2 L 1 λ 2 1 from Appendix A.
For consumer densities for which y < y * 1 , I set u * 2 = 0 and for y > y * 2 , I set u * 2 = 1. Figure 3A shows dependence of the singular strategy on consumer density.

It is proved in Appendix A that the singular strategy is not only evolutionarily stable but also convergence stable [START_REF] Eshel | Evolutionary and continuous stability[END_REF][START_REF] Eshel | Continuous stability and evolutionary convergence[END_REF]. This means that (1) the singular strategy is resistant to invasion of a small number of mutants that use a different strategy, and (2) it is also resistant to changes in the resident strategy (i. e., when the resident strategy is slightly perturbed, it will return to the convergence stable singular point). Strategies that are evolutionarily and convergence stable are thought as the end point of evolution [START_REF] Vincent | Evolutionary game theory, natural selection and darwinian dynamics[END_REF][START_REF] Dercole | Analysis of evolutionary processes[END_REF]. Now I will analyze the singular consumer strategy as a function of consumer density.

If resource 1 carrying capacity is low so that

L 1 < x * 1 (16)
(where x * 1 is given by ( 10)) consumers will include the alternative prey type in their diet (u * 2 = 1) at all consumer densities, because the critical thresholds y * 1 and y * 2 are negative. Because x * 1 is the critical threshold predicted by the optimal foraging theory, below which the alternative prey type is included in the consumers' diet, the present model predicts that at low environmental carrying capacities for the preferred prey type consumers will behave as generalists. This is because at such low carrying capacities the density of the more profitable prey type can never be higher than is the switching threshold x * 1 .

When the carrying capacity of the preferred resource is higher than is the switching density (i. e., the inequality in ( 16) is reversed), u * 2 is a piecewise linear function of the consumer density y (Figure 3A). In this case at low consumer densities consumers will behave as specialists feeding on the more profitable resource type only. As the number of consumers increases resource 1 density will decrease and consumers will start to feed on the alternative resource as well, and their preferences for this resource will increase with consumer numbers. When at high numbers consumers will feed on the alternative food items upon each encounter. This shows that when the consumer fitness reflects frequency and density dependence given by the ecological feedback, partial preferences do arise in the diet choice model. Now I will consider consumer population dynamics. Consumer dynamics along the stable branch of the prey isocline are obtained by substituting expression for x 1+ in the consumer equation in ( 11). This leads to the following consumer population growth equation

dy dt = ry ⎛ ⎝ 1 - K(u 2 ) 1 + 1 -c(u 2 )y ⎞ ⎠ ( 17
)
where r = e 1 h 1m and K(u 2 ) = 2 e 1 +(e 1 h 2 -e 2 h 1 )x 2 λ 2 u 2 H(u 2 )(e 1 -h 1 m) (Appendix B). Substituting for u 2 in model ( 17) the singular strategy (u * 2 ), describes consumer dynamics with optimal diet selection. For consumer densities for which the consumer singular strategy u * 2 is between 0 and 1, consumer population dynamics simplify to

dy dt = y e 2 -h 2 m h 2 .
Depending on the sign of e 2h 2 m, consumer density either increases or decreases when partial preferences occur. Because the above equation has no non-trivial population equilibrium it follows that no population equilibrium where consumers would show partial preferences for the alternative prey type exists. Consumer equilibria for consumer densities that are higher than y * 2 (i. e., when u * 2 = 1) or smaller than y * 1 (i. e., when u * 2 = 0) are given by the interior equilibrium of ( 17)

y * = K(u 2 )(2-K(u 2 )) c(u 2 )
where I substitute 0 or 1 for u 2 (these are given explicitly in Appendix B). Figure 3B shows one such trajectory of model ( 17) driven by the consumer singular strategy (solid lines). The corresponding resource 1 density is calculated from (13). Figure 3C shows the corresponding consumer singular strategy. These results predict that at low consumer densities where resource 1 density is near to its carrying capacity, consumers specialize on the more profitable resource type. As consumers increase in numbers, the preferred resource density decreases to the critical switching density (x * 1 ) predicted by the classical model of optimal foraging. From then on, consumer preferences for the alternative prey type start to increase, keeping the preferred resource density at the switching threshold for some time (i. e., times approx. between 5 and 14 in Figure 3B).

Because e 2mh 2 > 0 in Figure 3B, consumer density increases (it would decrease otherwise). Once the consumer preference reaches 1 (i. e., consumers attack the alternative resource upon each encounter), consumers and resources tend to a stable equilibrium.

This general pattern of population dynamics is also clearly shown in Figure 3D where the trajectory from panel B is projected on the consumer isocline (the solid dot denotes the population equilibrium). Thus, in the region of consumer densities where preferences for the less profitable resource type are between 0 and 1, consumers exhibit partial preferences. This is a new result because the classical diet choice model does not predict such a gradual shift in consumer preferences at the switching prey density [START_REF] Stephens | Foraging theory[END_REF]. Trajectories of model ( 11) with the optimal consumer foraging strategy are shown in Figure 3B as dashed lines. It is clear that they converge to the same equilibrium as trajectories of the single species model. Using a completely different approach [START_REF] Křivan | Optimal foraging and predator-prey dynamics[END_REF] (see formula (32) there) calculated consumer strategy at the switching threshold x 1 = x * 1 (i. e., when the classical diet choice model does not define consumer preference for the alternative prey type uniquely). It is an interesting observation that the formula for partial preferences obtained there coincides with the singular strategy u * 2 .

Discussion

In this article I have shown how to derive a frequency dependent fitness function from a frequency independent fitness function by using a time-scale argument. The time scale argument assumes that resource population dynamics run on a faster time scale when compared with consumer population dynamics. Having a frequency dependent fitness function allowed me to study optimal foraging strategies for multiple consumers and depletable resources. In particular I showed that partial preferences for alternative resources arise in the diet and patch models of optimal foraging.

The diet and patch choice models are two paradigms of the optimal foraging theory [START_REF] Macarthur | On optimal use of a patchy environment[END_REF][START_REF] Emlen | The role of time and energy in food preferences[END_REF]. They assume that per capita energy intake rate is a proxy for consumer fitness that is maximized by adaptive consumer foraging behavior. The diet choice model assumes that two or more resources are evenly distributed in the environment and consumption of these resources does not influence their numbers. This assumption is clearly violated in most natural systems. To comply with these assumptions, e. g., conveyor belts that bring food directly to consumers were used in experiments [START_REF] Krebs | Optimal prey selection in the great tit (Parus major )[END_REF][START_REF] Berec | Are great tits (Parus major) really optimal foragers?[END_REF]. So what is missing in the diet choice model is the link between consumer numbers and consumer strategy. This link defines an ecological feedback mechanism: consumer foraging strategy influences resource densities which, in turn, set the consumer strategy and density. To model this ecological feedback mechanism some authors considered population dynamics together with the optimal foraging models (e. g., [START_REF] Fryxell | Diet choice and predator-prey dynamics[END_REF][START_REF] Křivan | Optimal foraging and predator-prey dynamics[END_REF][START_REF] Fryxell | Individual behavior and community dynamics[END_REF][START_REF] Křivan | Dynamic ideal free distribution: Effects of optimal patch choice on predator-prey dynamics[END_REF][START_REF] Křivan | Adaptive foraging and flexible food web topology[END_REF][START_REF] Ma | Dynamic versus instantaneous models of diet choice[END_REF].

In these models it often happens that population dynamics tend to densities at which models of optimal foraging do not predict the optimal strategy uniquely. For example, in the patch model consumers will distribute between the two patches so that patch payoffs will equalize [START_REF] Křivan | Dynamic ideal free distribution: Effects of optimal patch choice on predator-prey dynamics[END_REF]. Under this situation there is no selection against mutants that can use any strategy, because individual fitness is independent of the strategy. Similarly, dynamical models of diet choice drive periodically density of the more profitable prey type to the threshold where optimal consumer diet is not uniquely predicted [START_REF] Křivan | Optimal foraging and predator-prey dynamics[END_REF]. In fact, this non-uniqueness causes fundamental problems for the existence of solutions to these population models [START_REF] Colombo | Selective strategies in food webs[END_REF][START_REF] Křivan | Optimal foraging and predator-prey dynamics[END_REF].

A simpler approach to study evolutionary stability is based on the concept of evolutionarily stable strategies (e. g., [START_REF] Smith | The logic of animal conflict[END_REF][START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF][START_REF] Cressman | Evolutionary Dynamics and Extensive Form Games[END_REF][START_REF] Vincent | Evolutionary game theory, natural selection and darwinian dynamics[END_REF]. However, to apply methods of static game theory, the fitness function must be frequency dependent, i. e., it must allow us to measure the effect of mutants on the resident strategy. In this article I derived two frequency dependent fitness functions for patch and diet choice models.

These fitness functions are based on the assumption that resource population dynamics run on a faster population time scale when compared with consumer population dynamics. Using these fitness functions I derived the corresponding evolutionarily stable strategies for the patch and the diet choice model. In both cases the corresponding ESS predicts partial preferences for resources to occur. In the case of the patch model separation of the resource and consumer time scales leads to consumer population growth that is described by the logistic equation for which the evolutionarily stable strategies were studied (e. g., [START_REF] Křivan | Habitat selection by two competing species in a two-habitat environment[END_REF][START_REF] Morris | Shadows of predation: habitat-selecting consumers eclipse competition between coexisting prey[END_REF][START_REF] Cressman | Migration dynamics for the ideal free distribution[END_REF]. These works show that the evolutionarily stable strategy corresponds with the IFD [START_REF] Fretwell | On territorial behavior and other factors influencing habitat distribution in birds[END_REF]. Similarly, when time scaling is applied to the diet choice model the resulting consumer population growth equation seems to be new.

Moreover, partial preferences arise for a range of consumer densities in both models.

This is a new prediction because partial preferences were not predicted by the optimal foraging theory [START_REF] Charnov | Optimal foraging: attack strategy of a mantid[END_REF]. In fact, my present analysis clearly shows that at low consumer densities consumers should specialize on the more profitable prey/patch type only. However, as consumer density increases, consumers also start to utilize the alternative, less profitable prey/patch type with increasing strength. This prediction should be easily tested using appropriate data on consumer preferences.

These predictions may have several consequences in population ecology. For example, if consumer preferences depend on consumer numbers, functional responses used in multiple species models should reflect this situation. Thus, they should depend not only on resource densities, but also on consumer densities (for a recent review of flexible foraging on the functional response see [START_REF] Abrams | Implications of flexible foraging for interspecific interactions: lessons from simple models[END_REF]. Such a dependency could lead to a more mechanistic explanation of the effect of consumer density on the functional response.

Appendix A: Evolutionary stability of the diet choice model

Fitness function ( 14) can be written as

W (ũ 2 , u 2 ) = Aũ 2 + B C ũ2 + D where A = 2e 2 h 1 x 2 λ 2 , B = 2h 1 e 1 L 1 λ 1 -e 1 H(u 2 )(1 -1 -c(u 2 )y), C = 2h 1 h 2 x 2 λ 2 , D = 2h 1 (1 + h 1 L 1 λ 1 ) -h 1 H(u 2 )(1 -1 -c(u 2 )y), c(u 2 ) = 4h 1 L 1 λ 2 1 a 1 (1 + h 1 L 1 λ 1 + h 2 x 2 u 2 λ 2 ) 2 , H(u 2 ) = 1 + h 1 L 1 λ 1 + h 2 x 2 u 2 λ 2 .
First, I calculate the singular strategy and study its evolutionary stability. The gradient of the fitness function with respect to the mutant strategy

∂W ∂ ũ2 (ũ 2 , u 2 ) = AD -BC (D + C ũ2 ) 2 .
At the singular strategy this gradient when evaluated at ũ2 = u 2 must be zero, i. e., AD -BC = 0. When I substitute expressions for A, B, C, D and after some simplification I obtain the following equation

1 - 4h 1 L 1 yλ 2 1 a 1 (1 + h 1 L 1 λ 1 + h 2 u 2 x 2 λ 2 ) 2 = 1 - 2h 1 ((e 1 h 2 -e 2 h 1 )L 1 λ 1 -e 2 ) (e 1 h 2 -e 2 h 1 )(1 + h 1 L 1 λ 1 + h 2 u 2 x 2 λ 2 ) . (18)
This equation posses a solution only provided the right hand side is between 0 and 1, i. e., when

e 2 λ 1 (e 1 h 2 -e 2 h 1 ) < L 1 < e 2 h 1 + e 1 h 2 + h 2 (e 1 h 2 -e 2 h 1 )u 2 x 2 λ 2 h 1 (e 1 h 2 -e 2 h 1 )λ 1 .
Solving equation ( 18) yields the singular strategy (15). Of course, u * 2 must be between 0 and 1 which holds for consumer densities that satisfy y * 1 ≤ y ≤ y * 2 where

y * 1 = a 1 e 1 h 2 (L 1 λ 1 (e 1 h 2 -e 2 h 1 ) -e 2 ) (e 2 h 1 -e 1 h 2 ) 2 L 1 λ 2 1 and y * 2 = a 1 h 2 (L 1 λ 1 (e 1 h 2 -e 2 h 1 ) -e 2 )(e 1 + (e 1 h 2 -e 2 h 1 )x 2 λ 2 ) (e 2 h 1 -e 1 h 2 ) 2 L 1 λ 2 1 .
If consumer density is too low (y < y * 1 ), consumers will feed on the more profitable prey type only while at high densities (y > y * 2 ) they will feed on the alternative prey type upon each encounter. In particular, the carrying capacity for the more profitable prey type must be high enough, i. e., L 1 > e 2 λ 1 (e 1 h 2e 2 h 1 ) is negative, i. e., 1 < K < 2. I remark that condition 1 < K is equivalent to

L 1 < e 1 + h 1 m + (-2e 2 h 1 + h 2 (e 1 + h 1 m))u 2 x 2 λ 2 h 1 (e 1 -h 1 m)λ 1
which is the condition for the interior equilibrium of the resource-consumer model (11), to be stable. For larger values of the environmental carrying capacities the interior equilibrium is unstable and a stable limit cycle arises in the Rosenzweig-MacArthur model ( 11) [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF].

The above analysis assumed fixed consumer preference for the alternative resource.

-m(1 + h 1 L 1 λ 1 + h 2 x 2 λ 2 )) L 1 (e 1 -h 1 m) 2 λ 2 1 .
Similarly when y < y * 1 , u 2 = 0 and population dynamics on the stable manifold are obtained by substituting u 2 = 0 to (17). These population dynamics have equilibrium y 0 eq = a 1 e 1 (e 1 L 1 λ 1m(1 + h 1 L 1 λ 1 )) L 1 (e 1h 1 m) 2 λ 2 1 .

Figures caption.

Figure 1: Solutions of the patch model (3) (solid line) where consumers follow the optimal foraging strategy given by (2), and solutions of the single-species model (7)(dashed line). Resource densities for model (7) are given by x i = L i (1y(t)λ i v i /a i ) with the optimal strategy v i given by ( 6). The left panel (A,B,C) assumes that con-

  sumer demographic parameters (e 1 = 0.015, e 2 = 0.01, m 1 = m 2 = 0.02) are much smaller when compared to resource parameters. This discrepancy causes resource population dynamics to run on a fast time scale when compared with consumer population dynamics. The right panel (D,E,F) assumes more similar time scales for both resource and consumer dynamics (e 1 = 0.15, e 2 = 0.1, m 1 = m 2 = 0.2). Other parameters used in simulations: L 1 = L 2 = 10, a 1 = 1.5, a 2 = 0.5, λ 1 = λ 2 = 1.

Figure 2 :

 2 Figure 2: Resource isocline of model (11). Only the solid part of the isocline is stable provided resource 1 population dynamics are fast when compared to consumer population dynamics. Arrows indicate direction of trajectories.

Figure 3 :

 3 Figure 3: Panel A shows consumer preference (15) for the alternative prey type in the diet choice model as a function of consumer density. Panel B compares a trajectory of the resource-consumer model (11) (dashed lines) where consumers follow predictions of the classical prey model (i. e., zero-one rule) with a trajectory of the reduced model (17) (solid lines). The corresponding trajectory for resources is given by (13). Panel C shows the corresponding consumer optimal strategy. Panel D shows the solid trajectory from panel B in the resource 1-consumer preference-consumer density phase space. Parameters: a 1 = 1.2, h 1 = h 2 = 0.2, e 1 = 0.2, e 2 = 0.05, λ 1 = λ 2 = 0.5, x 2 = 8, m = 0.2, L 1 = 10.

  Figure 1:

  For h 2 < e 2 /m consumer densities increase, while for larger handling times they decrease. No nontrivial equilibrium exists. If y ≥ y * 2 then u 2 = 1 and population dynamics on the stable manifold are obtained by substituting u 2 = 1 to (17). These population dynamics have equilibriumy 1 eq = a 1 (e 1 + (e 1 h 2e 2 h 1 )x 2 λ 2 )(e 1 L 1 λ 1 + e 2 x 2 λ 2

	Now I will study population dynamics (17) driven by the singular strategy. Substituting
	singular strategy (15) in population dynamics (17) I get
	dy dt	= y	(e 2 -h 2 m) h 2
	for consumer densities satisfying y * 1 < y(t) < y * 2 and (19).

Acknowledgments

The Institute of Entomology is funded by the Academy of Sciences of the Czech Republic (Z50070508). A part of this work was conducted while a Sabbatical Fellow at the National Institute for Mathematical and Biological Synthesis, an Institute sponsored by the National Science Foundation, the U.S. Department of Homeland Security, and the U.S. Department of Agriculture through NSF Award #EF-0832858, with additional support from The University of Tennessee, Knoxville.

for consumer partial preferences to arise. If the opposite inequality holds, fitness maximizes at u 2 = 1. I derived the singular strategy under the assumption that the right hand side of equation ( 18) is positive, because otherwise no singular solution exists.

Substituting the singular strategy to the right hand side of (18) leads to expression

This condition must hold together with other constraints for the singular strategy to be between 0 and 1.

Because AD -BC = 0 at the singular strategy, W (ũ 2 , u * 2 ) is independent of the mutant strategy ũ2 and, after some calculations, W (ũ 2 , u * 2 ) = e 2 /h 2 for all mutant strategies u 2 .

To study evolutionary stability of the singular strategy I need to check the local stability condition [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF] that asserts that mutants cannot spread in the resident population. The local stability condition requires that

for every u 2 = u * 2 in a neighborhood of the singular strategy

2 ) = 0 and with a help of some computer algebra package (e. g., Mathematica) it is easy to show that provided inequality (19) holds, g (u * 2 ) = 0 and

Thus, g(u 2 ) > 0 in a neighborhood of the singular strategy (u 2 = u * 2 ) and the local ESS condition (20) holds.

Second, I will prove that the singular strategy is also continuously stable. This follows from the fact [START_REF] Eshel | Evolutionary and continuous stability[END_REF] that

for consumer densities that satisfy (19).

Appendix B: Derivation of model (17).

Substituting x 1+ given by ( 13) in the right hand side of the consumer population growth equation ( 11) leads to

where

Extending the above fraction by expression

where r = e 1 h 1m, K = 2 e 1 +u 2 λ 2 x 2 (e 1 h 2 -e 2 h 1 )

There exists a non-zero equilibrium of the single species consumer model ( 21)

which is exactly the same as the interior consumer equilibrium of model ( 11). For y * to be an equilibrium, it must be positive, i. e., 0 < K < 2. This condition is equivalent to