Lotka–Volterra dynamics kills the Red Queen: population size fluctuations and associated stochasticity dramatically change host-parasite coevolution

Background Host-parasite coevolution is generally believed to follow Red Queen dynamics consisting of ongoing oscillations in the frequencies of interacting host and parasite alleles. This belief is founded on previous theoretical work, which assumes infinite or constant population size. To what extent are such sustained oscillations realistic? Results Here, we use a related mathematical modeling approach to demonstrate that ongoing Red Queen dynamics is unlikely. In fact, they collapse rapidly when two critical pieces of realism are acknowledged: (i) population size fluctuations, caused by the antagonism of the interaction in concordance with the Lotka-Volterra relationship; and (ii) stochasticity, acting in any finite population. Together, these two factors cause fast allele fixation. Fixation is not restricted to common alleles, as expected from drift, but also seen for originally rare alleles under a wide parameter space, potentially facilitating spread of novel variants. Conclusion Our results call for a paradigm shift in our understanding of host-parasite coevolution, strongly suggesting that these are driven by recurrent selective sweeps rather than continuous allele oscillations.


Background
The Red Queen from Lewis Carroll's tale 'Through the looking glass' is commonly used as a metaphor for selection-induced rapid evolution [1][2][3]. It is based on the observation that persistence in an environment with changing selective constraints requires ongoing adaptations to the encountered challenges [4]. Host-parasite coevolution with antagonistic and inter-dependent interactions represents one of the role models for such rapid evolutionary change [5,6]. For instance, an increase in host resistance reduces parasite fitness, thus immediately favoring parasite varieties with altered virulence and/or immune-evasion mechanisms. In turn, a novel parasite attack mechanism decreases host fitness, thus favoring host varieties with new counter-defenses. If the interaction persists, then it will lead to continuous parasite adaptations and host counter-adaptations. The rapid evolutionary dynamics associated with these interactions is very well documented in the literature, ranging from field studies on rabbits and their myxoma viruses [7], snails and their trematode parasites [8], Daphnia magna waterfleas and their bacterial parasites [9] to laboratorybased coevolution experiments between bacteria and their phages [10][11][12], the nematode Caenorhabditis elegans and bacterial parasites [13,14], or the red flour beetle Tribolium castaneum and its microsporidian parasite [15,16].
It is thus widely accepted that these interactions evolve fast and continuously. Yet, to date, the exact underlying selection dynamics are not always well understood. http://www.biomedcentral.com/1471-2148/ 13/254 These dynamics can generally be influenced by metapopulation structure and environmental variation [17,18]. Within a particular population and specific environmental context, two main alternatives are thought to be of prime importance: recurrent selective sweeps and negative frequency-dependent selection [5,6,[19][20][21]. Both alternatives are consistent with the above original definition of the Red Queen hypothesis by Van Valen [4], whereas, curiously, only the second alternative is referred to as Red Queen dynamics [5,6,20] . The two alternatives are closely related because both assume a selective advantage of a rare genotype, for example a novel host resistance variant. However, they differ fundamentally in the way in which the new variant originates and spreads within the population. The concept of recurrent selective sweeps (often termed arms race dynamics) consists of two steps: the de novo appearance of a beneficial allele (e.g., by mutation or immigration) and its subsequent spread through the population to fixation (i.e., the selective sweep). These sweeps occur repeatedly in host and parasite populations, usually each time with a new beneficial allele. They may only lead to fast changes in absolute time if at least one of the following factors applies: new alleles arise frequently, new alleles become immediately visible and thus selectable at the phenotypic level, the new alleles provide a high selective advantage, and/or the organisms have short generation times. This situation is best met in bacteriaphage interactions, which are usually characterized by large population sizes (i.e., high likelihood of the occurrence of favorable mutations), short generation times, and haploid genomes (i.e., new mutations are immediately expressed phenotypically) [11,[22][23][24] (but see also [25]).
In contrast, the dynamics for multicellular host systems are traditionally viewed to be determined by negative frequency-dependent selection leading to sustained oscillations of the same alleles (i.e., Red Queen dynamics [6,20]), but not to the fixation of single alleles. In this case, standing genetic variation is required, because the population sizes for these hosts are usually comparatively small, their generation times comparatively long, and their genomes diploid. As a consequence, recurrent selective sweeps are commonly thought to be rather slow in these systems. Instead, if standing genetic variation is available, then negative frequency-dependent selection can produce fast and continuous allele frequency changes even in these host systems. Such negative frequency-dependent dynamics seem to be present in some multicellular host systems, including the freshwater snail Potamopyrgus antipodarum [8,26] and the waterflea Daphnia magna [9].
Numerous theoretical models have been developed to study the underlying selection dynamics. Interestingly, the current models typically focus on evolutionary change (i.e., the rate of change in host and parasite allele frequencies in response to the type of interaction). These approaches have thus largely neglected ecological dynamics, which can have a huge impact on the evolutionary process. Population size fluctuations deserve particular attention in this context, because they are induced by reciprocal selection among the antagonists and, therefore, represent an inherent property of host-parasite coevolution -irrespective of additional environmental variation [7,10,[27][28][29][30]. Since selection is reciprocal, population size fluctuations should be coupled between the antagonists, and generally follow Lotka-Volterra dynamics [31,32]. Such demographic variations have the potential to affect the dynamics of host-parasite allele frequency changes by introducing two important effects. Firstly, the rising and falling population sizes produce bottlenecks where selection favours a particular allele. The favored allele may thus reach comparatively high frequencies during the bottleneck, possibly enhancing its spread in the subsequently expanding population. Secondly, the elevated stochasticity during the bottleneck may lead to a further increase and thus spread of the favored allele.
In this manuscript, we aim at understanding in how far Lotka-Volterra population size fluctuations and the associated stochastic effects influence the dynamics of allele frequency changes during host-parasite coevolution. While several previous theoretical models have applied the Lotka-Volterra dynamics to host-parasite coevolution (e.g., [33][34][35][36][37]), their influence on the evolutionary dynamics has not yet been systematically explored by comparison with a model with constant population size. Similarly, stochastic effects during host-parasite coevolution have only been considered in a few theoretical studies (e.g., [38,39]), yet, to our knowledge, with a single exception [40] under constant population size and not in combination with Lotka-Volterra dynamics. Hence, while the previous studies have independently utilised stochastic effects or Lotka-Volterra dynamics, a systematic analysis of the consequences of each of these factors, either alone or in combination, is as yet missing -in spite of their potential importance. The novelty of our study lies in bringing together these two aspects and comparing their influence to the traditional model, in which Lotka-Volterra dynamics and stochastic effects are excluded. More specifically, we here use the standard matching-alleles host-parasite interaction model to assess allele frequency dynamics in the presence versus absence of Lotka-Volterra oscillations for a stochastic versus an analogous deterministic model.

Methods
Based on the Lotka-Volterra equations [31,32], we address the population dynamics of interacting hosts and parasites. The host corresponds to the prey in the original model, and the parasite to the predator. The host http://www.biomedcentral.com/1471-2148/13/254 consumes a (constant) food supply F and reproduces at rate c 1 . Parasites infect hosts at rate c 2 , leading to elimination of a host and generation of an additional parasite. Parasites die at rate c 3 . The number of host and parasite individuals are given by H and P. In a stochastic system these interactions can be defined by the following reactions [41,42], Usage of these specific reactions facilitates tracking of each unit of the interacting antagonists and, thus, it allows a more precise characterization of the resulting dynamics. These reactions can also be used directly for exact stochastic simulations based on the Gillespie algorithm. They further provide a microscopic dynamics from which the deterministic Lotka-Volterra equations emerge in the limit of infinite population size [42], Host-parasite coevolution is modeled using the standard matching alleles model [6]. For this, we define two host and two parasite types, H 1 and H 2 for the host and P 1 and P 2 for the parasite. This is equivalent to a haploid system with two antagonists, each of which possesses two alleles at a single locus. The interaction according to the matching alleles model is described with the following six reactions, In the matching alleles model, the interactions between alternate hosts and parasites (H 1 , P 2 and H 2 , P 1 ) are without consequence and thus do not appear here. While the absence of these interactions is the standard assumption in the matching alleles model, allowing a small amount of these interactions does not change our results qualitatively (see Appendix). In the limit of infinite population size [42], we obtain a set of four coupled nonlinear differential equations, where the frequencies of H i and P i are given by h i and p i . The above equations consider interdependence of host and parasite demographies, allowing population sizes to vary in response to the interaction with the antagonist, consistent with the Lotka-Volterra model. The precise nature of the resulting oscillations in population size is determined by the parameters, most importantly by b.
As we are interested in the effects of population size variation induced by the Lotka-Volterra equations, we have to compare this to a scenario in which the population size is constant. Such constant population size models are common, e.g. the Wright-Fisher model or the Moran process. However, microscopically these models are distinct from the Lotka-Volterra equations considered above. Therefore, we used the above approach and enforced constant population size by resetting host and parasite population sizes to their initial values after every generation (N avg transition events, see Appendix), while relative allele frequencies were maintained. The dynamics were subsequently assessed for different average population sizes. To ensure comparability of allele frequency fluctuations across population sizes and evolutionary models, we rescaled the interaction parameters with N avg for the deterministic analogues of the considered stochastic scenarios (Appendix).

Results
Host-parasite coevolutionary dynamics are analyzed in the presence and absence of Lotka-Volterra dynamics. Figure 1 illustrates an exemplary result. All models initially produce oscillatory allele frequency changes, but only with Lotka-Volterra dynamics are these accompanied by changes of population size. As a consequence, changes in allele numbers are also more pronounced (top versus bottom in Figure 1). As the deterministic model allows for arbitrary small frequencies of each type, it formally never leads to allele fixation and thus produces continuous oscillations. In contrast, the corresponding stochastic models have absorbing states, making fixation possible. Interestingly, allele fixation appears to be substantially faster in the stochastic model that includes Lotka-Volterra fluctuations (top versus bottom panels, Figure 1). As such, it seems that these conditions favor rapid termination of the Red Queen oscillations.
We next analyze the impact of the average population size on this pattern. In the following, we focus on the time http://www.biomedcentral.com/1471-2148/13/254   , then simulations may initially produce allele oscillations as above and below. However, alleles usually spread to fixation (or go extinct) at a much faster rate. Bottom: Dynamics without Lotka-Volterra cycles, fixing the average population size of both species to N avg = 1000 by resetting it after every N avg reactions, while maintaining the ratio between the alleles. The 50 individual stochastic simulations now only rarely reach fixation. The figure illustrates the scenario where the rare host allele (H 1 ) is more likely to reach fixation than the frequent host allele (see Figure 3). This fixation probability decreases with increasing initial frequency (cf. Figure 3). The simulation parameters are a = 5, c = 2.5, b = 10/N avg = 0.01 with H 1 = 5%, H 2 = 95%, P 1 = 20%, P 2 = 80% as initial condition.
until one of the alleles from either of the antagonists has reached fixation in order to compare evolutionary rates across population sizes and models. In general, Lotka-Volterra dynamics cause a substantial increase in allele fixation rate ( Figure 2). Interestingly, in this case, allele fixation rates depend only weakly on average population size. Figure 1 suggests that this is because allele frequencies can become very small during the Lotka-Volterra demographic fluctuations. In contrast, average population size has a much stronger effect when it is artificially kept constant. Here, the time until allele fixation increases exponentially with increasing population size ( Figure 2). Figure 2 explores the time to fixation of any of the alleles in either the host or the parasite using a specific combination of initial allele frequencies (i.e., the rare host allele is present at 5%, the common at 95%, whereas the parasite alleles are at 20% and 80% respectively). How does this depend on the initial allele frequencies in both antagonists? For instance, the selective advantage of a rare allele is not only the result of its own frequency, but also determined by the abundance of the corresponding allele in the antagonist. Allele fixation rates were thus explored as a function of initial allele frequencies in the  . This is true across a relatively wide distribution of initial frequencies for the corresponding parasite allele. Interestingly, it even applies when the corresponding parasite allele has high initial frequencies (Figure 3, top left corner in top middle panel). This counterintuitive result can be explained by consideration of the dynamics that ensue from these initial conditions. In this particular case, the low initial frequency of host allele 1 means that host allele 2 is initially common, whereas the high initial frequency of parasite allele 1 means that parasite allele 2 is rare. High host allele 2 abundance then specifically favors parasite allele 2, which rapidly increases in frequency. The unexpected consequence of these starting conditions is that these two interacting alleles subsequently engage in highly pronounced frequency oscillations that show larger amplitudes than those observed for host and parasite alleles 1 (Figure 4). If during these oscillations low allele 2 frequencies coincide with a Lotka-Volterra bottleneck and associated stochastic effects, then host allele 2 has a very high likelihood of going extinct, resulting in fixation of host allele 1 (see Figure 4).
The results also highlight that the dynamics are usually determined by fixation of one of the host alleles (red colour is mainly found in middle rather than right panel of the top row of Figure 3). Note that the simulations are stopped as soon as either one of the host or one of the parasite alleles reaches fixation and, thus, the fixation probabilities of both host as well as both parasite alleles sum up to one. In our case, fixation of the host allele is more likely than fixation of the parasite allele because for our parameter combination and initial condition, it is usually the host that first experiences a Lotka-Volterra bottleneck and consequentially a drop in the frequency of one of the alleles (see also Figure 4). Nevertheless, if both the parasite and corresponding host allele are common, then it is the parasite allele that has a high probability of fixation ( Figure 3, top right).
The overall pattern looks different in the absence of Lotka-Volterra fluctuations (Figure 3, bottom row). Host allele fixation probability increases with its own high initial frequency and, at the same time, low initial abundance of the corresponding parasite allele (Figure 3 middle panel in bottom row). Parasite allele fixation is enhanced when both parasite and corresponding host alleles are initially common (Figure 3 bottom right). Under these conditions, fixation probabilities of both host and parasite alleles are almost identical at initially intermediate frequencies, most likely due to negative frequency dependent selection, as illustrated in Figure 1.      In this figure we depict the dynamics that occur at the initial conditions with a low H 1 frequency and a high P 1 frequency. In this particular case, the low initial frequency of H 1 means that H 2 is initially common, which in turn favours P 2 . This parasite allele thus rapidly increases in frequency, subsequently causing highly pronounced H 2 and P 2 frequency oscillations that show larger amplitudes than the interacting H 1 and P 1 alleles. If low H 2 frequencies coincide with a Lotka-Volterra bottleneck in the hosts, then the associated stochastic effects lead to a higher likelihood of H 2 going extinct, resulting in an overall higher fixation probability of H 1 . The top panel shows the average population dynamics, whereas the bottom panel shows the frequency changes for the indicated host and parasite alleles across the ten independent simulations. The vertical lines in the bottom panel denote the time points where the simulation is terminated due to a loss of an allele. Out of the 10 simulation runs 9 are stopped due to the allele H 2 going extinct and only one due to H 1 going extinct. The interaction parameters are a = 5, c = 2.5, b = 10/N avg = 0.01.

Discussion
Population size fluctuations represent an inherent property of host-parasite interactions. Unequivocal evidence for such interaction-dependent demographic variations was obtained from controlled host-parasite mesocosm experiments, for example with E. coli and its phage [10,30], Hydra hosts and its Hydramoeba parasite [43], house fly and its parasitic wasps [44,45], or azuki bean weevil Callosobruchus chinensis and its parasitoid Heteropilus prosopidis [46]. Similar observations were made under field conditions, for example for rabbits and their myxoma viruses [7], or red grouse and its nematode parasite [47]. Additional examples are summarized by [28] and [29]. As population size fluctuations produce regular bottlenecks, random genetic drift is likely to influence allele frequencies. Previous theoretical models, developed in a different context, strongly suggest that even large populations are influenced by such stochastic processes [48,49]. More generally, under natural conditions in a finite population, it is difficult to imagine that changes in population size do not affect evolutionary dynamics.
Consequently, an in-depth understanding of the evolution of host-parasite interactions should take account of the associated ecological processes based on Lotka-Volterra fluctuations.
Very few previous theoretical models on host-parasite interactions have considered Lotka-Volterra fluctuations [33][34][35][36][37]. These studies usually used a deterministic approach and thus excluded stochastic effects, which are most prominent during bottlenecks. Similarly, only few theoretical studies considered stochastic effects in this context [38,39], yet under constant population size, but not in combination with Lotka-Volterra dynamics. We are aware of only one study that looked at host-parasite coevolution in consideration of Lotka-Volterra interactions and stochasticity [40]. However, this study had a different focus, and thus, it did not include a systematic comparison to models without Lotka-Volterra cycles or without stochasticity. Consequently, the interaction of these two aspects for host-parasite coevolution is so far unexplored. At the same time, their relevance was demonstrated for evolution of only one of the antagonists, namely the parasite. For example, the probability of fixation of a beneficial mutation in a bacterial population was shown to be enhanced by periodical bottlenecks [50][51][52]. Similar results were obtained in a model that explored the effect of bottlenecks during pathogen transmission [53]. Our study explicitly evaluates the influence of both Lotka-Volterra fluctuations and stochastic effects on the dynamics of host-parasite interactions using a comparison to a model with constant population size and/or absence of stochasticity. As the demographic variations are an inherent property of such antagonistic interactions, their influence should apply across a wide range of environmental http://www.biomedcentral.com/1471-2148/13/254 conditions and thus be of general relevance for our understanding of host-parasite coevolution.
Based on our approach, we obtained evidence that both Lotka-Volterra fluctuations and associated stochastic effects significantly affect the course and pace of coevolutionary adaptations. In particular, both factors facilitate selective sweeps (i.e., the spread and fixation of an allele). Most impressively, this effect appears to be independent of average population size ( Figure 2) and occur at a substantially faster rate (Figure 3, left column). Moreover, allele frequency changes are not exclusively due to drift, which should favor fixation of initially frequent alleles and loss of initially rare alleles. In contrast, our results indicate that initially rare host alleles can spread to fixation across a relatively wide range of conditions ( Figure 3, top middle panel). Rare parasite alleles may not necessarily go extinct, but have a certain likelihood of spreading contingent on the frequency of the corresponding host allele (Figure 3, top right panel). Based on these results we propose that selective sweeps rather than oscillatory negative frequency-dependent selection may represent the main driving force during host-parasite coevolution.
Recurrent selective sweeps have been repeatedly suggested to determine coevolutionary dynamics for parasite or host systems with large population sizes such as bacterial hosts or microbial parasites, where novel mutations are frequent and often directly exposed to selection because of a haploid genetic system. If these selective dynamics also apply to multicellular host and parasite systems, then two contrasting effects may be expected on the coevolutionary process. On the one hand, these systems usually have much smaller population sizes, facilitating spread of alleles in spite of the often diploid genetic system. On the other hand, continuous coevolution may become difficult because it is usually assumed that small population size results in a reduced likelihood of the occurrence of advantageous novel mutations [6]. However, the latter assumption may not always be true. It is possible that new alleles become available for example by frequent immigration or a high rate of gene duplication. These processes may further favor the formation of novel genotypes if they act in combination with recombination and/or mutation. Interestingly, the possible impact of gene duplications is usually not addressed in theoretical work on host-parasite coevolution, even though such duplications are known to be common in almost all organisms [54][55][56] and often affect genes of relevance for the interaction such as virulence genes in parasites [57,58] or immunity genes in animal hosts [54,59,60].
Several additional factors may favor ongoing coevolution. One of these is founded on a more complex genetic architecture underlying host-parasite co-adaptation that consists of several interacting loci across the genome (e.g., [61][62][63][64]). In this case, allele fixation at one locus may still permit maintenance of variation at other loci, which could then mediate evolutionary responses to the antagonist. Yet another possibility may depend on metapopulation structure, consisting of coevolutionary hot and cold spots and migration among demes, as evidenced for flax and its parasitic rust fungus [65] or the above mentioned snail-trematode interaction [17]. Such an interconnected network could then maintain allelic diversity across the entire metapopulation, even if alleles become temporarily fixed within single demes. Moreover, environmental gradients or perturbations are known to influence hostparasite coevolutionary dynamics [66]. They could similarly prevent loss of alleles, even if the coevolutionary interaction itself would specifically favour only one of the alleles. Obviously, the above processes act in combination with each other in natural systems. Therefore, it is indeed conceivable that recurrent selective sweeps shape longlasting coevolutionary dynamics even in multicellular host systems.

Conclusion
In conclusion, decades of empirical efforts have tried to demonstrate the presence of Red Queen dynamics during host-parasite coevolution. This has led to most ingenious experiments which repeatedly and independently confirmed negative frequency dependence as a driving force [8,9,26,67,68] and such a trend continues to date [21,69]. These studies yielded impressive evidence that parasite abundance typically increases first and, once the host evolves a defense mechanism, it decreases again. However, sustained allele frequency oscillations of a particular allele, as predicted by numerous theoretical models assuming constant population size in the absence of any stochastic effects, have not been reported. We here propose that Lotka-Volterra population size fluctuations and the associated stochastic effects represent an inherent property of host-parasite interactions that can lead to rapid fixation of alleles, even those initially rare, thus preventing sustained oscillations. Consequently, Lotka-Volterra population size fluctuations have the potential to stop the Red Queen -unless novel variants are introduced into the population and/or additional selective constraints maintain allelic diversity. In retrospect, our findings may not be entirely unexpected. However, to date, they have not yet been directly demonstrated using a systematic analysis approach, as implemented here. More importantly, they are generally neglected in the numerous current empirical studies on the topic, in spite of their potential importance. They clearly deserve specific attention in future theoretical and empirical work aimed at an improved understanding of host-parasite coevolutionary dynamics. http://www.biomedcentral.com/1471-2148/13/254

Relating stochastic and deterministic dynamics
Stochastic models are often developed starting from deterministic formulations [42]. Since the same deterministic formulation can be the limiting case of many individual based models, this procedure may be problematic. Instead, beginning from a stochastic, individual based description and then calculating the deterministic analogue will provide only a single direct link between the two approaches and allows for their direct comparison.
We consider a haploid system involving two antagonistic pairs, two alleles in hosts and parasites each. Firstly, all possible changes are written in the form of simple chemical reactions. In our particular case we have eight such possible reactions. We denote the two hosts and the two parasites by H 1 and H 2 and P 1 and P 2 respectively. Thus we have, For instance, a parasite 2 individual dies with the ratec. From these rate reactions, we obtain the transition rates of the system. Depending on the number of individuals of the different types namely n = {n H 1 , n H 2 , n P 1 , n P 2 }, we write the rates as, T(n H 1 + 1, n H 2 , n P 1 , n P 2 |n) =μ n H 1 N avg T(n H 1 , n H 2 − 1, n P 1 , n P 2 + 1|n) = 2b n H 2 N avg n P 2 N avg T(n H 1 , n H 2 , n P 1 − 1, n P 2 |n) =c n P 1 N avg T(n H 1 , n H 2 , n P 1 , n P 2 − 1|n) =c n P 2 N avg where the reaction rates, have been corrected by each reactions combinatorial possibility [70,71] and N avg is the average population size which we consider to be the same for the hosts as well as the parasites (the difference in the average population size can be interpreted as the ratio between the growth rate of hosts and the death rate of parasites). Using these rates, we can write down deterministic differential equations for the change in the average number of a certain type, e.g.
Introducing rescaled reaction rates, μ =μ N avg , b = 2b N avg and c =˜c N avg , we obtain In the limit of a large population size we recover the mean field approximation or the population level model [71], whereḣ 1 = d n H 1 dtN avg and the frequencies of H i and P i are given by h i and p i . In the same way, we can derive deterministic differential equations for the frequencies of the other types,

Stochastic simulations
The Gillespie algorithm gives an exact numerical solution of the Master equation algorithm with the transition rates as defined in Eqs. 6. Since the population size is not constrained, this simulation method includes a stochastic analogue of the Lotka-Volterra cycles. We computationally remove the Lotka-Volterra cycles by culling the population of each species after N avg transitions have taken place. During the N avg transitions the types within a species can evolve to different frequencies. But in the end they are reset to sum up to N avg while maintaining the relative abundances. The Gillespie method is discrete in the number of individuals but continuous in time. The unit of time is the same as in the deterministic system.
Alternatively we can consider a small amount ε of interactions between the otherwise independent Lotka-Volterra interactions. This is then represented by the following set of differential equations, Even for this case, including Lotka-Volterra interactions causes a faster extinction of the Red Queen cycles involving all four types. As an example we provide simulation results where in addition to similar parameters as in Figure 2 we add a ε = 0.1b ( Figure 5). Although the fixation time is elevated as compared to the case with no interactions (Figure 2), they are still not comparable to the extremely high fixations times observed when Lotka-Volterra dynamics is excluded.