Host-parasite coevolution in populations of constant and variable size

The matching-allele and gene-for-gene models are widely used in mathematical approaches that study the dynamics of host-parasite interactions. Agrawal and Lively (Evolutionary Ecology Research 4:79–90, 2002) captured these two models in a single framework and numerically explored the associated time discrete dynamics of allele frequencies. Here, we present a detailed analytical investigation of this unifying framework in continuous time and provide a generalization. We extend the model to take into account changing population sizes, which result from the antagonistic nature of the interaction and follow the Lotka-Volterra equations. Under this extension, the population dynamics become most complex as the model moves away from pure matching-allele and becomes more gene-for-gene-like. While the population densities oscillate with a single oscillation frequency in the pure matching-allele model, a second oscillation frequency arises under gene-for-gene-like conditions. These observations hold for general interaction parameters and allow to infer generic patterns of the dynamics. Our results suggest that experimentally inferred dynamical patterns of host-parasite coevolution should typically be much more complex than the popular illustrations of Red Queen dynamics. A single parasite that infects more than one host can substantially alter the cyclic dynamics.


Background
The antagonistic interaction between hosts and their parasites are of particular interest in ecology and evolution, because they are ubiquitous and usually associated with high selection pressure that affects numerous life history traits. Because of the negative effect of parasites on host fitness, the study of these interactions is of central importance in biomedical [1,2], agricultural [3,4] and species conservation research [5,6]. The exact dynamics of the two coevolving populations are usually evaluated with the help of mathematical models. Among these, models including an explicit genetic description of host-parasite interaction, such as genefor-gene (GfG) and matching-allele (MA) models, are particularly widespread. Genetic interaction is usually *Correspondence: traulsen@evolbio.mpg.de 1 Max Planck Institute for Evolutionary Biology, August-Thienemann-Str. 2, 24306, Plön, Germany Full list of author information is available at the end of the article incorporated by taking into account the current understanding of resistance-infectivity patterns in biological systems. The gene-for-gene (GfG) model was proposed by Flor [7] to capture disease resistance patterns in plants.
Here, a host individual carrying a resistance gene can recognize parasites harboring the corresponding avirulence product and trigger a defense response averting the infection [8]. Inspired by self-nonself recognition in immune systems [9], the matching-allele (MA) model was introduced to reflect host-pathogen interactions in animals. In this case, parasites carrying a certain allele can only invade host individuals with the corresponding allele. By combining predictive power of mathematical modeling and their connection to empirical data, these models successfully served to understand key evolutionary problems. To mention only the most important examples, these models were used to assess the Red Queen hypothesis for the evolution of sexual reproduction [10], the maintenance of genetic diversity by parasite-mediated selection [11], and the role of the cost of resistance/virulence in coevolution [12,13].
Agrawal and Lively [14] developed a general model that interpolates between a pure matching-allele model and a pure gene-for-gene model, as a single parameter is tuned between 0 and 1. This model was introduced for haplotypes of two loci with mutation and recombination. Variance in host and parasite allele frequency was plotted as an evaluation of the time discrete dynamics. The highly dynamical aspects of matching-allele models were observed across most of the MA-GfG continuum. Agrawal and Lively showed that cyclic dynamics of host and parasite genotypes is observed not only in the MA model, but also in all the intermediate models and in the GfG model. This finding indicates that the cyclic Red Queen dynamics does not hinge upon the use of a particular model for host parasite interactions. However, this study was computational and only performed for particular parameter sets due to the complexity of the model. Instead of tackling the dynamics from an analytical perspective to allow for general statements for all parameter sets, subsequent theoretical approaches have increased complexity of the assumed interaction in order to increase the biological realism, for instance by defining a multilocus model that deals with various combinations of MA loci and GfG loci [15].
In this paper, we take the existing and currently widely applied modeling approach [14] and develop a new mathematical framework that allows for an analytical characterization of the involved dynamics. Based on this framework, the aim of our study is to improve our understanding of host-parasite coevolution by investigating both the impact of different types of interactions and the consequence of interaction-dependent population size changes. We simplify the model of Agrawal and Lively and focus on a single locus to keep interaction among loci from interfering with the conclusion, in particular the differences between the GfG model [16] and the MA model [17]. We use the assumptions of Agrawal and Lively [14] inspired by Parker [13] to connect the two popular models by a single parameter, but also provide an alternative, linear interpolation in the discussion. To enhance clarity, we focus on a system with two host genotypes and two parasite genotypes and use their interaction to characterize the evolutionary dynamics of the two populations involved in this reaction. In addition, we depart from the usual assumption of constant population size and apply the Lotka-Volterra equations to acknowledge inter-dependent population dynamics during host-parasite coevolution. To compare the dynamics with a model assuming constant population density, we apply the Replicator Dynamics [18][19][20][21] with the same interaction matrix between hosts and parasites. While the dynamics between the two models is different, it seems to be crucial to understand both constant as well as changing population size, as there are biological examples for both of them.
We conducted a linear stability analysis at the interior fixed point of the resulting nonlinear dynamical system, which indicates critical differences in dynamical patterns between the models of host-parasite coevolution. Either with constant or with changing population size, the population densities oscillate with a single frequency in a pure MA model. In a model deviating from MA, a second oscillation frequency arises with changing population density, but not with constant population size.

Model and Results
We consider haploid hosts and parasites with two alleles on a single locus. Hence, there are two host types and two parasite types that are denoted by H 1 , H 2 , P 1 , and P 2 , respectively. In the simplest case, each parasite type can only infect the corresponding host type. Hence, no host/parasite type is superior to the other. This case corresponds to the matching-allele model, which under the assumption of constant population density is equivalent to the evolutionary game of matching pennies [21,22]. In a GfG model, the virulent parasite P 2 can potentially infect both hosts, the one with susceptible allele H 1 and the one with resistance allele H 2 . Yet, the avirulent parasite P 1 can only infect the susceptible host H 1 , as the host H 2 with the resistance allele can prevent infection by P 1 . Thus, there is an advantage to the virulent parasite and the resistant host. To maintain the different types in the population, intrinsic costs of virulence and resistance have been suggested [23]. Figure 1 illustrates the fitness of the two parasites on each host for the MA and the GfG model and also for two intermediate cases, where the parasite P 2 can "partially" infect the host H 1 .
We simplified the model of Agrawal and Lively [14] by regarding only one locus. The interactions between hosts and parasites can be expressed with two matrices (corresponding to a bi-matrix game in evolutionary game theory). We are interested in the role of population size changes, determined by the birth and death rates of host and parasite. We assume that the interaction between hosts and parasites affects the death rates of hosts and the birth rates of parasites. The birth rate of hosts and the death rates of parasites are chosen either as constant parameters (such that the population size changes) or in a way that ensures that the total population sizes of hosts and parasites remain constant, see below.
For the parasite, we assume that the interactions with the hosts increase birth rates according to the matrix ( 1 ) Fig. 1 Illustration of parasite fitness. Fitness of avirulent parasite P 1 and virulent parasite P 2 on the two hosts H 1 and H 2 for the matching-allele model (α = 0, top), the gene-for-gene model (α = 1, bottom), and two intermediate models (α = 1/3 and α = 2/3). Gray areas represent the fitness reduction for P 2 due to the cost of virulence κ = 1/2, which is ακσ in H 2 (Eq. (1b)), hence, σ/2 in GfG model. In H 1 the fitness reduction for P 2 due to the cost of virulence is α 2 κσ For example, the birth rate of parasite P 1 in interactions with H 1 is σ , while it is 0 when P 1 interacts with H 2 (see Appendix A.3 for a generalization of these assumptions). The maximum virulence of the parasite is given by σ . The parameter κ describes the cost for the parasite virulence, as usually assumed in the GfG model. This model interpolates between the MA and the GfG model as the parameter α is varied between 0 (MA) and 1 (GfG). From the interaction matrix Eq. (1), the birth rates for the parasites are given by For the host, we assume that the interactions increase the death rate according to the matrix where the parameter γ describes the cost for the host resistance. This leads to the host death rates We assume a large population size and focus on the change in population densities. The population densities of the two host and two parasite types are given by h 1 , h 2 , p 1 , and p 2 , respectively. The population dynamics of the hosts and parasites are driven by their respective birth and death rates and can be captured by a set of differential equationsḣ where b h is the birth rate of both hosts, and d p is the death rate of both parasites. We will choose the host birth rate b h and parasite death rate d p in two distinct ways. Our first approach assumes constant values for b h and d p , which leads to a host/parasite population that is changing in size. This corresponds to the competitive Lotka-Volterra dynamics. The second approach focuses on relative abundances of host and parasite alleles and implies a normalization of the two population sizes. This corresponds to the Replicator Dynamics for an evolutionary role game between hosts and parasites, which implies two constant population sizes in our context.

Changing population size induced by interactions
With constant host birth rate b h and parasite death rate d p , inserting the host parasite interactions Eqs. (2) and (6) into the dynamical system Eq. (5) leads tȯ This model, which we analyze in detail below, results in changes in the population sizes of both hosts and parasites. In particular, the changes are caused by the antagonistic interactions between the hosts and the parasite -as a consequence of the Lotka-Volterra relationship.

Constant population size
To obtain a model of constant population size that is comparable to the one described above, we retain the interaction matrices and adjust the host birth rate and parasite death rate to maintain the population size. Requiring constant h 1 + h 2 and constant p 1 + p 2 impliesḣ 1 +ḣ 2 = 0 as well asṗ 1 +ṗ 2 = 0. This leads to The normalization h 1 + h 2 = 1 implies that a single equation for h 1 is sufficient to describe the dynamics for the host. Similarly, due to the normalization p 1 + p 2 = 1 the parasite dynamics are fully captured by tracking p 1 .
Inserting the dynamical host birth and parasite death rates in the dynamical system Eq. (2), the equations become identical to the Replicator Dynamics (RD) [19,21,24], While the death rates of the host still depend on the parasites and the birth rates of the parasites still depend on the hosts, the dynamics of this system is in general less complex than in the case of changing population size, as it is only two-dimensional.

Population dynamics
To obtain first information about the population dynamics, we calculated the trajectories of the system numerically for a particular set of parameters. In addition, we identify the fixed points of the differential equations and study their stability to gain insight into the coevolutionary dynamics for all parameter sets. More specifically, we can use a linear stability analysis of the unique interior fixed point to infer the dynamical patterns arising in this system [16,25]. Finally, we also assess constants of motion.

Numerical solution of the dynamics
To illustrate the differences in the population dynamics described in Eqs. (6) and (8), we show numerical solutions side by side in Fig. 2.
The dynamics in models with constant host and parasite population sizes resemble the common Red Queen pattern. Under changing population sizes the system is uncoupled into two independent host-parasite pairs in a pure MA model. As the model deviates from the MA model with increasing α, the dynamics becomes more complex, since the four population densities of the types P 1 , P 2 , H 1 , and H 2 are coupled.

Stability of boundary fixed points
The fixed points of the system are the points where all population sizes remain constant in time,ḣ 1 =ḣ 2 =ṗ 1 = p 2 = 0. The position of the fixed points and their stability change with changing parameters.
For the Lotka-Volterra dynamics, a trivial fixed point is (h 1 , h 2 , p 1 , p 2 ) = (0, 0, 0, 0) where both the hosts and parasites are absent, cf. Eq. (6). Additionally, extinction of one host and the associated parasite leads to two further fixed points, (h 1 , h 2 , p 1 , In gene-for-gene-like models, α > 0, the susceptible host H 1 and the virulent P 2 can coexist in the absence of H 2 and (1−ακ) ). The opposite case, coexistence between H 2 and P 1 in the absence of H 1 and P 2 is not possible, as our host-parasite interaction model assumes that the birth rate of P 1 is zero in the absence of H 1 . A linear stability analysis of the Lotka-Volterra model shows that all boundary fixed points are unstable for αγ < σ . That is, if the cost of resistance αγ (which is scaled by the amount of GfG influence) is less than the maximum host fitness reduction caused by infection σ , then all host and parasite types will coexist. The Replicator Dynamic system, Eq. (8), has four fixed points at the boundaries, each is reflecting fixation of one host and one parasite: 1). A linear stability analysis reveals that all these fixed points are unstable.

Stability of the interior fixed point
In addition to the boundary fixed points, the system has a unique fixed point in the interior. In the Lotka-Volterra system, we obtain a non-trivial fixed point of the four dimensional dynamical system described in Eq. (6) when αγ < σ . This fixed point, where all types coexist, is given by Fig. 2 Example of population dynamics based on the Lotka-Volterra equations (left) and the Replicator Dynamics (right). While the dynamics on the right side resembles the common Red Queen pattern, the left side is more complex. In a pure matching-allele model (top), the plot on the left shows two independent sets of Lotka-Volterra dynamics, one for H 1 and P 1 (blue and red solid lines, correspondingly) and a second one for H 2 and P 2 (blue and red dotted lines). As the model deviates from MA model with increasing α (rows 2-4) more complicated dynamics arise, since the four population densities of H 1 , H 2 , P 1 , and P 2 are coupled (parameters γ = 0.005, κ = 0.5, and σ = 0.01 for both Lotka-Volterra and Replicator Dynamics. Host birth rate b h = 1.5 and parasite death rate d p = 1.0 in the Lotka-Volterra case. Initial population densities h 1 = p 1 = 150, For αγ > σ , the resistant host is always disadvantageous because of the high cost of the resistance allele (γ ). Consequently, extinction of H 2 and P 2 then becomes a stable fixed point. For αγ < σ , h * 1 and h * 2 increase linearly with parasites' death rate d p , while p * 1 and p * 2 increase linearly with hosts' birth rate b h . A linear stability analysis of the interior fixed point (see Appendix A.1 for details) shows that the equilibrium is neutrally stable. Close to the interior fixed point, the system exhibits undamped oscillations. More specifically, the four eigenvalues of the Jacobi-matrix are two distinct pairs of complex conjugates without real parts. This means there are two distinct oscillation frequencies in the system, where measures the ratio between the two oscillation frequencies. This ratio decreases when we move away from the MA interaction model. For α 1 we find In particular, for the MA model both oscillation frequencies collapse into a single one. However, all solutions for α > 0 exhibit both of the frequencies (Fig. 3).
For the Replicator Dynamics system, in which the population size is constant, the non-trivial fixed point of Eq. (8) is given by A linear stability analysis shows that the interior fixed point is again neutrally stable, as the two eigenvalues are a pair of purely imaginary, complex conjugated numbers when αγ < σ (see Appendix A.2 for details). Hence, there is only one characteristic oscillation frequency of the dynamical system at the fixed point, l has a maximum value, σ/(4π), in the pure matching-allele model (α = 0). Close to the matching-allele model, α 1 the oscillation frequency decreases with increasing α as The solutions around the fixed point exhibit the oscillation frequency described by Eq. (14). The trajectories are closed circles as shown on the right side of Fig. 3.

Disentangling evolutionary and ecological dynamics
To clarify the ecological effect on the dynamics, particularly at the interior fixed point, we derive the dynamics of the host and parasite population sizes, h = h 1 + h 2 and p = p 1 + p 2 , and the relative abundance of H 1 and P 1 in the population, x = h 1 /h and y = p 1 /p, from Eq. (2). According to Eq. (2) the differential equations for the population sizes of hosts h and parasites p arė assuming fully general interaction matrices. The differential equations for relative abundances of H 1 and P 1 arė (see Eq. (40) in Appendix A.3) and the geometric mean of host and parasite population size h * · p * , i.e., (calculated from Eq. (33) in Appendix A.3). Thus, one of the oscillations results purely from ecological interactions, while the other one arises from the combination of ecology and evolution in our system.

Constants of motion
The system with constant population size has a constant of motion (Eq. (10.22) in [21]) given by Due toL = 0, we obtain sustained oscillations for any initial condition, even far away from the interior fixed point Eq. (13) The case of changing population size is more intricate. In the case of a matching allele model α = 0, the two equations decouple and we have two independent Lotka-Volterra systems with sustained oscillations, characterized by the two constants of motion (24b) While we do not find a constant of motion for the general case of α > 0, particular initial conditions can lead to invariants. If the initial condition fulfills and (25a) which corresponds to a two-dimensional subspace, then there are two constants that remain invariant over time, (26b) Note that with the condition Eq. (25a) the ratio p 1 /p 2 remains constant and with the condition Eq. (25b), the ratio h 1 /h 2 remains constant. This shows that the nature of the dynamics in this case does not only depend on the choice of parameters, but also on the initial state of the system, which in principle leads to a further complication for the corresponding experimental systems.

Short overview
Host-parasite interactions are acknowledged as a driving evolutionary force promoting biological diversity and sexual reproduction [10,11], with the MA and GfG model being the most popular models to describe the genetic interaction for coevolving hosts and parasites [26][27][28][29][30][31][32]. Despite a number of important insights provided within their framework, the generality of findings often suffers from the complexity of the models employed and, as a consequence, the difficulty to fully understand them analytically [33].
In this study, we present a very general yet parsimonious model of host-parasite coevolution spanning from MA to GfG with either constant or interaction-driven changing population size. Derived analytical solutions revealed that the coevolution dynamics differs qualitatively between the models with constant and changing population sizes. Apart from the pure MA situation, the well known Red Queen dynamics with closed trajectories is only observed in models with constant population size. This implies that the patterns of host-parasite dynamics to be expected in real biological systems can be much more intricate than suggested by the most popular theoretical models.

Main results and analytical solution
Our study is based on a simplification of the model suggested by Agrawal and Lively [14] that explores a continuum between the MA and GfG models. We study the model in the context of haplotypes with a single locus, but relax the restriction to constant population size. With a coevolutionary system of two host and two parasite types we achieved an analytical characterization across the entire parameter space. To study ecological effects caused by the victim-exploiter interaction [34] between hosts and parasites, we consider models with changing population size aside of models with constant population size. Under the assumption of constant population size, the dynamics in MA and GfG models appear to be very similar, both showing sustained oscillations with only one oscillation frequency. Yet, introducing changing population size according to the Lotka-Volterra equations, we obtain distinct patterns of the population dynamics. For changing population sizes, a single oscillation frequency is present only in the MA model. An additional oscillation frequency arises for all other points on the MA-GfG continuum in that case. In other words, changing population size leads to a much more complex dynamics in GfG-like models, but not in the pure MA model.
In the stochastic model analyzed in [29], the analysis of allele fixation time for the MA model revealed that Lotka-Volterra dynamics in combination with the associated stochastic effects quickly break down the Red Queen circle. As the dynamics in GfG-like models take a completely different nature with changing population size, the influence of Lotka-Volterra dynamics on the Red Queen circle is yet unclear and remains to be assessed in more detail in the future, especially as our current analysis did not take stochastic effects into account.

Generality of results
To test the generality of our findings we additionally analyzed the interaction matrix suggested by Parker [13] (Eq. (37)). There, a factor that denotes the fitness reduction of the avirulent parasite encountering the resistant host and an advantage of the virulent parasite meeting the resistant host are assumed in addition. These two parameters together with the costs of resistance and virulence determine whether the model is MA or GfG. Again we obtain two distinct oscillation frequencies for the population dynamics with changing population sizes in GfG-like models (the ratio is shown in Eq. (38) in the Appendix A.3).
Despite the convincing biological relevance of the interaction matrix elements in [14], they do not change monotonically on the MA-GfG continuum, e.g., with a cost of virulence κ > 0.5, M p 21 in Eq. (1) first increases then decreases as α increases from 0 to 1. As an alternative interpolation, we therefore also considered interaction matrices that describe a linear transition from MA to GfG model, such that The analysis in Appendix A.4 shows that our conclusion also holds for the linear interpolation. One should keep in mind that both MA and GfG models and even the intermediate models proposed by Parker, Agrawal & Lively, or us are only a small subset of the possible models for host-parasite interaction. An observation that will hold for any such model is that as long as the population sizes are kept constant, the population dynamics follows a closed circle with a single oscillation frequency. However, with changing population size a second oscillation frequency arises when the model become GfG-like, which can lead to much more intricate dynamics. For a pure MA model or an inverse MA model (where the diagonal instead of the off-diagonal matrix elements are zero), there still is only one oscillation frequency (see Eq. (36) in Appendix A.3).

Impact of eco-evo feedback in genetically explicit models
In the last two decades it has been realized that evolutionary changes can be faster than previously thought and, thus, occurring on the same time-scale as ecological interactions, especially in case of coevolving hosts and parasites [35][36][37][38]. Population dynamics can influence the pace of coevolution via so called eco-evolutionary feedbacks, or even give rise to a new type of coevolutionary dynamics as we showed in our study. Interestingly enough, a comprehensive part of the theoretical studies on eco-evolutionary feedbacks is conducted within the framework of evolutionary game theory and adaptive dynamics [21,39]. In contrast to our model, these approaches usually do not include an explicit definition of genetic interaction between the species, which limits their application for interpreting patterns of genetic variability in natural populations [40]. Rapid changes in genetic composition may lead to perturbation in host demography and disease dynamics, as was observed for the myxoma virus epidemic in Australian populations of European rabbit [41]. Genetic adaptation can improve overall population fitness and "buffer" the unfavorable impact of pathogens (evolutionary rescue) [42]. However, population perturbations may constrain adaptability, for example, via enhancing inbreeding, affecting trait heritabilities and disturbing allele composition irrespective of natural selection [43][44][45][46]. Thus, models accounting simultaneously for the genetic basis of host-parasite interaction and associated population dynamics may be necessary to fully understand ongoing coevolution among species and the effect it would have on genetic diversity. We are aware of only a few such models [29,[47][48][49][50][51], and most of them confirm that ecological parameters can have a very strong effect on coevolution.

Implications for maintenance of genetic diversity
Numerous field studies identified the presence of comprehensive heritable variation in resistance-infectivity patterns for plant and animal populations and their respective pathogens, suggesting that coevolution acts to maintain genetic diversity [3,11,[52][53][54][55]. However, already the first studies, which attempted to explain such variation by cycling dynamics, encountered the problem of stability. This is especially true for the GfG model, as a parasite with the virulent allele would be quickly fixed, unless having a cost of virulence [3,12,56]. In addition to the cost, other factors have been examined for their potential role in maintaining variation, including epidemiological feedback [51,57], spatial structure [48,49,58,59], genetic drift [60], diffuse multi-species coevolution [61], models with multiple alleles and multiple loci [16,60,62]. Several studies proposed that multiple factors need to act jointly for long-term coexistence of multiple resisto-and infectotypes [33]. The view of a multifactorial basis of the maintenance of diversity creates an additional challenge for theoretical and empirical studies to disentangle them. As opposed to that, Tellier [34] presented a simple GfG framework showing that the general condition for stability is the presence of direct frequency-dependent selection (where fitness of an allele declines with increasing frequency of that allele itself ). In this context, the distinction is made between direct frequency dependence and indirect frequency-dependent selection where fitness is mediated by the frequency of the corresponding antagonist. Direct frequency-dependent selection can be introduced in the model by incorporation of epidemiological or ecological factors ( [32], Table 1). If we introduce a direct frequency-dependent element by applying competitive Lotka-Volterra equations or the concept of empty spaces [63] (implying the existence of a carrying capacity) into our model, the neutrally stable interior fixed point becomes stable. Instead of forming tori or moving along closed circles, the deterministic trajectory spirals inwards. In this case, the oscillation of allele frequencies lasts longer in stochastic simulations, hence the polymorphic state is more stable.
The stability analysis derived the condition for coexistence αγ < σ , suggesting that departing from the GfG end of the continuum would increase a range of parameters at which the oscillation of allele frequencies is maintained. Therefore, patterns of "partial" infectivity by a virulent parasite are more likely to result in cycling dynamics compared to a pure GfG situation. Agrawal and Lively [14] came to the same conclusion by evaluating computational simulations. This reinforces the importance of exploring dynamics for intermediate points on the MA-GfG continuum, especially as experimental studies provide some examples of such types of interaction [64]. In contrast to [34] and many other studies [14,16,59], our model is implemented on a continuous time-scale and, therefore, covers host and parasite systems with overlapping generations. Interestingly, it has been proposed that models with discrete generations would favor coevolutionary cycling by synchronizing ecological and epidemiological processes [51], while in [34] the condition for stable cycling is more restrictive for discrete generations when compared to the continuous model.

Conclusions
In summary, we have shown that only a small and possibly biased subset of possible host-parasite interaction dynamics is captured by the mathematical models that assume fixed population size or particular genetics for the interaction, such as the MA model. Even in a simple model that allows for a full analytical description, the dynamics can vary substantially between subsequent coevolutionary cycles. We demonstrate analytically that the complex dynamics found for changing population sizes is not a result of choosing a particular interaction matrix. The complex pattern is not limited to the set of models considered here, but rather a general property of models beyond fixed population size. Our findings highlight the importance of the interconnectedness between coevolution and population dynamics, and its potential role in understanding the generation and maintenance of genetic variation. Our model provides a solid framework that can be extended to more realistic (i.e., usually more complex) scenarios in the future, including the specification of alternative virulence components and their differential effect on host fitness, or consideration of population carrying capacity with its limiting effect on growth rate.

A.1 Stability of the interior fixed point in the Lotka-Volterra dynamics
In order to analyse the system at the interior fixed point (h * 1 , h * 2 , p * 1 , p * 2 ), we first linearise the system around this point. For general points (h 1 , h 2 , p 1 , p 2 ), the linearised system is given by by the Jacobian matrix J(h 1 , h 2 , p 1 , The eigenvalues of this matrix determine linear stability at the fixed point [25]. If there is at least one eigenvalue with positive real part, the point would be unstable. If all eigenvalues have negative real parts, the point would be stable. In our case, the four eigenvalues are 1,2 = ±i b h d p and (29) Except the term √ αγ − σ , the remaining factors in in 3,4 are positive. For αγ > σ , allele H 1 is always beneficial.
Consequently, the fixed point is unstable as one of the eigenvalues 3 or 4 is positive. For αγ < σ , the fixed point is a center with neutral stability as all eigenvalues are purely imaginary. Only the case of αγ < σ is of further interest in this manuscript, as the result is straightforward in the opposite case.

A.3 Stability of the interior fixed point for general interaction matrices
The appearance of the second oscillation frequency at the interior fixed point in gene-for-gene-like models with changing population sizes does not depend on the exact choice of the interaction matrices in Eq. (2). To show this, we recalculate the interior fixed point and apply linear stability analysis on interaction matrices of a general form,