## Abstract

Host condition is often likely to influence parasite virulence. Furthermore, condition may often be correlated with host density, and therefore, it is important to understand the role of density-dependent virulence (DDV). We examine the consequences of DDV to the evolution of parasites in both seasonal and non-seasonal environments. In particular, we consider seasonality in host birth rate that results in a fluctuating host density and therefore a variable virulence. We show that parasites are selected for lower exploitation, and therefore lower transmission and virulence as the strength of DDV increases without seasonality. This is an important insight from our models; DDV has the opposite effect on the evolution of parasites to that of higher baseline mortality. Our key result is that although seasonality does not affect the evolution of virulence in classical models, with DDV parasites in seasonal environments are predicted to evolve to be more acute. This suggests that in more seasonal environments wildlife disease is likely to be more rather than less virulent if DDV is widespread.

## 1. Introduction

Parasites by definition cause host damage, and understanding how parasite virulence evolves is key to disease management [1,2]. As a consequence, there is a substantial theoretical literature that has explored the effects of epidemiological characteristics on the evolution of virulence, defined in this literature as parasite-induced mortality [1,3–5]. Recent empirical work has emphasized that virulence is likely to be context dependent, with the lethality of infection contingent on factors such as host nutritional status, oxygen availability or temperature [6–8]. Clearly, a key factor in determining this context is host density leading to the potential for density-dependent virulence (DDV).

Host density can influence disease incidence and severity in a number of ways [9–11]. Poor condition can alter the likelihood of transmission [9], however, the effect is system specific if present at all [10]. The field survey of Lively *et al.* [12] in *Impatiens capensis* infected with the rust *Puccinia recondita*, for example, revealed that the proportion of plants that were infected was not related to density, but that the effect of infection on plant growth was more severe under high host density (i.e. DDV was present). Furthermore, a recent *in situ* experiment [13] reported DDV in oomycete-infected seedlings of the Neotropical tree, *Sebastiana longicuspis*; and in general, DDV seems to be widespread in plant disease [10,11]. In animals, several papers report findings that relate host starvation to disease-induced mortality [7,14,15]. As a whole, this work emphasizes the importance of condition to virulence in a diverse range of systems.

Despite being increasingly recognized in empirical studies, the impact of DDV on disease dynamics has received relatively little theoretical attention. Lively [16] developed a population dynamic model with density-dependent regulation of birth where the density-dependent component was amplified due to infection. His results showed that an infection appearing benign at low density can be effectively castrating at the higher equilibrium density. This emphasizes the importance of considering parasite evolution in its ecological context, which will often include population cycles. Seasonality, a cause of such cycles, is well known to have a profound impact on human and wildlife disease epidemiology [17–20] but has rarely been considered in evolutionary models. The evolution of the sensitivity to seasonality [21] and how seasonality may contribute to the evolution of the parasite life cycle have been examined [22]. Since seasonality may lead to density fluctuations, an exploration of the role of DDV in infectious disease systems subject to fluctuating host abundance will provide novel insight into the evolution of virulence. We therefore develop evolutionary models where DDV acts to increase disease-induced mortality as host density increases and where hosts are subject to a fluctuating birth rate. We apply modern evolutionary game theory (adaptive dynamics; [23,24]) to assess the evolutionary epidemiology of this system as the amplitude of seasonality increases.

## 2. Methods and results

### (a) Model framework

We develop a model for the density of susceptibles, *X*, infecteds, *Y* and recovereds (immune), *Z*, in which immunity can wane (an SIRS framework). This is represented by the following nonlinear ordinary differential equations,
2.1
2.2
2.3where the total host density, *H*, is the sum of the densities of susceptibles, infecteds and recovereds (*H* = *X* + *Y* + *Z*). All hosts reproduce at rate *a*, with host self-regulation through a crowding parameter *q*, which is related to carrying capacity, *K*, as *K* = (*a* − *b*)/*aq*. All offspring are born susceptible. Hosts die at natural death rate *b*. Transmission is a mass action process between susceptible and infected types, with transmission coefficient *β*. High host density, associated with increased scarcity of resources, can lead to increases in disease-induced mortality [12] and therefore, we assume infected hosts are subject to both some baseline virulence rate, *α*, and a density-dependent component, *c**α**H*, where *c* scales the impact of density on virulence. Infected hosts recover to immunity at rate *γ*, while recovered hosts lose immunity at rate *μ*.

The parasite is maintained at endemic levels when the host-only equilibrium () = ((*a* − *b*)/*aq*, 0, 0) becomes unstable. Analysis of the eigenvalues shows that the host-only equilibrium loses stability when,
2.4therefore, increasing the strength of DDV makes it harder for the parasite to invade, requiring a greater transmission rate.

### (b) Parasite evolution

In an evolutionary context, adaptation of a virulence trait is likely to be based on its positive correlation with disease transmission particularly for obligate parasites (see [25,26] for experimental support and [4,27,28] for theoretical application). The mechanism underlying this correlation assumes that an increase in parasite replication rate will enhance transmission but also lead to host damage. We therefore assume that an increase in transmission, *β*, in the parasite is associated with an increase in baseline virulence, *α* (i.e. a trade-off between *β* and *α* such that *α* = *α*(*β*)) and that total virulence is amplified as density increases through DDV. Evolutionary models based on the maximization of an infection's basic reproduction number, , may not be accurate when density dependence is present in epidemiological terms [2]. The recent framework of adaptive dynamics [23,24] which explicitly relates ecological dynamics to evolutionary dynamics is more appropriate since we are investigating the implication of DDV, a density-dependent epidemiological term at the ecological scale, to the evolution of parasites. This method assumes a separation of evolutionary and ecological time scales—an assumption not always valid for epidemiological systems [29]. However, simulations where mutations are allowed to occur before the ecological attractor has been reached, have often confirmed the qualitative robustness of results to a relaxation of this assumption [30]. Taking this approach, we assume a rare mutant with a slightly different phenotypic value of the adaptive trait (here virulence) attempts to invade a resident population at equilibrium.

The success of the mutant depends on its invasion fitness, given by,
2.5where denotes the mutant transmission rate and *X*, *H* are the equilibrium densities at the endemic steady state of the SIRS system, which depends on the resident trait *β*. If *r* > 0 then the mutant parasite will invade to coexist with or replace the resident. Through a series of mutation–substitution events, the population will evolve in the direction of the fitness gradient until it reaches the vicinity of the singularity, *β**, where the fitness gradient is infinitesimal.

Singularities therefore satisfy,
2.6The evolutionary outcome for a singularity depends on two criteria: evolutionary stability (ES, whether the strategy is a local fitness maximum or minimum), requiring for a fitness maximum that,
2.7and convergence stability (CS, whether the strategy is locally attracting or repelling) [23], requiring for a local attractor that,
2.8
2.9To examine how the singular transmission rate, *β**, varies for changes in other parameters (and therefore how the singular value of baseline virulence, *α**, varies) we choose a trade-off such that the singular point is a continuously stable strategy (CSS), i.e. ES and CS, which is an attracting endpoint for evolutionary trajectories (see figure 1 legend for the exact form of the trade-off). By equation (2.7) this requires *α*″(*β*) > 0, corresponding to accelerating costs of transmission.

Figure 1*a* shows the behaviour as the strength of DDV, *c*, is increased, indicating that the parasite decreases investment in transmission as the density-dependence increases. We have repeated this analysis for a wide range of parameter values and find the pattern to be qualitatively robust to parameter and trade-off choice (though accelerating costs through virulence of increased transmission are required). Furthermore, it can be shown that the decrease in CSS virulence follows directly from the increase in susceptible density, see electronic supplementary material, appendix A and figure 1.

Investment in transmission as functions of the other demographic and epidemiological parameters can be understood from considering the fitness gradient, equation (2.6), which locates the singularity position at the trait satisfying *α*′(*β**) = *X*(*β**)/(1 + *cH*(*β**)). Higher values of *α*′(*β**) correspond to higher evolved transmission and virulence. The term *X*(*β**)/(1 + *cH*(*β**)) which fixes the singularity location is composed of a numerator and denominator, reflecting two opposing selective forces. According to the numerator parameter values that decrease (increase) the susceptible density select for lower (higher) transmission, but simultaneously, according to the denominator, those same values may lead to a decrease (increase) in total host density which selects for higher (lower) transmission because of its effect on DDV. Figure 2*a*–*e* demonstrates that the balance of these selective forces favours an increase in CSS virulence for an increase in all parameters except the birth rate, *a*.

Where an evolutionary singularity has convergence stability but not evolutionary stability, the population will find itself at a fitness minimum, undergo disruptive selection and branch into two distinct strains. In this model, however, despite the additional density dependence, the conditions for this cannot be met and branching cannot occur (see electronic supplementary material, appendix B for details and electronic supplementary material, appendix C for branching in a related model).

### (c) Model incorporating seasonality

To induce fluctuations in the population dynamics, we introduce a seasonally forced reproduction rate of intensity *δ* so that *a* = *a*(*t*).

We present results based on a sinusoidal forcing function *a*(*t*) = *a*_{0}(1+*δ* sin(2*π**t*/*ε*)), following previous studies [31,32], where is the amplitude and *ɛ* is the period of the seasonal forcing which has an average value of *a*_{0}. *δ* = 1 is the maximum amplitude that can be attained without a negative reproduction rate while *δ* = 0 corresponds to a constant reproduction rate (no seasonality).

To emphasize the importance of DDV, let us first consider the standard SIRS model without DDV with *a* = *a*(*t*). This corresponds to our model when *c* = 0. The population dynamics will converge to either a periodic attractor or a chaotic attractor [32]. When the dynamics are periodic the invasion fitness of the system over a period, of length *T*, from *P*_{0} to *P*_{1} [33], where *P*_{1} = *P*_{0} + *T*, is
2.10with denoting the mutant parameter value. We show in electronic supplementary material, appendix D that the singularity position is unaffected by seasonality for the classical SIRS model with periodic dynamics (it can be similarly shown for chaotic dynamics) and this finding has been verified by simulation of the evolutionary process (see White *et al.* [34] for an explanation of the simulation process). The proof in electronic supplementary material, appendix D rests on the fact that the average density of susceptibles over a period is constrained to be constant with respect to the magnitude of seasonality, *δ*, and the period of the forcing, *ε*,
2.11The CSS remains invariant under seasonality when other parameters such as transmission are forced and this applies also to other models beyond the standard SIRS such as the SEIRS (see electronic supplementary material, appendix E). Therefore, a key result is that in infectious disease systems that do not include DDV, the amplitude and period of seasonality will not alter the evolutionary endpoint of virulence.

### (d) Population dynamics of the SIRS model with density-dependent virulence under seasonality

Under DDV (*c* > 0), average susceptible density is no longer conserved. Now,
2.12and the average susceptible density depends on the average total host density which itself may depend on the magnitude of seasonality, *δ*. Owing to complexity of the problem, we cannot analytically find the average population densities needed to show the effect of seasonality on a CSS. Instead, the average densities are calculated numerically. We use the trade-off of §2*b* (see caption of figure 1) between baseline virulence and transmission so that the singularity, when *δ* = 0, is a CSS. Numerically calculated total population densities over a period show that both average susceptible density and average total host density decrease as *δ* increases (figure 3*a*). We have verified that this pattern is consistent across parameter space for periods within reasonable bounds (we assume that this is less than 100 times, and more than one hundredth of, the infectious period) and for parameters that maintain infection within the host population for .

### (e) Parasite evolution under seasonality

For our seasonal model, the invasion fitness is,
2.13and applying the singularity condition implies,
2.14
2.15where substitution of equation (2.12) into equation (2.14) results in equation (2.15). Thus, the value of the CSS, *α**, is fixed by the average total host density over the period. Referring back to the population dynamics shown in figure 3*a*, since average total population density decreases with increasing amplitude of seasonality (and does so across the parameter space) the right-hand side of equation (2.15) is increasing with *δ*. We therefore expect that the CSS value will increase with the amplitude of seasonality. To verify this, we numerically locate the position of the CSS using equation (2.14) for various values of *δ* (figure 3*b*) which confirms that selection favours an increased investment in transmission as the amplitude of seasonality increases. The CSS values are also confirmed using simulations of the adaptive dynamics process (see White *et al.* [34] for an explanation of the simulation process). Furthermore, the increase in virulence that evolves as a result of seasonality in the SIRS system with DDV, depends on the nature of the infection with chronic infection or long-lived hosts (low *γ*, low *b*) as opposed to acute infection or short-lived hosts (high *γ*, high *b*) having the greatest percentage increase in virulence under seasonality (figure 3*c*,*d*). Seasonality therefore has a larger impact on virulence evolution (in the SIRS model with DDV) as the infective period increases.

## 3. Discussion

We have developed a series of models that examine the impact of DDV on the evolution of parasites in both seasonal and non-seasonal environments. Our results may be important in a wide range of natural disease systems because an interaction between host stress and cost of infection, captured here by DDV, is likely to be common [9,11]. We show firstly that parasites are selected to lower exploitation, and therefore lower transmission and virulence as the strength of DDV increases. Our second key result is that seasonal forcing of systems incorporating DDV leads to an increase in exploitation, and therefore higher transmission and virulence, with the impact of seasonality more pronounced for longer infection periods. As such our models emphasize both the importance of context dependence and seasonality in the evolution of parasites.

Under the assumptions of the trade-off hypothesis high transmission follows from intense host exploitation, the fitness gain of which must be counted against a short infectious period associated with severe parasite-induced mortality [4,27]. Theory has shown that intermediate parasite virulence is often favoured when the ensuing trade-off is saturating [3,35,36]. The virulence that evolves is that which optimizes parasite fitness and thus under the trade-off hypothesis it can depend directly on host exploitation strategy, demographic and epidemiological rates. Here, with the inclusion of DDV, optimal fitness also depends indirectly on these rates through their effect on total host density. We find that decreased virulence is selected for as the strength of DDV increases. This is in contrast with increasing natural mortality (and recovery) that leads to the evolution of higher host exploitation (here and in the models of May & Anderson [37], Kakehashi & Yoshinaga [38], Williams & Day [5] and Choo *et al.* [39]). This is an important insight from our models; DDV has the opposite effect on the evolution of parasites to that of higher baseline mortality. It emphasizes the importance of ecological feedbacks in understanding the evolution of virulence. Other theoretical studies on the evolution of virulence [40,41] have found that frequency-dependent transmission leads to qualitatively similar results to the density-dependent case. Here, however, qualitatively different results were found to be possible in our model when the transmission is frequency dependent (results not included). In particular, CSS virulence was found to increase when the trade-off between virulence and transmission was weakly convex.

When seasonality in host birth is included in classical epidemiological models with DDV, we show that the result is an increase in evolved virulence. This outcome is rooted in the effect of the DDV interaction on the average population densities over the seasonal period. In a range of classic baseline models (e.g. SIRS, SEIRS and their variants) lacking this extra density dependence, the average susceptible population density in the fitness expression over a period of any type of forcing (and for a multiplicity of forced terms) remains constant with respect to the magnitude of seasonality. Hence, deterministically, the evolution of virulence is unchanged by seasonality except when density dependence acting on virulence or recovery breaks the invariance of the average susceptible population density. This is consistent with the results of Williams & Dye [42] where *R*_{0} for the classical SIR model with seasonal transmission was found to be independent of the amplitude of seasonality.

The effect of seasonality on the evolution of virulence was previously associated with the intuition that periods of low density ought to purge those parasites with more virulent strategies. This is a stochastic argument based on the most virulent parasites showing the greatest vulnerability to stochastic extinction during the periods of low population density. The deterministic effect, hitherto ignored, when DDV is present, is a marked increase in evolved virulence. Jokela *et al.* [7] discuss the intuition when they examine how anoxic conditions and starvation of hosts affect virulence of two closely related trematode parasites, *Rhipidocotyle campanula* and *Rhipidocotyle fennica*, of the freshwater clam, *Anodonta piscanalis* (sampled from temperate lakes in Finland). Seasonal starvation and anoxia in these Finnish lakes were implicated in stress-dependent virulence which was found to be prevalent in all populations. Contrary to the hypothesis that seasonality should lead to lower exploitation strategies through a purging of more virulent strains, high virulence was observed in multiple populations of *Rhipidocotyle campanula*. The authors suggest that this may be due to eradication of the infection from the clam population in periods of high stress and repopulation from refugia of infection in a separate host species. Here, our theoretical results offer an alternative explanation, since acute virulence is predicted by our models under the similar conditions of seasonal DDV.

In a related model assuming evolving baseline virulence but this time with non-evolving DDV, evolutionary branching has the potential to occur (see electronic supplementary material, appendix C). The requisite condition for branching of two strategies as proposed in Heino *et al.* [43] (see also [44]), namely a feedback dimension (to *per capita* growth of infecteds) of two, is present in this model. This is reminiscent of branching found in an SI type model when density-dependent natural mortality is included (Svennungsen & Kisdi [45] which is an extension of Pugliese [46]; see also [47]). The result extends the range of disease models for which branching has been demonstrated; among which already are counted models incorporating density-dependent death [45,46], superinfection [48,49], ‘specialist’ parasitism [50–52] and selective predation [53]. When DDV evolves together with baseline virulence in our main model, however, and both population feedbacks have adaptive traits, mutual invasibility becomes impossible and no branching can occur.

To conclude, our models emphasize the importance of context dependency and seasonality in the evolution of infectious disease. In particular, we highlight the need for both experimentalists and theoreticians in evolutionary epidemiology to allow for the possibility of a density-dependent virulence, as theory predicts its presence may significantly alter the evolutionary outcome of such systems.

- Received October 16, 2012.
- Accepted October 31, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.