Open Access

Host-parasite coevolution in populations of constant and variable size

  • Yixian Song1,
  • Chaitanya S Gokhale2,
  • Andrei Papkou3,
  • Hinrich Schulenburg3 and
  • Arne Traulsen1Email author
BMC Evolutionary Biology201515:212

Received: 19 December 2014

Accepted: 21 August 2015

Published: 29 September 2015



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.


Matching-allele Gene-for-gene Lotka-Volterra equation Replicator dynamics. Red Queen hypothesis Stability analysis


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 gene-for-gene (GfG) and matching-allele (MA) models, are particularly widespread. Genetic interaction is usually 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 multi-locus 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 [1821] 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 \(\mathcal {H}_{1}\), \(\mathcal {H}_{2}\), \(\mathcal {P}_{1}\), and \(\mathcal {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 \(\mathcal {P}_{2}\) can potentially infect both hosts, the one with susceptible allele \(\mathcal {H}_{1}\) and the one with resistance allele \(\mathcal {H}_{2}\). Yet, the avirulent parasite \(\mathcal {P}_{1}\) can only infect the susceptible host \(\mathcal {H}_{1}\), as the host \(\mathcal {H}_{2}\) with the resistance allele can prevent infection by \(\mathcal {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 \(\mathcal {P}_{2}\) can “partially” infect the host \(\mathcal {H}_{1}\).
Fig. 1

Illustration of parasite fitness. Fitness of avirulent parasite \(\mathcal {P}_{1}\) and virulent parasite \(\mathcal {P}_{2}\) on the two hosts \(\mathcal {H}_{1}\) and \(\mathcal {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 \(\mathcal {P}_{2}\) due to the cost of virulence κ=1/2, which is α κ σ in \(\mathcal {H}_{2}\) (Eq. (1b)), hence, σ/2 in GfG model. In \(\mathcal {H}_{1}\) the fitness reduction for \(\mathcal {P}_{2}\) due to the cost of virulence is α 2 κ σ

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
$$ \begin{array}{ll} &\begin{array}{cc} \qquad\;\;\mathcal{H}_{1} & \,\quad\qquad\mathcal{H}_{2} \end{array}\\ \mathcal{M}^{p} = \begin{array}{l} \mathcal{P}_{1}\\ \mathcal{P}_{2} \end{array}& \left(\begin{array}{cc} \sigma & 0\\ \alpha (1 - \alpha \kappa) {\sigma} & (1- \alpha \kappa) {\sigma} \end{array} \right). \end{array} $$
For example, the birth rate of parasite \({\mathcal {P}_{1}}\) in interactions with \({\mathcal {H}_{1}}\) is σ, while it is 0 when \({\mathcal {P}_{1}}\) interacts with \({\mathcal {H}_{2}}\) (see Appendix A.3 Stability of the interior fixed point for general interaction matrices 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
$$\begin{array}{@{}rcl@{}} b_{\mathcal{P}_{1}}\!\! &=& \!\!\mathcal{M}_{11}^{p} h_{1} +\mathcal{M}_{12}^{p} h_{2}= \sigma \; h_{1}\\ b_{\mathcal{P}_{2}}\!\! &=& \!\!\mathcal{M}_{21}^{p} h_{1} +\mathcal{M}_{22}^{p} h_{2}= \alpha(1-\alpha\kappa) \sigma\; h_{1} + (1-\alpha\kappa)\sigma \; h_{2} \end{array} $$
For the host, we assume that the interactions increase the death rate according to the matrix
$$ \begin{array}{ll} &\begin{array}{cc} \quad\mathcal{P}_{1} & \qquad\qquad\qquad\mathcal{P}_{2} \end{array}\\ \mathcal{M}^{h} = \begin{array}{l} \mathcal{H}_{1}\\ \mathcal{H}_{2} \end{array}& \left(\begin{array}{cc} -\sigma & - \alpha (1 - \alpha \kappa) \sigma\\ - \alpha \gamma & (1 - \alpha \gamma) (1 - (1 - \alpha \kappa) \sigma)-1 \end{array} \right), \end{array} $$
where the parameter γ describes the cost for the host resistance. This leads to the host death rates
$$\begin{array}{@{}rcl@{}} d_{\mathcal{H}_{1}} &=& \mathcal{M}_{11}^{h} p_{1} +\mathcal{M}_{12}^{h} p_{2}= -\sigma \; p_{1} -\alpha (1-\alpha\kappa) \sigma \; p_{2}\\ d_{\mathcal{H}_{2}} &=& \mathcal{M}_{21}^{h} p_{1} +\mathcal{M}_{22}^{h} p_{2}\\&=& -\alpha\gamma \; p_{1} + ((1-\alpha\gamma) (1-(1-\alpha\kappa)\sigma)-1) \; p_{2} \end{array} $$
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
$$\begin{array}{*{20}l} \dot{h_{1}} &= h_{1} (b_{h}+d_{\mathcal{H}_{1}}) \\ \dot{h_{2}} &= h_{2} (b_{h}+d_{\mathcal{H}_{2}}) \end{array} $$
$$\begin{array}{*{20}l} \dot{p_{1}} &= p_{1} (b_{\mathcal{P}_{1}} - d_{p}) \\ \dot{p_{2}} &= p_{2} (b_{\mathcal{P}_{2}} - d_{p})\, , \end{array} $$

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 to
$$\begin{array}{*{20}l} \dot{h_{1}} &= h_{1} \left(b_{h} -p_{1} \sigma - p_{2} \alpha (1-\alpha\kappa) \sigma \right) \\ \dot{h_{2}} &= h_{2} \left(b_{h} -p_{1} \alpha\gamma - p_{2} \left((1-\alpha\gamma) (1-\alpha\kappa)\sigma+\alpha\gamma \right)\right) \end{array} $$
$$\begin{array}{*{20}l} \dot{p_{1}} &= p_{1} (\sigma h_{1} - d_{p}) \\ \dot{p_{2}} &= p_{2} \left(\sigma \left(h_{1} \alpha(1-\alpha\kappa) + h_{2} (1-\alpha\kappa)\right) - d_{p} \right)\, . \end{array} $$

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 \(\dot {h_{1}}+\dot {h_{2}} = 0 \) as well as \( \dot {p_{1}}+\dot {p_{2}} = 0\). This leads to
$$\begin{array}{*{20}l} b_{h} &=-\frac{h_{1} d_{\mathcal{H}_{1}} + h_{2} d_{\mathcal{H}_{2}}}{h_{1}+h_{2}} \end{array} $$
$$\begin{array}{*{20}l} d_{p} &=\frac{p_{1} b_{\mathcal{P}_{1}} + p_{2} b_{\mathcal{P}_{2}}}{p_{1}+p_{2}}\, . \end{array} $$
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],
$$\begin{array}{*{20}l} \dot{h_{1}} &= h_{1} (1-h_{1}) (d_{\mathcal{H}_{1}}-d_{\mathcal{H}_{2}}) \end{array} $$
$$\begin{array}{*{20}l} \dot{p_{1}} &= p_{1} (1-p_{1}) (b_{\mathcal{P}_{1}}-b_{\mathcal{P}_{2}})\, . \end{array} $$

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.
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 \(\mathcal {H}_{1}\) and \(\mathcal {P}_{1}\) (blue and red solid lines, correspondingly) and a second one for \(\mathcal {H}_{2}\) and \(\mathcal {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 \(\mathcal {H}_{1}\), \(\mathcal {H}_{2}\), \(\mathcal {P}_{1}\), and \(\mathcal {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, h 2=p 2=50)

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 \(\mathcal {P}_{1}\), \(\mathcal {P}_{2}\), \(\mathcal {H}_{1}\), and \(\mathcal {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, \(\dot {h_{1}} =\dot {h_{2}} = \dot {p_{1}} = \dot {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},p_{2})=(\frac {d_{p}}{ \sigma }, 0, \frac {b_{h}}{\sigma }, 0)\) and \((h_{1},h_{2},p_{1},p_{2})=(0, \frac {d_{p}}{ \sigma (1-\alpha \kappa)}, 0,\frac {b_{h}}{\alpha \gamma (1-\sigma)+\sigma (1-\alpha \kappa (1-\alpha \gamma))})\). In gene-for-gene-like models, α>0, the susceptible host \(\mathcal {H}_{1}\) and the virulent \(\mathcal {P}_{2}\) can coexist in the absence of \(\mathcal {H}_{2}\) and \(\mathcal {P}_{1}\), \((h_{1},h_{2},p_{1},p_{2})=(\frac {d_{p}}{\alpha \sigma (1-\alpha \kappa)}, 0, 0, \frac {b_{h}}{\alpha \sigma (1-\alpha \kappa)})\). The opposite case, coexistence between \(\mathcal {H}_{2}\) and \(\mathcal {P}_{1}\) in the absence of \(\mathcal {H}_{1}\) and \(\mathcal {P}_{2}\) is not possible, as our host-parasite interaction model assumes that the birth rate of \(\mathcal {P}_{1}\) is zero in the absence of \(\mathcal {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: (h 1,p 1)=(0,0), (h 1,p 1)=(1,0), (h 1,p 1)=(0,1), (h 1,p 1)=(1,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
$$\begin{array}{*{20}l} h_{1}^{*} &= \frac{1}{\sigma} d_{p} \\ h_{2}^{*} &= \frac{1}{\sigma}\frac{1-\alpha (1- \alpha\kappa)}{1 -\alpha\kappa }d_{p} \end{array} $$
$$\begin{array}{*{20}l} p_{1}^{*} &= \frac{1}{\sigma}\frac{\sigma(1-\alpha)(1-\alpha\kappa) + \alpha\gamma(1-\sigma(1-\alpha\kappa)}{\sigma (1-\alpha \gamma) (1-\alpha \kappa)+\alpha \gamma (1-\alpha (1-\alpha \kappa))} b_{h}\\ \\ p_{2}^{*} &= \frac{1}{\sigma}\frac{(\sigma -\alpha \gamma)}{\sigma (1-\alpha \gamma) (1-\alpha \kappa)+\alpha \gamma (1-\alpha (1-\alpha \kappa))}b_{h}\, . \end{array} $$
For α γ>σ, the resistant host is always disadvantageous because of the high cost of the resistance allele (γ). Consequently, extinction of \(\mathcal {H}_{2}\) and \(\mathcal {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 Stability of the interior fixed point in the Lotka-Volterra dynamics 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,
$$ \frac{1}{2\pi}\sqrt{b_{h} d_{p}} \quad\text{and}\quad \frac{m}{2\pi}\sqrt{b_{h} d_{p}}\, , $$
$$ {\fontsize{8.7}{6}\begin{aligned} m=&\,\frac{\sqrt{\sigma (1-\alpha) (1-\alpha \kappa)+\alpha \gamma (1-\sigma(1-\alpha\kappa))} \sqrt{1-\alpha (1-\alpha \kappa)}}{\sqrt{\sigma} \sqrt{\sigma (1-\alpha \gamma) (1-\alpha \kappa)+\alpha \gamma (1-\alpha (1-\alpha \kappa)) }}\\&\times\sqrt{\sigma-\alpha \gamma } \end{aligned}} $$
measures the ratio between the two oscillation frequencies. This ratio decreases when we move away from the MA interaction model. For α1 we find
$$ m\approx 1-\alpha\left(1+\frac{\gamma}{2\sigma}\right)\, . $$
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).
Fig. 3

Trajectories close to the interior fixed points (black points) on the h 1p 1 plane (dark green solid lines both for LV and RD equations) and the h 2p 2 plane (light green dashed lines LV only). The black crosses mark the initial conditions. The black rectangle represent a special set of initial condition while the black solid/dashed lines show the corresponding trajectories. With Replicator Dynamics the h 1p 1 trajectory is a closed circle. With Lotka-Volterra dynamics, the trajectories are closed circle when the initial conditions fulfill Eq. (25) (black lines). For the closed circles (black in LV and green in RD) the initial host population densities, h 1 and h 2 are 5 % above the corresponding fixed point, while the parasite population densities are 5 % beneath the fixed point. Except for α=0 (MA) the green trajectories with LV resemble tori instead of closed circles, an implication for two oscillation frequencies. To show the shift of the interior fixed point as α increases from 0 to 1, the trajectories are plotted all in the same coordinate system at the bottom

For the Replicator Dynamics system, in which the population size is constant, the non-trivial fixed point of Eq. (8) is given by
$$\begin{array}{*{20}l} h_{1}^{*}&= \frac{1-\alpha \kappa }{ 2 -\alpha \kappa-\alpha(1-\alpha\kappa)} \end{array} $$
$$\begin{array}{*{20}l} p_{1}^{*}&= \frac{\alpha \gamma (1-\sigma (1-\alpha \kappa))+\sigma (1-\alpha) (1-\alpha \kappa)}{\sigma ((1-\alpha \gamma (1-\alpha \kappa))+(1-\alpha) (1-\alpha \kappa))}\, . \end{array} $$
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 Stability of the interior fixed point in the Replicator Dynamics for details). Hence, there is only one characteristic oscillation frequency of the dynamical system at the fixed point,
$$ l=\frac{\sqrt{(1-\alpha \kappa) (1-\alpha (1-\alpha \kappa)) (\alpha \gamma (1-\sigma (1-\alpha \kappa))+(1-\alpha) \sigma (1-\alpha \kappa))}}{\sqrt{(1+(1-\alpha) (1-\alpha \kappa)) (1-\alpha \gamma (1-\alpha \kappa)+(1-\alpha) (1-\alpha \kappa))}} \frac{\sqrt{\sigma -\alpha \gamma}}{2\pi} $$
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
$$ l\approx \frac{\sigma}{4\pi}-\frac{\alpha}{16\pi}(2+\gamma+2\kappa)\sigma \, . $$

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 \(\mathcal {H}_{1}\) and \(\mathcal {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 are
$$\begin{array}{*{20}l} \dot{h} &= h(p\,\, f\!(x,y)+b) \end{array} $$
$$\begin{array}{*{20}l} \dot{p} &= p(h\,\, g\!(x,y)-d) \end{array} $$
$$\begin{array}{*{20}l} f(x,y) &= \mathcal{M}_{11}^{h} x y+\mathcal{M}_{12}^{h} x (1-y)+\mathcal{M}_{21}^{h} (1-x) y+\mathcal{M}_{22}^{h} (1-x) (1-y)\\ g(x,y) &= \mathcal{M}_{11}^{p} y x + \mathcal{M}_{12}^{p} y (1 - x) + \mathcal{M}_{21}^{p} (1 - y) x + \mathcal{M}_{22}^{p} (1 - y) (1 - x), \end{array} $$
assuming fully general interaction matrices. The differential equations for relative abundances of \(\mathcal {H}_{1}\) and \(\mathcal {P}_{1}\) are
$$\begin{array}{*{20}l} \dot{x} &= p x(1-x) ((\mathcal{M}_{11}^{h}-\mathcal{M}_{21}^{h}) y + (\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h})(1-y)) \end{array} $$
$$\begin{array}{*{20}l} \dot{y} &= h y(1-y) ((\mathcal{M}_{11}^{p}-\mathcal{M}_{21}^{p})x+ (\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p})(1-x))\, , \end{array} $$

If f(x,y) and g(x,y) are constant in time Eq. (16) yield simple Lotka-Volterra dynamics, while Eq. (18) result in Replicator Dynamics with rescaled time if the population sizes are kept constant.

At the interior fixed point one of the oscillation frequencies, \(\sqrt {b_{h} d_{p}}/(2\pi)\), results solely from Lotka-Volterra dynamics. The other oscillation frequency
$$ \frac{\sqrt{b_{h} d_{p}}}{2 \pi} m = \frac{\sqrt{b_{h} d_{p}}}{2 \pi} \frac{\sqrt{(\mathcal{M}_{11}^{h}-\mathcal{M}_{21}^{h}) (\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h}) (\mathcal{M}_{11}^{p}-\mathcal{M}_{21}^{p}) (\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p})}}{\sqrt{(\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h})} \sqrt{ (\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p})}} $$
(see Eq. (36) in Appendix A.3 Stability of the interior fixed point for general interaction matrices) is the product of the oscillation frequency with constant population size
$$ l=\frac{\sqrt{(\mathcal{M}_{11}^{h} - \mathcal{M}_{21}^{h}) (\mathcal{M}_{22}^{h} - \mathcal{M}_{12}^{h}) (\mathcal{M}_{11}^{p} - \mathcal{M}_{21}^{p}) (\mathcal{M}_{22}^{p} - \mathcal{M}_{12}^{p})}}{2\pi\sqrt{(\mathcal{M}_{11}^{h} + \mathcal{M}_{22}^{h}- \mathcal{M}_{12}^{h} -\mathcal{M}_{21}^{h}) (\mathcal{M}_{12}^{p} + \mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p} - \mathcal{M}_{22}^{p})}} $$
(see Eq. (40) in Appendix A.3 Stability of the interior fixed point for general interaction matrices) and the geometric mean of host and parasite population size \(\sqrt {h^{*}\cdot p^{*}}\), i.e.,
$$ \frac{m\sqrt{b_{h} d_{p}}}{2 \pi} =l\sqrt{h^{*}\cdot p^{*}}\, , $$
$$\begin{array}{*{20}l} h^{*}&=\frac{d (\mathcal{M}_{11}^{p}-\mathcal{M}_{12}^{p}-\mathcal{M}_{21}^{p}+\mathcal{M}_{22}^{p})}{\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}} \end{array} $$
$$\begin{array}{*{20}l} p^{*}&=\frac{b (\mathcal{M}_{12}^{h}+\mathcal{M}_{21}^{h}-\mathcal{M}_{11}^{h}-\mathcal{M}_{22}^{h})}{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} \end{array} $$

(calculated from Eq. (33) in Appendix A.3 Stability of the interior fixed point for general interaction matrices). 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
$$ \begin{aligned} {\mathcal{L}}=&+\left(\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h}\right)\ln p_{1} + \left(\mathcal{M}_{21}^{h}-\mathcal{M}_{11}^{h}\right)\ln (1-p_{1})\\ &-\left(\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p}\right)\ln h_{1} -\left(\mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p}\right)\ln(1-h_{1}) \\ =&+\left(\alpha \gamma (1-\sigma (1-\alpha \kappa))+(1-\alpha) \sigma (1-\alpha \kappa)\right)\ln p_{1} \\&+ \left(\sigma -\alpha \gamma\right)\ln (1-p_{1})+ (1-\alpha \kappa)\sigma \ln h_{1} \\ &+ (1-\alpha (1-\alpha \kappa))\sigma \ln(1-h_{1}). \end{aligned} $$

Due to \(\dot {\mathcal {L}}=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
$$\begin{array}{*{20}l} {\mathcal{L}}_{1} &= b_{h} \ln p_{1} - \sigma p_{1} + d_{p} \ln h_{1} - \sigma h_{1} \end{array} $$
$$\begin{array}{*{20}l} {\mathcal{L}}_{2} &= b_{h} \ln p_{2} - \sigma p_{2} + d_{p} \ln h_{2} - \sigma h_{2}\, . \end{array} $$
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
$$\begin{array}{*{20}l} \frac{h_{1}}{h_{2}} &= \frac{\mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p}}{\mathcal{M}_{11}^{p}-\mathcal{M}_{21}^{p}} \quad\text{and} \end{array} $$
$$\begin{array}{*{20}l} \frac{p_{1}}{p_{2}} &= \frac{\mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h}}{\mathcal{M}_{11}^{h}-\mathcal{M}_{21}^{h}} \end{array} $$
which corresponds to a two-dimensional subspace, then there are two constants that remain invariant over time,
$$\begin{array}{*{20}l} {\mathcal{L}}_{1} &= b_{h} \ln p_{1}+\mathcal{M}_{11}^{h} p_{1}+\mathcal{M}_{12}^{h} p_{2} + d_{p} \ln h_{1} \\ &\quad- \mathcal{M}_{11}^{p} h_{1} - \mathcal{M}_{12}^{p} h_{2} \end{array} $$
$$\begin{array}{*{20}l} {\mathcal{L}}_{2} &= b_{h} \ln p_{2}+\mathcal{M}_{21}^{h} p_{1} +\mathcal{M}_{22}^{h} p_{2}+ d_{p} \ln h_{2} \\ &\quad- \mathcal{M}_{21}^{p} h_{1} - \mathcal{M}_{22}^{p} h_{2}. \end{array} $$

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 [2632]. 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 Stability of the interior fixed point for general interaction matrices).

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, \(\mathcal {M}_{21}^{p}\) 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
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \,\quad\mathcal{P}_{1} & \qquad\qquad\mathcal{P}_{2} \end{array}\\ \mathcal{M}^{h} = \begin{array}{l} \mathcal{H}_{1}\\ \mathcal{H}_{2} \end{array}& \left(\begin{array}{cc} -\sigma & - \alpha (1 - \kappa) \sigma\\ - \alpha \gamma & - \alpha \gamma - (1 - \alpha \kappa) \sigma \end{array} \right) \end{array} \end{array} $$
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \,\qquad\mathcal{H}_{1} & \qquad\quad\mathcal{H}_{2} \end{array}\\ \mathcal{M}^{p} = \begin{array}{l} \mathcal{P}_{1}\\ \mathcal{P}_{2} \end{array}& \left(\begin{array}{cc} {\sigma} & 0\\ \alpha (1 - \kappa) {\sigma} & (1- \alpha \kappa) {\sigma} \end{array} \right). \end{array} \end{array} $$

The analysis in Appendix A.4 Linear interpolation between MA and GfG models 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 Stability of the interior fixed point for general interaction matrices).

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 [3538]. 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 [4346]. 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, 4751], 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, 5255]. 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.


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,p 2)=
$$\begin{array}{@{}rcl@{}} \left(\begin{array}{cccc} b_{h}-p_{1} \sigma -p_{2} \alpha (1-\alpha \kappa) \sigma & 0 & -h_{1} \sigma & -h_{1} \alpha (1-\alpha \kappa) \sigma\\ 0 & b_{h}-p_{1} \alpha \gamma +p_{2} ((1-\alpha \gamma) (1-(1-\alpha \kappa) \sigma)-1) & -h_{2} \alpha \gamma & h_{2} ((1-\alpha \gamma) (1-(1-\alpha \kappa) \sigma)-1) \\ p_{1} \sigma & 0 & h_{1} \sigma -d_{p} & 0 \\ p_{2} \alpha (1-\alpha \kappa) \sigma & p_{2} (1-\alpha \kappa) \sigma & 0 & (h_{2} (1-\alpha \kappa)+h_{1} \alpha (1-\alpha \kappa)) \sigma -d_{p} \\ \end{array} \right)\, . \end{array} $$
At the interior fixed point \((h^{*}_{1},h^{*}_{2},p^{*}_{1},p^{*}_{2})\), we have \(J(h^{*}_{1},h^{*}_{2},p^{*}_{1},p^{*}_{2})=\)
$$ \left(\begin{array}{cccc} 0 & 0 & - d_{p} & d_{p} \alpha (\alpha \kappa -1) \\ 0 & 0 & \frac{d_{p} \alpha \gamma (\alpha (\alpha \kappa -1)+1)}{(\alpha \kappa -1) \sigma} & \frac{d_{p} (\alpha (\alpha \kappa -1)+1) (\alpha \gamma +(\alpha \gamma -1) (\alpha \kappa -1) \sigma)}{(\alpha \kappa -1) \sigma} \\ \frac{b_{h} (\alpha \gamma +(\gamma \alpha +\alpha -1) (\alpha \kappa -1) \sigma)}{\alpha \gamma (\alpha (\alpha \kappa -1)+1)+(\alpha \gamma -1) (\alpha \kappa -1) \sigma} & 0 & 0 & 0 \\ \frac{b_{h} \alpha (1-\alpha \kappa) (\sigma -\alpha \gamma)}{\alpha \gamma (\alpha (\alpha \kappa -1)+1)+(\alpha \gamma -1) (\alpha \kappa -1) \sigma} & \frac{b_{h} (1-\alpha \kappa) (\sigma -\alpha \gamma)}{\alpha \gamma (\alpha (\alpha \kappa -1)+1)+(\alpha \gamma -1) (\alpha \kappa -1) \sigma} & 0 & 0 \\ \end{array} \right)\, . $$
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
$$\begin{array}{*{20}l} \Lambda_{1,2} &= \pm i \sqrt{b_{h} d_{p}} \quad\text{and} \\ \Lambda_{3,4} &= \pm \frac{ \sqrt{b_{h} d_{p}}\sqrt{\sigma (1-\alpha) (1-\alpha \kappa)+\alpha \gamma (1-\sigma(1-\alpha\kappa))} \sqrt{1-\alpha (1-\alpha \kappa)}}{\sqrt{\sigma} \sqrt{\sigma (1-\alpha \gamma) (1-\alpha \kappa)+\alpha \gamma (1-\alpha (1-\alpha \kappa)) }}\sqrt{\alpha \gamma -\sigma }, \end{array} $$

Except the term \(\sqrt {\alpha \gamma -\sigma }\), the remaining factors in in Λ 3,4 are positive. For α γ>σ, allele \(\mathcal {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.2 Stability of the interior fixed point in the Replicator Dynamics

For the system with constant population size, the Jacobian matrix in general is J(h 1,p 1)=
$$ \left(\begin{array}{cc} (1-2 \text{h1}) (\alpha (\gamma -\sigma (1-\text{p1}) (\gamma +(-\gamma \alpha -\alpha +1) \kappa +1))-2 \text{p1} \sigma +\sigma) & \sigma \text{h1} (1-\text{h1}) \left(-\kappa (\gamma +1) \alpha^{2}+(\gamma +\kappa +1) \alpha -2\right) \\ \sigma \text{p1} (1-\text{p1}) (-(1-\alpha) \kappa \alpha -\alpha +2) & \sigma (1-2 \text{p1}) (\text{h1} (-\kappa (1-\alpha) \alpha -\alpha +2)+\alpha \kappa -1) \\ \end{array} \right)\, . $$
At the interior fixed point \((h^{*}_{1},p^{*}_{1})\), the matrix is given by \(J(h^{*}_{1},p^{*}_{1})=\)
$$ \left(\begin{array}{cc} 0 & \frac{(\alpha \kappa -1) \left((\gamma +1) \kappa \alpha^{2}-(\gamma +\kappa +1) \alpha +2\right) (\alpha (\alpha \kappa -1)+1) \sigma }{((\alpha -1) \kappa \alpha -\alpha +2)^{2}} \\ -\frac{((\alpha -1) \kappa \alpha -\alpha +2) (\alpha \gamma -\sigma) (\alpha \gamma +(\gamma \alpha +\alpha -1) (\alpha \kappa -1) \sigma)}{\left((\gamma +1) \kappa \alpha^{2}-(\gamma +\kappa +1) \alpha +2\right)^{2} \sigma} & 0 \\ \end{array} \right)\, . $$
The eigenvalues are
$$ \Lambda_{1,2} = \mp i \frac{\sqrt{(1-\alpha \kappa) (1-\alpha (1-\alpha \kappa)) (\alpha \gamma (1-\sigma (1-\alpha \kappa))+(1-\alpha) \sigma (1-\alpha \kappa))}}{\sqrt{(1+(1-\alpha) (1-\alpha \kappa)) (1-\alpha \gamma (1-\alpha \kappa)+(1-\alpha) (1-\alpha \kappa))}}\sqrt{\sigma -\alpha \gamma} \, . $$

For α γ<σ, the eigenvalues are purely imaginary, hence, the fixed point is a neutral center.

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,
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \quad\mathcal{P}_{1} & \quad\mathcal{P}_{2} \end{array}\\ \mathcal{M}^{h} = \begin{array}{l} \mathcal{H}_{1}\\ \mathcal{H}_{2} \end{array}& \left(\begin{array}{cc} \mathcal{M}_{11}^{h} & \mathcal{M}_{12}^{h}\\ \mathcal{M}_{21}^{h} & \mathcal{M}_{22}^{h} \end{array} \right) \end{array} \end{array} $$
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \quad\mathcal{H}_{1} & \!\quad\mathcal{H}_{2} \end{array}\\ \mathcal{M}^{p} = \begin{array}{l} \mathcal{P}_{1}\\ \mathcal{P}_{2} \end{array}& \left(\begin{array}{cc} \mathcal{M}_{11}^{p} & \mathcal{M}_{12}^{p}\\ \mathcal{M}_{21}^{p} &\mathcal{M}_{22}^{p} \end{array} \right)\, . \end{array} \end{array} $$
The interior fixed point for our host parasite system with Lotka-Volterra dynamics (Eq. (2)) is then
$$\begin{array}{*{20}l} h_{1}^{*} &= \frac{\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p}}{\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}} d_{p}\\ &\\ h_{2}^{*} &=\frac{\mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p}}{\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}} d_{p} \end{array} $$
$$\begin{array}{*{20}l} p_{1}^{*} &= \frac{\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h}}{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} b_{h} \\ &\\ p_{2}^{*} &= \frac{\mathcal{M}_{21}^{h}-\mathcal{M}_{11}^{h}}{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} b_{h} \, . \end{array} $$
The Jacobian matrix at any defined point is J(h 1,h 2,p 1,p 2)=
$$ \left(\begin{array}{cccc} {b_{h}}+\mathcal{M}_{11}^{h} {p_{1}}+\mathcal{M}_{12}^{h} {p_{2}} & 0 & {h_{1}} \mathcal{M}_{11}^{h} & {h_{1}} \mathcal{M}_{12}^{h} \\ 0 & {b_{h}}+\mathcal{M}_{21}^{h} {p_{1}}+\mathcal{M}_{22}^{h} {p_{2}} & {h_{2}} \mathcal{M}_{21}^{h} & {h_{2}} \mathcal{M}_{22}^{h} \\ \mathcal{M}_{11}^{p} {p_{1}} & \mathcal{M}_{12}^{h} {p_{1}} & -{d_{p}}+{h_{1}} \mathcal{M}_{11}^{p}+{h_{2}} \mathcal{M}_{12}^{p} & 0 \\ \mathcal{M}_{21}^{p} {p_{2}} & \mathcal{M}_{22}^{p} {p_{2}} & 0 & -{d_{p}}+{h_{1}} \mathcal{M}_{21}^{p}+{h_{2}} \mathcal{M}_{22}^{p} \\ \end{array} \right)\, . $$
At the interior fixed point \((h^{*}_{1},h^{*}_{2},p^{*}_{1},p^{*}_{2})\), we now have \(J(h^{*}_{1},h^{*}_{2},p^{*}_{1},p^{*}_{2})=\)
$$ \left(\begin{array}{cccc} 0 & 0 & \frac{\mathcal{M}_{11}^{h} (\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p}){d_{p}} }{\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}} & \frac{\mathcal{M}_{12}^{h} (\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p}){d_{p}}}{\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}-\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}} \\ 0 & 0 & \frac{\mathcal{M}_{21}^{h} (\mathcal{M}_{11}^{p}-\mathcal{M}_{21}^{p}){d_{p}} }{\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}} & \frac{\mathcal{M}_{22}^{h} (\mathcal{M}_{11}^{p}-\mathcal{M}_{21}^{p}){d_{p}}}{\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p}} \\ \frac{\mathcal{M}_{11}^{p}(\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h}){b_{h}} }{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} & \frac{\mathcal{M}_{12}^{p} (\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h}) {b_{h}}}{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} & 0 & 0 \\ \frac{\mathcal{M}_{21}^{p}(\mathcal{M}_{21}^{h}-\mathcal{M}_{11}^{h}){b_{h}} }{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} & \frac{\mathcal{M}_{22}^{p} (\mathcal{M}_{21}^{h}-\mathcal{M}_{11}^{h}) {b_{h}}}{\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h}} & 0 & 0 \\ \end{array} \right)\, . $$
There are four eigenvalues
$$\begin{array}{*{20}l} \Lambda_{1,2} &= \pm i \sqrt{b_{h} d_{p}} \quad\text{and} \\ \Lambda_{3,4} &= \pm i \frac{\sqrt{b_{h} d_{p}} \sqrt{(\mathcal{M}_{11}^{h}-\mathcal{M}_{21}^{h}) (\mathcal{M}_{12}^{h}-\mathcal{M}_{22}^{h}) (\mathcal{M}_{11}^{p}-\mathcal{M}_{21}^{p}) (\mathcal{M}_{12}^{p}-\mathcal{M}_{22}^{p})}}{\sqrt{ (\mathcal{M}_{11}^{h} \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} \mathcal{M}_{21}^{h})} \sqrt{(\mathcal{M}_{11}^{p} \mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p} \mathcal{M}_{21}^{p})}}\, . \end{array} $$
It is often assumed that (i) \(\mathcal {M}_{11}^{h}<\mathcal {M}_{21}^{h}\leq 0\) (\(\mathcal {H}_{2}\) is beneficial if there is only \(\mathcal {P}_{1}\) in the population), (ii) \(\mathcal {M}_{22}^{h}<\mathcal {M}_{12}^{h}\leq 0\) (\(\mathcal {H}_{1}\) is beneficial if there is only \(\mathcal {P}_{2}\) in the population), (iii) \(\mathcal {M}_{11}^{p}>\mathcal {M}_{21}^{p}\geq 0\) (\(\mathcal {P}_{1}\) is beneficial if there is only \(\mathcal {H}_{1}\) in the population), and (iv) \(\mathcal {M}_{22}^{p}>\mathcal {M}_{12}^{p}\geq 0\) (\(\mathcal {P}_{1}\) is beneficial if there is only \(\mathcal {H}_{1}\) in the population). With these minimal assumptions the eigenvalues are purely imaginary, i.e., the interior fixed point is a neutrally stable center. The ratio between the eigenvalues, which determines the oscillation frequencies at the center, differs in different interaction models. For example, in Parker [13] the interaction matrices for haploid types are
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \qquad \; \quad\mathcal{P}_{1} & \qquad\qquad\qquad\mathcal{P}_{2} \end{array}\\ \mathcal{M}^{h} = \begin{array}{l} \mathcal{H}_{1}\\ \mathcal{H}_{2} \end{array}& \left(\begin{array}{cc} -\sigma & -(1-\kappa) \sigma\\ -\gamma -\sigma (1-\tau) & -\sigma (\alpha -\kappa +1)-\gamma \end{array} \right) \end{array} \end{array} $$
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \qquad\mathcal{H}_{1} & \;\quad\qquad\mathcal{H}_{2} \end{array}\\ \mathcal{M}_{i,j}^{p} = \begin{array}{l} \mathcal{P}_{1}\\ \mathcal{P}_{2} \end{array}& \left(\begin{array}{cc} \sigma & \sigma (1-\tau)\\ (1-\kappa) \sigma & \sigma (1+\alpha -\kappa) \end{array} \right)\, , \end{array} \end{array} $$
where the notations a, c, k, t, and s in [13] are changed to α, γ, κ, τ, and σ, respectively. According to [13], the fitness of the “narrowly virulent pathogen” \(\mathcal {P}_{1}\) is reduced by a factor τ by interacting with the resistant host \(\mathcal {H}_{2}\); a fitness penalty κ (the cost of virulence) is inflicted on the “broadly virulent pathogen” \(\mathcal {P}_{2}\) independent of which host it exploits; α the “advantage of adapted pathogens on resistant host” measures a special advantage of \(\mathcal {P}_{2}\) on \(\mathcal {H}_{2}\); a fitness penalty γ (the cost of resistance) is paid by the resistant host \(\mathcal {H}_{2}\). When τ=κ=α=1 and γ=0 the fitnesses conform to the pattern of pure MA model. When τ=1 and α=0 the fitnesses revert to a pure GfG pattern. The ratio between the two oscillation frequencies at the interior fixed point is
$$ m= \frac{\sqrt{\kappa} \sqrt{\alpha \sigma +\gamma} \sqrt{\alpha -\kappa +\tau} \sqrt{\sigma \tau -\gamma }}{\sqrt{\sigma} \sqrt{\alpha -\kappa \tau +\tau} \sqrt{\sigma (\alpha -\kappa \tau +\tau)+\gamma \kappa }}\, . $$

The ratio is 1 for pure MA model. With a set of parameter used in [13], α=0.33, γ=0, κ=0.05, and σ=τ=1 the ratio is about 0.1.

The same method can be applied for the system with constant population size. In that case, the interior fixed point expressed by the general interaction matrices elements is
$$\begin{array}{*{20}l} h_{1}^{*} &= \frac{\mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p}}{\mathcal{M}_{11}^{p}+ \mathcal{M}_{22}^{p}-\mathcal{M}_{12}^{p} - \mathcal{M}_{21}^{p}} \end{array} $$
$$\begin{array}{*{20}l} p_{1}^{*} &= \frac{\mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h}}{\mathcal{M}_{11}^{h}+ \mathcal{M}_{22}^{h}-\mathcal{M}_{12}^{h} -\mathcal{M}_{21}^{h}}\, , \end{array} $$
while \(h_{2}^{*}=1-h_{1}^{*}\) and \(p_{2}^{*}=1-p_{1}^{*}\). The eigenvalues of the Jacobian matrix at the interior fixed point are
$$\begin{array}{*{20}l} \Lambda_{1,2} &= 0 \quad\text{and} \\ \Lambda_{3,4} &= \pm i \sqrt{-\frac{(\mathcal{M}_{11}^{h} - \mathcal{M}_{21}^{h}) (\mathcal{M}_{22}^{h} - \mathcal{M}_{12}^{h}) (\mathcal{M}_{11}^{p} - \mathcal{M}_{21}^{p}) (\mathcal{M}_{22}^{p} - \mathcal{M}_{12}^{p})}{(\mathcal{M}_{11}^{h} + \mathcal{M}_{22}^{h}- \mathcal{M}_{12}^{h} -\mathcal{M}_{21}^{h}) (\mathcal{M}_{11}^{p} + \mathcal{M}_{22}^{p} - \mathcal{M}_{12}^{p} - \mathcal{M}_{21}^{p})}} \end{array} $$

Hence, there only is one oscillation frequency at the interior fixed point in models with constant population size, regardless of the specific assumption for the interaction matrices.

A.4 Linear interpolation between MA and GfG models

Alternatively to the models of [14] and [13], one could also use a linear interpolation between MA and gene-for-gene model, where the matrix elements linearly spans over the values of the two models as a single parameter α varies between 0 and 1
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \,\quad\mathcal{P}_{1} & \qquad\qquad\mathcal{P}_{2} \end{array}\\ \mathcal{M}^{h} = \begin{array}{l} \mathcal{H}_{1}\\ \mathcal{H}_{2} \end{array}& \left(\begin{array}{cc} -\sigma & - \alpha (1 - \kappa) \sigma\\ - \alpha \gamma & - \alpha \gamma - (1 - \alpha \kappa) \sigma \end{array} \right) \end{array} \end{array} $$
$$\begin{array}{*{20}l} \begin{array}{ll} &\begin{array}{cc} \;\qquad\mathcal{H}_{1} & \quad\qquad\mathcal{H}_{2} \end{array}\\ \mathcal{M}^{p} = \begin{array}{l} \mathcal{P}_{1}\\ \mathcal{P}_{2} \end{array}& \left(\begin{array}{cc} {\sigma} & 0 \\ \alpha (1 - \kappa) {\sigma} & (1- \alpha \kappa) {\sigma} \end{array} \right)\, . \end{array} \end{array} $$
The fixed point with Lotka-Volterra dynamics is then
$$\begin{array}{*{20}l} h_{1}^{*} &= \frac{1}{\sigma} d_{p}\\ \\ h_{2}^{*} &=\frac{1-\alpha (1-\kappa)}{\sigma (1-\alpha \kappa)} d_{p} \end{array} $$
$$\begin{array}{*{20}l} p_{1}^{*} &=\frac{\alpha (\gamma -\sigma)+\sigma }{\sigma (\alpha \gamma (1-\alpha (1-\kappa))+\sigma(1-\alpha \kappa))} b_{h} \\ \\ p_{2}^{*} &=\frac{\sigma -\alpha \gamma }{\sigma (\alpha \gamma (1-\alpha (1-\kappa))-\sigma (1-\alpha \kappa))} b_{h}\, , \end{array} $$
and the eigenvalues of the Jacobian matrix at this point are
$$ {\fontsize{9}{6} \begin{aligned} \Lambda_{1,2} &= \pm i \sqrt{b_{h} d_{p}} \quad\text{and} \\ \Lambda_{3,4} &= \!\pm i \sqrt{b_{h} d_{p}} \frac{\sqrt{1-\alpha (1-\kappa)} \sqrt{(\sigma -\alpha \gamma) (\alpha \gamma +(1-\alpha) \sigma)}}{\sqrt{\sigma} \sqrt{\alpha \gamma (1-\alpha (1-\kappa))+\sigma (1-\alpha \kappa)}}\, . \end{aligned}} $$
As long as α γ<σ, the ratio m=Λ 3,4/Λ 1,2 increases with increasing cost of virulence κ, while m decreases with increasing α. For α1 we find
$$ m\approx 1-\frac{\alpha (\gamma +2 (1-\kappa) \sigma)}{2 \sigma }\, . $$

Hence, there are always two distinct oscillation frequencies at the interior fixed point in gene-for-gene-like models with changing population size.

With Replicator Dynamics, the interior fixed point is
$$\begin{array}{*{20}l} h_{1}^{*} &= \frac{1-\alpha\kappa}{2-\alpha} \end{array} $$
$$\begin{array}{*{20}l} p_{1}^{*} &= \frac{\alpha\gamma+\sigma(1-\alpha)}{\sigma(2-\alpha)}\, , \end{array} $$
while \(h_{2}^{*}=1-h_{1}^{*}\) and \(p_{2}^{*}=1-p_{1}^{*}\). The eigenvalues of the Jacobian matrix at the interior fixed point are
$$ {\fontsize{9}{6} \begin{aligned} \Lambda_{1,2} &= 0 \quad\text{and} \\ \Lambda_{3,4} &= \pm i \frac{\sqrt{1-\alpha+\alpha^{2} \kappa (1-\kappa)} \sqrt{(\sigma -\alpha \gamma) (\alpha \gamma +(1-\alpha) \sigma)}}{2-\alpha } \end{aligned}} $$
Hence, there is only one oscillation frequency l=Λ 3/(i2π) at the interior fixed point in models with constant population size. As long as α γ<σ, the oscillation frequency l decreases with α and increases with γ and σ, while l increases with κ until κ reaches the value 1/2, then l decreases as κ increases from 1/2 to 1. For α1,
$$ l\approx \frac{1}{2\pi} \left(\frac{\sigma}{2}-\frac{\alpha\sigma}{4} \right)\, . $$



YS, CSG, and AT acknowledge generous funding by the Max Planck Society. AP was funded through grants of the German Science Foundation to HS (DFG grants SCHU 1415/8 and SCHU 1415/9 within the German priority program SPP1399 on host-parasite coevolution). AP was additionally supported by the International Max-Planck Research School (IMPRS) for Evolutionary Biology.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

Max Planck Institute for Evolutionary Biology
New Zealand Institute for Advanced Study, Massey University
Department of Evolutionary Ecology and Genetics, University of Kiel


  1. Woolhouse MEJ, Webster JP, Domingo E, Charlesworth B, Levin BR. Biological and biomedical implications of the co-evolution of pathogens and their hosts. Nat Genet. 2002; 32(4):569–77.View ArticlePubMedGoogle Scholar
  2. Woolhouse MEJ, Haydon DT, Antia R. Emerging pathogens: the epidemiology and evolution of species jumps. Trends Ecol Evol. 2005; 20(5):238–44. doi:10.1016/j.tree.2005.02.009.View ArticlePubMedGoogle Scholar
  3. Van der Plank JE. Disease Resistance in Plants, 2nd revised edition edition. Orlando: Academic Press Inc; 1984.Google Scholar
  4. Gladieux P, Byrnes EJ, Aguileta G, C Fisher M, Heitman J, Giraud T. Epidemiology and evolution of fungal pathogens in plants and animals. In: Genetics and Evolution of Infectious Disease. London: Elsevier: 2011l. p. 59–132. Scholar
  5. Altizer S, Harvell D, Friedle E. Rapid evolutionary dynamics and disease threats to biodiversity. Trends Ecol Evol. 2003; 18(11):589–96. doi:10.1016/j.tree.2003.08.013.View ArticleGoogle Scholar
  6. Thompson RCA, Lymbery AJ, Smith A. Parasites, emerging disease and wildlife conservation. Int J Parasitol. 2010; 40(10):1163–70. doi:10.1016/j.ijpara.2010.04.009.View ArticlePubMedGoogle Scholar
  7. Flor HH. The complementary genetic systems in flax and flax rust. Adv Genet. 1956; 8:29–54.View ArticleGoogle Scholar
  8. Jones JDG, Dangl JL. The plant immune system. Nature. 2006; 444(7117):323–9.View ArticlePubMedGoogle Scholar
  9. Grosberg RK, Hart MW. Mate selection and the evolution of highly polymorphic self/nonself recognition genes. Science. 2000; 289:2111–4.View ArticlePubMedGoogle Scholar
  10. Lively CM. A review of red queen models for the persistence of obligate sexual reproduction. J Heredity. 2010; 101(suppl 1):13–20. doi:10.1093/jhered/esq010. ArticleGoogle Scholar
  11. Lively CM, Apanius V. Genetic diversity in host-parasite interactions. In: Ecology of Infectious Diseases in Natural Populations. Cambridge, UK: Cambridge University Press: 1995. p. 421–49.Google Scholar
  12. Leonard KJ. Selection pressures and plant pathogens. Ann NY Acad Sci. 1977; 287:207–22. doi:10.1111/j.1749-6632.1977.tb34240.x.View ArticleGoogle Scholar
  13. Parker MA. Pathogens and sex in plants. Evol Ecol. 1994; 8(5):560–84. doi:10.1007/BF01238258.View ArticleGoogle Scholar
  14. Agrawal A, Lively CM. Infection genetics: gene-for-gene versus matching-alleles models and all points in between. Evol Ecol Res. 2002; 4:79–90.Google Scholar
  15. Agrawal AF, Lively CM. Modelling infection as a two-step process combining gene-for-gene and matching-allele genetics. Proc R Soc B: Biol Sci. 2003; 270(1512):323–34.View ArticleGoogle Scholar
  16. Tellier A, Brown JKM. Polymorphism in multilocus host-paraiste coevolutionary interactions. Genetics. 2007; 177:1777–90.PubMed CentralView ArticlePubMedGoogle Scholar
  17. Sardanyés J, Solé RV. Matching allele dynamics and coevolution in a minimal predator–prey replicator model. Phys Lett A. 2008; 372:341–50.View ArticleGoogle Scholar
  18. Zeeman EC. Population dynamics from game theory. Lecture Notes Math. 1980; 819:471–97.View ArticleGoogle Scholar
  19. Taylor PD, Jonker L. Evolutionarily stable strategies and game dynamics. Math Biosci. 1978; 40:145–56.View ArticleGoogle Scholar
  20. Hofbauer J, Schuster P, Sigmund K. A note on evolutionary stable strategies and game dynamics. J Theor Biol. 1979; 81:609–12.View ArticlePubMedGoogle Scholar
  21. Hofbauer J, Sigmund K. Evolutionary Games and Population Dynamics. Cambridge, UK: Cambridge University Press; 1998.View ArticleGoogle Scholar
  22. Traulsen A, Claussen JC, Hauert C. Coevolutionary dynamics: From finite to infinite populations. Phys Rev Lett. 2005; 95:238701.View ArticlePubMedGoogle Scholar
  23. Leonard KJ. Stability of equilibria in a gene-for-gene coevolution model of host-parasite interactions. Phytopathology. 1994; 84:70–7.View ArticleGoogle Scholar
  24. Schuster P, Sigmund K. Replicator dynamics. J Theor Biol. 1983; 100:533–8.View ArticleGoogle Scholar
  25. Strogatz S. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Cambridge, Massachusetts: Perseus Books; 1994.Google Scholar
  26. Frank SA. Specificity versus detectable polymorphism in host–parasite genetics. Proc R Soc London. Series B: Biol Sci. 1993; 254(1341):191–7. doi:10.1098/rspb.1993.0145. ArticleGoogle Scholar
  27. Otto SP, Michalakis Y. The evolution of recombination in changing environments. Trends Ecol Evol. 1998; 13(4):145–51. doi:10.1016/S0169-5347(97)01260-3.View ArticlePubMedGoogle Scholar
  28. Lively CM. The maintenance of sex: host–parasite coevolution with density-dependent virulence. J Evol Biol. 2009; 22(10):2086–93. doi:10.1111/j.1420-9101.2009.01824.x.View ArticlePubMedGoogle Scholar
  29. Gokhale CS, Papkou A, Traulsen A, Schulenburg H. Lotka-Volterra dynamics kills the Red Queen: population size fluctuations and associated stochasticity dramatically change host-parasite coevolution. BMC Evol Biol. 2013; 13:254.PubMed CentralView ArticlePubMedGoogle Scholar
  30. Luijckx P, Fienberg H, Duneau D, Ebert D. A Matching-Allele Model Explains Host Resistance to Parasites. Current Biol. 2013; 23(12):1085–8.View ArticleGoogle Scholar
  31. Clay K, Kover PX. The red queen hypothesis and plant/pathogen interactions. Ann Rev Phytopathol. 1996; 34(1):29–50. doi:10.1146/annurev.phyto.34.1.29. PMID: 15012533. arXiv: View ArticleGoogle Scholar
  32. Brown JKM, Tellier A. Plant-parasite coevolution: Bridging the gap between genetics and ecology. Ann Rev Phytopathol. 2011; 49(1):345–67. doi:10.1146/annurev-phyto-072910-095301. PMID: 21513455. arXiv: ArticleGoogle Scholar
  33. Bergelson J, Dwyer G, Emerson JJ. Models and data on plant-enemy coevolution. Ann Rev Genet. 2001; 35(1):469–99. doi:10.1146/annurev.genet.35.102401.090954.View ArticlePubMedGoogle Scholar
  34. Tellier A, Brown JKM. Stability of genetic polymorphism in host-parasite interactions. Proc R Soc B: Biol Sci. 2007; 274(1611):809–17. doi:10.1098/rspb.2006.0281.View ArticleGoogle Scholar
  35. Hendry AP, Kinnison MT. Perspective: The pace of modern life: Measuring rates of contemporary microevolution. Evolution. 1999; 53(6):1637. doi:10.2307/2640428.View ArticleGoogle Scholar
  36. Thompson JN. Rapid evolution as an ecological process. Trends Ecol Evol. 1998; 13(8):329–32. doi:10.1016/S0169-5347(98)01378-0.View ArticlePubMedGoogle Scholar
  37. Hairston NG, Ellner SP, Geber MA, Yoshida T, Fox JA. Rapid evolution and the convergence of ecological and evolutionary time. Ecol Lett. 2005; 8(10):1114–27. doi:10.1111/j.1461-0248.2005.00812.x.View ArticleGoogle Scholar
  38. Schoener TW. The newest synthesis: Understanding the interplay of evolutionary and ecological dynamics. Science. 2011; 331(6016):426–9. doi:10.1126/science.1193954.View ArticlePubMedGoogle Scholar
  39. Dieckmann U. Adaptive dynamics of pathogen-host interactions. In: Adaptive Dynamics of Infectious Diseases: in Pursuit of Virulence Management. Cambridge Studies in Adaptive Dynamics. Cambridge, UK: Cambridge University Press: 2002. p. 39–59. doi:10.1017/CBO9780511525728.006.View ArticleGoogle Scholar
  40. Day T. Modelling the ecological context of evolutionary change: Déjàvu or something new? In: Beisner K, Cuddington BE, editors. Ecological Paradigms Lost. Theoretical Ecology Series. Chap. 13 - Modelling the ecological context of evolutionary change. Burlington: Academic Press: 2005. p. 273–309. Scholar
  41. Fenner F, Fantini B. Biological Control of Vertebrate Pests. The History of Myxomatosis–an Experiment in Evolution. Oxfordshire: CABI Publishing; 1999.Google Scholar
  42. Gomulkiewicz R, Holt RD. When does evolution by natural selection prevent extinction? Evolution. 1995; 49(1):201. doi:10.2307/2410305.View ArticleGoogle Scholar
  43. O’Brien SJ, Evermann JF. Interactive influence of infectious disease and genetic diversity in natural populations. Trends Ecol Evol. 1988; 3:254–9.View ArticlePubMedGoogle Scholar
  44. Lande R. Genetics and demography in biological conservation. Science. 1988; 241(4872):1455–60.View ArticlePubMedGoogle Scholar
  45. Gomulkiewicz R, Houle D. Demographic and genetic constraints on evolution. Am Naturalist. 2009; 174(6):218–29. doi:10.1086/599011.View ArticleGoogle Scholar
  46. Saccheri I, Hanski I. Natural selection and population dynamics. Trends Ecol Evol. 2006; 21(6):341–7. doi:10.1016/j.tree.2006.03.018.View ArticlePubMedGoogle Scholar
  47. Frank SA. Ecological and genetic models of host-pathogen coevolution. Heredity. 1991; 67(1):73–83.View ArticlePubMedGoogle Scholar
  48. Frank SA. Coevolutionary genetics of plants and pathogens. Evol Ecol. 1993; 7(1):45–75.View ArticleGoogle Scholar
  49. Gandon S, Capowiez Y, Dubois Y, Michalakis Y, Olivieri I. Local adaptation and gene-for-gene coevolution in a metapopulation model. Proc R Soc B: Biol Sci. 1996; 263(1373):1003–1009. doi:10.1098/rspb.1996.0148.View ArticleGoogle Scholar
  50. Quigley BJZ, García López D, Buckling A, McKane AJ, Brown SP. The mode of host-parasite interaction shapes coevolutionary dynamics and the fate of host cooperation. Proc R Soc B: Biol Sci. 2012; 279(1743):3742–748. doi:10.1098/rspb.2012.0769.View ArticleGoogle Scholar
  51. Ashby B, Gupta S. Parasitic castration promotes coevolutionary cycling but also imposes a cost on sex. Evolution. 2014; 68(8):2234–44. doi:10.1111/evo.12425.PubMedGoogle Scholar
  52. Thompson JN, Burdon JJ. Gene-for-gene coevolution between plants and parasites. Nature. 1992; 360(6400):121–5. doi:10.1038/360121a0.View ArticleGoogle Scholar
  53. Carius HJ, Little TJ, Ebert D. Genetic variation in a host-parasite association: potential for coevolution and frequency-dependent selection. Evolution. 2001; 55(6):1136–45.View ArticlePubMedGoogle Scholar
  54. Wilfert L, Jiggins FM. Host-parasite coevolution: genetic variation in a virus population and the interaction with a host gene. J Evol Biol. 2010; 23(7):1447–55. doi:10.1111/j.1420-9101.2010.02002.x.View ArticlePubMedGoogle Scholar
  55. Luijckx P, Fienberg H, Duneau D, Ebert D. Resistance to a bacterial parasite in the crustacean daphnia magna shows mendelian segregation with dominance. Heredity. 2012; 108(5):547–51. doi:10.1038/hdy.2011.122.PubMed CentralView ArticlePubMedGoogle Scholar
  56. Jayakar SD. A mathematical model for interaction of gene frequencies in a parasite and its host. Theoretical Popul Biol. 1970; 1(2):140–64.View ArticleGoogle Scholar
  57. May RM, Anderson RM. Epidemiology and genetics in the coevolution of parasites and hosts. Proc R Soc B: Biol Sci. 1983; 219(1216):281–313. doi:10.1098/rspb.1983.0075.View ArticleGoogle Scholar
  58. Thrall PH, Burdon JJ. Host-pathogen dynamics in a metapopulation context: The ecological and evolutionary consequences of being spatial. J Ecol. 1997; 85(6):743–53. doi:10.2307/2960598.View ArticleGoogle Scholar
  59. Thrall PH, Burdon JJ. Evolution of gene-for-gene systems in metapopulations: the effect of spatial scale of host and pathogen dispersal. Plant Pathol. 2002; 51(2):169–84. doi:10.2307/2960598.View ArticleGoogle Scholar
  60. Salathé M, Scherer A, Bonhoeffer S. Neutral drift and polymorphism in gene-for-gene systems. Ecol Lett. 2005; 8:925–32.View ArticleGoogle Scholar
  61. Karasov TL, Kniskern JM, Gao L, DeYoung BJ, Ding J, Dubiella U, Lastra RO, Nallu S, Roux F, Innes RW, Barrett LG, Hudson RR, Bergelson J. The long-term maintenance of a resistance polymorphism through diffuse interactions. Nature. 2014; advance online publication doi:10.1038/nature13439.
  62. Sasaki A. Host-parasite coevolution in a multilocus gene-for-gene system. Proc R Soc B: Biol Sci. 2000; 267(1458):2183–188. doi:10.1098/rspb.2000.1267.View ArticleGoogle Scholar
  63. Hauert C, Holmes M, Doebeli M. Evolutionary games and population dynamics: maintenance of cooperation in public goods games. Proc R Soc B. 2006; 273:2565–70.PubMed CentralView ArticlePubMedGoogle Scholar
  64. García-Arenal F, Fraile A. Trade-offs in host range evolution of plant viruses. Plant Pathol. 2013; 62:2–9. doi:10.1111/ppa.12104.View ArticleGoogle Scholar


© Song et al. 2015