Evolution of class-structured populations in periodic environments.

What is the influence of periodic environmental fluctuations on life-history evolution? We present a general theoretical framework to understand and predict the long-term evolution of life-history traits under a broad range of ecological scenarios. Specifically, we investigate how periodic fluctuations affect selection when the population is also structured in distinct classes. This analysis yields time-varying selection gradients that clarify the influence of the fluctuations of the environment on the competitive ability of a specific life-history mutation. We use this framework to analyse the evolution of key life-history traits of pathogens. We examine three different epidemiological scenarios and we show how periodic fluctuations of the environment can affect the evolution of virulence and transmission as well as the preference for different hosts. These examples yield new and testable predictions on pathogen evolution, and illustrate how our approach can provide a better understanding of the evolutionary consequences of time-varying environmental fluctuations in a broad range of scenarios.

Many organisms experience periodic fluctuations of their environment. These fluctuations may be driven by abiotic variations of the environment at different time scales (e.g., diurnal and seasonal variability), or by the dynamics of biotic interactions between organisms (e.g., predator–prey or host–parasite limit cycles). In such periodically changing environments, selection on life‐history traits is likely to fluctuate over time, but we currently lack a good understanding of the feedback between periodic environmental dynamics and long‐term phenotypic evolution (Barraquand et al., 2017).

A good measure of selection in periodic environments should tell us whether, on average over one period of the fluctuation, a mutation increases or decreases in frequency. But how should we compute this average fitness when selection may vary both in time but also among different classes of individuals? Floquet theory provides an answer to this question through the computation of the invasion fitness of a rare mutant in the periodic environment produced by the wild type (Klausmeier, 2008; Metz et al., 1992; Meszéna et al., 2005; Metz, 2008). However, the analysis based on Floquet theory is numerical and yields little biological insight. It only provides a good understanding of evolutionary dynamics when a single class of individuals is needed to describe the mutant dynamics. In this case the invasion fitness is simply the average, over one period, of the per‐capita growth rate of a rare mutant (Cornet et al., 2014; Donnelly et al., 2013; Ferris & Best, 2018; Gandon, 2016; Kremer & Klausmeier, 2013; Metz, 2008; Pigeault et al., 2018). In class‐structured populations, however, the lack of an analytical expression for the invasion fitness hampers the biological interpretation of the results obtained with Floquet analysis. In this paper, we fill this gap and provide a new method to analyse selection in class‐structured populations subject to periodic environmental fluctuations.

In constant environments, it has been shown that the direction of selection should depend on the relative abundance of each class as well as the productivity of the focal organism in each class, so that we need to keep track of both the quantity and the quality of different classes (Gandon, 2004; Lehmann & Rousset, 2014; Lion, 2018a; Rousset, 2004; Taylor, 1990; Taylor & Frank, 1996). Our approach extends this idea to periodic environments and allows us to derive, using only the standard weak‐selection assumption, an expression of the selection gradient in terms of the quantity and quality of classes, which are now time‐dependent variables (see Fig. 1 for a graphical summary). With this approach, the selection gradients in periodic and constant environments are directly comparable and conceptually similar.

Figure 1

Graphical summary of our approach. Environmental fluctuations are used to evaluate the quality and quantity of individuals in the different classes at a given time. This information can then be used to calculate the selection gradient.

We first provide a general description of eco‐evolutionary dynamics in class‐structured, polymorphic populations, then turn to the dynamics of a mutant invading a resident population. We show that, under weak selection (that is, for mutations of small phenotypic effects), a separation of time scales argument can be used to derive the selection gradient in periodically varying environments. To illustrate the potential use of this approach, we focus on the evolution of pathogen life‐history traits (such as transmission and virulence) in three different epidemiological scenarios when there is periodic variation in the availability of susceptible hosts. The focus on pathogens is not restrictive, and the method can be applied to a variety of life cycles. Evolutionary epidemiology, however, provides a very natural framework in which to think about potentially complex eco‐evolutionary feedbacks.

Eco‐Evolutionary Dynamics

We consider a focal population composed of K different “classes” of individuals. For instance, these different classes may correspond to distinct developmental stages of the organism (e.g., young and old, male or female), different immune states, or different locations in a spatially structured environment. Because we are interested in evolution we also assume that the population is composed of M different “genotypes.”

The life cycle is defined by a matrix of average transition rates, , which refer to the net production of class‐k individuals by class‐j individuals, averaged over all genotypes (Lion, 2018b, 2018a). These transitions can be due to reproduction, mortality, maturation, or dispersal depending on the biological context. Crucially, these rates can vary with a change in the environment which is referred to as (Lion, 2018b; Metz et al., 1992; Metz, 2008). These environmental variations may be driven by density‐dependent effects caused by changes in population densities , by frequency‐dependent effects caused by changes in the frequencies of the different types, but also by changes in extrinsic variables (such as the density of a resource) which we refer to as , so that . Table 1 gives a summary of the main notations.

Table 1

Definition of the main mathematical symbols

Mathematical symbol	Description
n^k(t)	Density of individuals in class k (1≤k≤K)
n(t)=∑_kn^k(t)	Total density of individuals
f^k(t)=n^k(t)/n(t)	Frequency of individuals in class k (with respect to the total population)
v^k(t)	Individual reproductive value for class k
n(t)	Vector of class densities n^k(t)
f(t)	Vector of class frequencies f^k(t)
v(t)	Vector of individual reproductive values v^k(t)
e(t)	Vector of external ecological variables
fik(t)	Frequency of genotype i within class k (1≤i≤M)
p(t)	Vector of genotype frequencies fik (dimension K×M)
f_i(t)=∑_kfik(t)f^k(t)	Global frequency of genotype i
f∼_i(t)=∑_kv^k(t)fik(t)f^k(t)	Reproductive‐value‐weighted frequency of genotype i
rikj	Per‐capita transition rate of genotype‐i individuals from class j to class k
r¯^kj=∑_irikjfij	Average per‐capita transition rate at from class j to class k
R¯	Matrix of average per‐capita transition rates r¯^kj.
R_w	Matrix of per‐capita transition rates for the resident type rwkj.

The average transition rates depend on the transition rates of the various genotypes, , and on the frequencies of genotype i in class j. Thus, we have Note that the vital rates may themselves be time‐dependent, hence we make time an explicit argument of the transition rates . Following Lion (2018b, 2018a), this yields the following eco‐evolutionary dynamics:

where is the fraction of individuals in class j at time t, and is the matrix of average per‐capita transition rates between different classes of individuals. Note that if the environment depends on extrinsic variables (e.g., the density of a resource or a predator), one needs to specify the dynamics of to complete the characterisation of the dynamical system (2). Thus, the eco‐evolutionary dynamics are described by the equations of system (2), plus the equations needed to describe the dynamics of extrinsic variables.

DYNAMICS OF MUTANT FREQUENCIES

Now suppose that, for simplicity, we only have two types in the population: a resident wild type (w) and a mutant (m). The change in the global frequency of the mutant, , can then be decomposed as follows (Day & Gandon, 2006; Lion & Gandon, 2016; Lion, 2018a; Osnas et al., 2015). The first line represents the average effect of within‐class selection, given by the genetic variance within class j, times the difference in transition rates from class j to all other classes. The average is taken over the class distribution . In contrast, the second line represents the effect of gene flow between classes, which depends on the relative contribution of the different classes when there is variation in genotype frequencies (differentiation) among classes. Importantly, this second term conflates both the effect of selection (which can shape the differentiation ) and of purely demographic processes (what Grafen (2015) termed “passive changes” in allele frequencies), so that any estimation of selection based on equation (3) may be biased by the existence of intrinsic differences in qualities between classes (Gardner, 2015; Grafen, 2015; Lion, 2018a).

Thus, although there can be value in explicitly tracking the dynamics of the genetic differentiation between classes (see e.g., Berngruber et al. 2013, 2015; Lion & Gandon 2016), it is for our purpose more convenient to use an alternative measure of the mutant frequency, which weights the mutant frequency in class k by the quantity () and quality () of individuals in class k at time t. Specifically, we use the reproductive value of an individual in class k at time t as the measure of quality . As previously shown (Lion, 2018a), the dynamics of this weighted average frequency can then be written as In contrast to equation (3), equation (5) describes the net effect of selection on the change in mutant frequency, without the confounding variations due to “passive changes” (Gardner, 2015; Grafen, 2015; Lion, 2018a). It thus provides a more convenient measure of selection where the overall change in mutant frequency is driven by the sum of the effects of the mutation on the rates weighted by the frequency of class j and the individual reproductive value of class k. Crucially, both the quantity and the quality of the different classes can change in a fluctuating environment, so, unlike classical equilibrium theory we need to characterise these fluctuations to understand and predict life‐history evolution. In order to characterise these fluctutations, we will use a weak selection assumption to decouple the ecological fluctuations from the evolutionary dynamics.

Weak Selection Approximation

In this section, we use a weak‐selection approximation of equation (5) to derive an expression of the selection gradient using dynamical reproductive values. Our approach is based on a separation of time scales, which readily occurs when the mutant phenotype, , is close to the resident phenotype , so that , where ε is small. Under this separation of time scales, the population will settle on a population dynamical attractor, such as a fixed point or a limit cycle, and we can compute the quality and quantity of each class at any given time on the attractor. The fixed point case corresponds to the classical theory developed for constant environments (Otto & Day, 2007; Rousset, 2004; Taylor, 1990), and the limit cycle case corresponds to an extension of this theory to periodic population dynamics, which we now present.

SELECTION GRADIENT

In Appendix A, we show that the class densities , class frequencies , and individual reproductive values are “fast” variables. This means that, under weak selection, the environment settles on a periodic attractor which is well approximated by the attractor of a monomorphic resident population, . On this monomorphic attractor, we only need the densities and the external variables to fully characterise the environment experienced by the mutant type.

In contrast, the weighted mutant frequency is a slow variable. To first‐order in ε, we can write the dynamics of as where the instantaneous selection gradient is given by Note that, in Equation (6), we have replaced the class variances in Equation (5) by the population variance . This is because the differences will typically be and therefore this substitution will only contribute an extra term to equation (6).

As a result of periodic fluctuations on the fast time scale, the value and sign of the selection gradient may fluctuate. However, under weak selection, these fast fluctuations can be averaged out. This is known as the averaging principle (see, e.g., Cai & Geritz 2020), which allows us to approximate the dynamics of the mutant frequency on the slow time scale by the solution of the so‐called averaged system where is the selection gradient averaged over one period of the resident attractor (i.e. for any given τ on the periodic attractor).

The selection gradient (9) takes the form of a sum of marginal selective effects (giving the influence of the evolving trait on between‐class transitions) weighted by the time‐varying quantity and quality of different classes and is therefore reminiscent of the expressions obtained in equilibrium class‐structured populations. In fact, for constant environments, the class frequencies, reproductive values and environment are constant, and therefore equation (9) exactly reduces to the classical expression of the selection gradient for class‐structured equilibrium populations (Otto & Day, 2007; Rousset, 2004; Taylor, 1990). For periodic environments, equation (9) provides an analytical approximation, for weak selection, of the invasion fitness of a mutant typically calculated as a Floquet exponent (Appendix B). Note that this first‐order approximation gives information on the direction of selection and its potential evolutionary endpoints, but not on the evolutionarily stability of these singularities, for which a numerical computation of the Floquet exponent is still needed (Appendix B).

In equation (9), the class frequencies and reproductive values are calculated on the attractor of the monomorphic resident population. In the next two sections, we show how these quantities can be calculated using dynamical equations of the resident population.

DYNAMICS OF CLASS FREQUENCIES

Following Lion (2018a), the class frequencies follow the dynamics: where is the per‐capita growth rate of the monomorphic resident population and is the matrix of resident transition rates . If we focus on the frequency of class j we have: Hence, the change in frequency of class j depends on the current class frequencies and on the relative contributions of all the classes to class j. Note that we need to account for the overall growth rate of the whole population. This is because we are monitoring the change in class frequencies, not class densities. The sum of class densities can increase (i.e. when ) but the sum of class frequencies must remain equal to 1 at all times (i.e. ).

Mathematically, we can calculate the class frequencies at time t by integrating equation (10) forward in time from an appropriate initial condition (i.e. an initial condition that leads to a biologically relevant periodic attractor after some time). For instance, we can use an initial condition where all the are equal to indicate that there is no a priori information on the relative quantities of the different classes. In other words, the quantity of class j at time t depends on the past trajectory of the population (Figure 2, top panels).

Figure 2

Calculating the class frequencies and reproductive values on the resident periodic attractor. Top panels: The class frequencies are calculated forward in time using equation (11) (left). Integrating from to , this gives dynamics that settle on a periodic attractor (right). In this figure, we assume that both classes are initially equifrequent (i.e. ) which amounts to saying that we have no a priori information about their relative abundance. – Bottom panels: Next, we insert the values into equation (13) to calculate the individual reproductive values. This is a backward process, because we need the values at t to compute the values at (left). We therefore integrate equation (12) using the final condition , which means that there is initially (i.e. at time ) no a priori information on the relative quality of the two classes. We can then choose a period of the attractor (e.g. the gray zone) and calculate the values and for any t in the period.

DYNAMICS OF INDIVIDUAL REPRODUCTIVE VALUES

Similarly, the individual reproductive values, collected in the row vector , follow the dynamics (Lion, 2018a): If we focus on the individual reproductive value in class j, we have: Mathematically, equation (12) is the adjoint of equation (10), which generalises the fact that, in constant environments, class frequencies and individual reproductive values are respectively right and left eigenvectors of the transition matrix of the resident population. Biologically, gives the relative quality of class k at time t. As in equation (10), we take into account the change in total population size (through ) so that and are co‐normalised at all times (i.e. . In other words, the reproductive value of an average individual, is equal to 1 at all times, and measures the contribution to the future of the population of an individual sampled in class k at time t, relative to a randomly sampled individual.

Hence, measuring the “quality” of a class at time t depends on the descendants and thus on the “future” of the population. In contrast with class frequencies, we compute the individual reproductive values by integrating equation (12) backward in time, from a terminal condition at a time in the distant future (Figure 2, bottom panels). Typically, we set for all classes at time (Barton & Etheridge, 2011; Lion, 2018a). This terminal condition indicates that there is no a priori information on the relative qualities of the different classes of individuals at that point in time and we use between‐class transitions on the periodic attractor to acquire information on these relative qualities, by integrating the system (12) backward in time.

A GENERAL RECIPE TO STUDY EVOLUTION IN PERIODIC ENVIRONMENTS

Equation (9) is the central result of our article and leads to a general recipe to study the evolution of life‐history traits in class‐structured populations experiencing periodic environmental fluctuations. The main steps of this method are summarised in Box 1. We start by computing the dynamics of the class densities in the resident monomorphic population, which allows us to calculate the class frequencies at a given time on the periodic attractor (i.e. in the gray zone in Fig. 2). We then use this information to calculate the dynamics of the individual reproductive values, which depend on the class frequencies through the computation of the transition rates . Although in general these solutions can only be computed numerically, we show in the next section that algebraic manipulations of equations (9), (11) and (13) can be used to shed light on the interplay between fluctuating environmental dynamics and selection. We do so by applying our general method to the evolution of pathogen life‐history traits when the density of hosts fluctuates.

Box 1: How to study life‐history evolution in periodically fluctuating environments?

We detail below the different steps allowing us to derive the selection gradient driving life‐history evolution in a class‐structured population in a periodic environment.

Step 1 ‐ Formalise the description of the life cycle of the focal organism in a monomorphic resident population. This is captured in the dynamical system: where is the matrix of between‐class transition rates in the resident population, together with the dynamical equations for the extrinsic variables in .

Step 2 ‐ Determine which life‐history traits are under selection (and the trade‐off between these traits). Indeed, multiple life‐history traits may be involved in the life‐cycle and it is crucial to be very explicit about the constraints acting on the evolving traits.

Step 3 ‐ Derive the selection gradient using equation (9). This step may often yield insight through the decomposition of selection into biologically meaningful quantities (see section 4 for examples). But in general, moving to steps 4 and 5 is needed to predict the influence of periodic fluctuations on life‐history evolution.

Step 4 ‐ Use equation (A) and (11) to determine the periodic attractor of the resident population. This will yield the forward dynamics of class frequencies.

Step 5 ‐ Use equation (12) and the numerical solution derived in step 4 to solve the backward dynamics of individual reproductive values (using a final condition where individual reproductive values are all equal to 1).

Step 6 ‐ The results of steps 4 and 5 can then be plugged into the selection gradient (step 3) to obtain a mathematical expression (or a numerical computation) of the selection gradient on the trait. This can be used to identify evolutionary singularities, the evolutionarily stability of which can be checked using Floquet analysis (Appendix B).

Pathogen Evolution in Periodic Environments

Many pathogens have to cope with environmental fluctuations (Altizer et al., 2006; Martinez, 2018). In this section, we use the approach in Box 1 to study pathogen evolution under three distinct epidemiological scenarios corresponding to three different pathogen life cycles. In all scenarios the pathogen can be present in two distinct classes of hosts (A and B) and the environmental fluctuations are captured by a periodic function that gives the probability of production of susceptible hosts at time t. Hence, using our general terminology, the forcing function causes periodic fluctuations in the vector (collecting the densities of susceptible hosts), which in turn drives fluctuations in the vector (collecting the densities of infected hosts in classes A and B). The general approach outlined in Box 1 can be applied to any periodic function, but for simplicity we consider a smooth version of a step function with minimum 0, maximum , and mean (Figure 3).

Figure 3

Environmental fluctuation. We model fluctuations in the environment with a smoothed step function with minimum 0, maximum and mean . This is captured by the following function: . We typically use and . Note that corresponds to a constant environment with value .

SCENARIO 1: THE CURSE OF THE PHARAOH HYPOTHESIS

The claim that long‐lived pathogen propagules could select for higher pathogen virulence has often been presented as the “Curse of the Pharaoh hypothesis” (Bonhoeffer et al., 1996; Gandon, 1998). Previous theoretical analyses focused mainly on temporally constant environments. Here we want to analyse the influence of fluctuations on the availability of susceptible hosts on the evolution of virulence for pathogens with free‐living stages.

Let us assume that there is a single class of susceptible hosts S but two different classes of the pathogens, which are (i) the infected host () and (ii) the propagule stage (which lives outside the infected host but that we still denote for consistency with the general framework). Thus, the vector of population densities is and is the extrinsic environmental variable . The epidemiological dynamics are given by where refers to the periodic influx of susceptible hosts in the population. A susceptible host becomes infected upon contact with the propagules, which occurs at rate . Infected hosts die at rate , where is the background mortality rate and represents virulence, and produce propagules at rate . These propagules die at a rate . Next, we follow the steps presented in Box 1 to explore how higher propagule survival (e.g. lower values of ) can affect the evolution of pathogen virulence in the infected host.

Step 1

We can use equations (14) to derive the matrix which captures the transition rates between classes A and B in a monomorphic resident population:

Step 2

At this stage it is important to specify the constraints acting on the underlying trait z that controls the evolution of pathogen virulence . As in previous studies, we assume that an increase in virulence is associated with an increase in pathogen transmission rate and in particular on the production of propagules (so that the derivatives of and with respect to z are both positive; Bonhoeffer et al. 1996; Day & Gandon 2006).

Step 3

We use equation (7) to obtain the instantaneous selection gradient on the trait z: The first term between brackets represents the gain in fitness if a pathogen in class A invests in the production of propagules (weighted by the individual reproductive value of propagules). The second term between brackets accounts for the loss in fitness if a pathogen in class A dies (weighted by the individual reproductive value ). The selection gradient on z is obtained by integrating equation (16) over one period of the fluctuation which yields: An analysis of the dynamics of class frequencies and individual reproductive values is required to better understand selection on virulence and transmission.

Steps 4 and 5

For our life cycle, equation (12), which gives the dynamics of reproductive values, can be written as: Fluctuations in the availability of susceptible hosts cause fluctuations of the reproductive values which can be obtained numerically from (18), using the solution of system (14) to evaluate and . But equation (18) also yields a useful expression for the average of the ratio of individual reproductive values (Appendix S1): The right‐hand side of equation (19) is the ratio of individual reproductive values in the absence of fluctuations, so that we see that the ratio fluctuates around its equilibrium value in a constant environment.

Step 6: Equation (19) can be used to rewrite the selection gradient (17) which yields: where is the class reproductive value at time t and refers to the gradient of selection in the absence of fluctuations. This expression shows that, as pointed out by earlier studies (Bonhoeffer et al., 1996; Day & Gandon, 2006), the mortality rate of the propagule has no effect on the evolutionary stable virulence in a constant environment. However, with periodic fluctuations, we will see that the mortality rate of the propagule does affect the evolutionary outcome. Equation (20) is particularly insightful because it shows how periodic fluctuations of the environment can affect the evolution of the pathogen: because , periodic fluctuations will affect the direction of selection only if the temporal covariance between and is non‐zero. In words, this means that, if infected hosts tend to be more abundant at times when propagules are relatively more valuable (e.g. the covariance is positive), the ES virulence will be higher than in a constant environment because the pathogen then reaps fitness benefits from increased propagule production.

Numerical simulations show that this covariance is expected to be positive when is low and negative when is large (Figures 4a, 4b, 4d; see the SOM for an attempt to understand the sign of this covariance). In contrast with the analysis of Bonhoeffer et al. (1996), we thus expect fluctuations to alter the predictions regarding the influence of on the evolution of virulence and transmission rates. As shown in Figure 4c, higher rates of survival (i.e., lower values of ) tend to select for higher virulence and transmission (as in the Curse of the Pharaoh hypothesis). However, this effect is non‐monotonic: when gets very low (i.e. when propagule live very long), the covariance vanishes because and thus , and therefore the ES virulence is the same as predicted in a constant environment. Note that Figure 4c also shows that, as expected, the ESS predicted from the selection gradient using time‐dependent reproductive values is consistent with the value predicted from a more standard Floquet analysis (dashed line).

Figure 4

Scenario 1: The Curse of Pharaoh. (a) Parametric plot of versus for and . The grey diamond gives the equilibrium for a constant environment. The slope of the regression line (dashed) is proportional to the temporal covariance between and . (b) Same as (a) but with . (c) Predicted ESS as a function of the mortality rate of propagules, using the reproductive‐value‐based approach (dots) and the Floquet exponent (dashed line). The dotted line gives the prediction of the corresponding equilibrium model, . (d) Temporal covariance between the class reproductive value and the ratio of individual reproductive values as a function of for , where is the ESS in the absence of fluctuations. Parameters: , , , , , , .

Conclusion of Scenario 1

In a constant environment the longevity of pathogen propagules does not affect the long‐term evolution of pathogen virulence. In contrast, we show that periodic fluctuations of the environment (i.e., fluctuations in the availability of susceptible hosts) can select for higher (when the mortality rate of propagules is low) or lower (when the mortality rate of propagules is high) pathogen virulence. The selection gradient given in equation (20) allows us to capture the effect of fluctuations through a single temporal covariance which measures the deviation from the selection gradient in a constant environment.

SCENARIO 2: HOST PREFERENCE

Many pathogens can infect several host species. Intuitively, which host the pathogen should prefer will depend on the relative qualities of the hosts. But what would be the best strategy when the qualities or abundances of the different host species fluctuate? In our second scenario, we investigate the effect of periodic fluctuations in the abundance of different host species on the evolution of the pathogen's preference strategy.

We consider an epidemiological model with contact transmission (no free‐living propagules) and assume that the pathogen can exploit two different hosts (A or B). When a pathogen enters a susceptible host (resp. ), we assume that the infection is successful with probability (resp. ). With these assumptions, the dynamical system becomes:

where is the force of infection. The two infection routes differ by their epidemiological parameters, so that one infection route may be more contagious or virulent than the other. Furthermore, we assume that the production of susceptible hosts is periodic, and such that and , where is the probability of production of B hosts (Figure 3).

From equations (21), we derive the transition matrix

We assume that there is a trade‐off between and such that . Hence, if , infection is only possible on A hosts, while if , both host classes are equally susceptible to infection. The trait z can thus be interpreted as measuring preference towards A hosts. For simplicity, we assume that the pathogen's virulence is higher in host A, but its transmissibility is independent of the host (i.e. , but ).

Based on these assumptions on the life cycle, a naive prediction could be that the pathogen should always prefer the “good” host B, in which it enjoys a longer lifespan because the pathogen is less virulent on this host. However, the optimal strategy depends on the relative availability of the two classes of hosts, which can fluctuate over time. To better understand the selective pressures on the preference trait, we use our approach to derive the selection gradient at time t and obtain: The terms have a simple interpretation as the expected reproductive output of a pathogen propagule at time t through class k. Thus, the direction of selection is determined by whether this reproductive output is larger through class A than through class B. Furthermore, potential evolutionary endpoints satisfy the balance condition which simply states that selection halts when average reproductive outputs are the same in both classes.

To better understand the impact of periodic fluctuations in host availability on the evolution of the preference trait, we now need to numerically calculate the dynamics of the reproductive values and densities of susceptible hosts. A detailed discussion of these dynamics is given in the SOM, but we show in Figure 5a,b that the period of has a strong impact on the dynamics of the difference in reproductive output . When the period is small (Figure 5a), D fluctuates rapidly around a mean that is close to its value in a constant environment (which is zero). In contrast, for large periods (Figure 5b), better tracks the environmental fluctuations (it is minimal when there are only B hosts and maximal when there are only A hosts) and its mean is greater than zero.

Figure 5

Scenario 2: host preference. (a) Dynamics of for and . The dashed line indicates the mean value of D, the dotted line corresponds to the equilibrium value of D in the constant environment, which is zero, and the grey line gives the dynamics of . (b) Same as (a) but with . (c) Predicted ESS as a function of the period of using the reproductive‐value‐based approach (dots) and the Floquet exponent (dashed line). The dotted line gives the prediction of the corresponding equilibrium model, . (d) Mean value of the difference in reproductive output as a function of the period for . Parameters: , , , , , .

Step 6

Using the dynamics of the ecological variables, it is possible to calculate the ES strategy for different values of T (Appendix S2) and Figure 5c shows that the ESS increases with the period of fluctuations. For small periods, the ecological variables fluctuate rapidly and the ESS is close to the prediction of the constant model: host preference is biased towards the B hosts because pathogens infecting these hosts have a higher reproductive value. In contrast, when the period of fluctuations is large, host preference is biased towards the A hosts even though their reproductive value is always lower (Figure S2). To understand this counter intuitive result it is important to note that, with the periodic fluctuation we consider, the mean gradient of selection is simply the average over a half‐period dominated by host A and a half‐period dominated by host B. When the period is large, the epidemiological dynamics reach an endemic equilibrium within each half‐period, and this yields: where and are the densities of the susceptible hosts at the endemic equilibrium of the corresponding single‐class models. It is possible to show that and . Solving then yields the evolutionary stable strategy which predicts that, indeed, this preference strategy is biased towards host A which suffers more from the infection (i.e., ).

Conclusion of Scenario 2

In a constant environment the pathogen evolves a preference for the host which suffers less from the infection because it prolongs the duration of infection. In contrast, we show that slow fluctuations in the abundance of the two hosts can select for the opposite strategy where pathogens evolve a preference for the host that suffers more from the infection. The selection gradient given in equation (23) allows us to shed some light on the effect of temporal fluctuations on pathogen evolution.

SCENARIO 3: IMPERFECT VACCINES

The use of imperfect vaccines may affect the evolution of pathogen virulence and transmission. These evolutionary consequences have been studied by Gandon et al. (2001, 2003) when the coverage of vaccination does not fluctuate in time. In our third epidemiological scenario, we ask how periodic fluctuations in vaccination coverage may affect the evolution of virulence, building on a recent study by Walter & Lion (2021).

We consider the same epidemiological dynamics as in Scenario 2 (equation (21)), but assume that hosts are inoculated at birth with an imperfect vaccine at a rate that fluctuates periodically (Figure 3). This yields a fluctuating influx of A hosts that are unvaccinated () and B hosts that are vaccinated ().

Steps 1 and 2

As in scenario 2, the transition matrix is given by equation (22). We now assume that , and consider that the vaccine can either act by reducing the transmissibility of hosts B () or by decreasing virulence (). This corresponds to the anti‐transmission and anti‐virulence vaccines introduced by Gandon et al. (2001) (noted r 3 and r 4 respectively in that paper). Finally, we assume that, as in Scenario 1, the trait under selection is the pathogen strategy of host exploitation, z, and that transmission and virulence both depend on z.

With these assumptions, the selection gradient at time t takes the form where and are the class reproductive values in class A and B respectively, and Note that has a useful intepretation: the denominator gives the quality of an “adult” pathogen in class k, while the numerator is the expected quality of a pathogen propagule and therefore quantifies the reproductive value of an “offspring” pathogen. So gives a measure of how valuable reproduction is compared to survival in class k at any given time.

Steps 4, 5 and 6: In general, the densities and reproductive values have complex periodic dynamics (Figures S3 and S4). However, using the dynamics of reproductive values, it is possible to analytically show (Appendix S3) that, in the resident population on its periodic attractor: The term on the right‐hand‐side is one over the basic reproduction ratio of a pathogen when only class k is present, and corresponds to the equilibrium value in a model with constant vaccination coverage.

This useful result allows us to rewrite the average selection gradient as where we have neglected the covariances that arise when taking the mean. Extensive numerical simulations show that this approximation fits very well the prediction of a Floquet analysis (Walter & Lion, 2021).

A full analysis is beyond the scope of this paper, and we refer the reader to the more complete study by Walter & Lion (2021). Nonetheless, is is interesting to see that the selection gradient takes the form of a weighted sum of the selection gradients in class A and B respectively, exactly as in constant environments. For a constant vaccination coverage, the ES strategy of host exploitation is a weighted mean of the optima in class A and in class B, respectively given by the zeros of and . For our anti‐transmission vaccine () , and therefore vaccination has no effect on the optimum. Equation (29) shows that the same holds true for periodic environment, and this is confirmed by numerical calculations of the Floquet exponents (Walter & Lion, 2021). For a vaccine that reduces virulence (), however, , and the position of the ESS is determined by a single variable, which is the mean value, over one period, of the class reproductive value in the resident population.

Figure 6b shows that the class reproductive value closely tracks the environmental fluctuation when T is large (it is close to 1 when , so that only A hosts are produced, and close to zero otherwise), whereas for short periods it quickly fluctuates around its value in a constant environment (Figure 6a). An interesting consequence is that, when the period of fluctuations increases, increases (Figure 6d), which selects for lower virulence compared to a scenario with constant vaccination coverage (Figure 6c). For large periods, converges towards 1/2 (the mean of ), which allows the ES virulence to be analytically calculated (Walter & Lion (2021), appendix S3). Note that, in the latter figure, the slight quantitative discrepancy between the prediction of equation (29) and the Floquet analysis is due to the fact that we have neglected the covariances .

Figure 6

Scenario 3: imperfect vaccines and virulence. (a) Dynamics of the class reproductive value for and . The dashed line indicates the mean value, the dotted line the value in a constant environment, and the grey line gives the dynamics of . (b) Same as (a) but with . (c) Predicted ESS as a function of the period of using the reproductive‐value‐based approach (dots) and the Floquet exponent (dashed line). The dotted line gives the prediction of the corresponding equilibrium model, . (d) Mean value of the class reproductive value as a function of the period for . Parameters: , , , , , , .

Conclusion of Scenario 3

In a constant environment, vaccination with an imperfect vaccine that affects the within‐host growth of the pathogen can select for higher virulence. Here we show that temporal fluctuations in the proportion of vaccinated hosts can mitigate this effect. Importantly, the deviation from the prediction in the absence of fluctuations can be captured by a single quantity, which is the average, over one period, of the class reproductive value (see the equation of the selection gradient (29)).

Discussion

We present a theoretical framework to analyse evolution in periodically fluctuating environments in class‐structured populations, and use it to study the evolution of pathogen traits in three epidemiological scenarios. In Scenario 1, we revisit the “Curse of Pharaoh” hypothesis and show that, while propagule longevity is not predicted to affect the evolution of virulence when the environment is constant (Bonhoeffer et al., 1996; Day & Gandon, 2006), fluctuations in the density of susceptible hosts can strongly alter the predictions on the effect of propagule longevity on the evolution of virulence. In Scenario 2, we show that periodic fluctuations in the availability of two types of hosts can bias pathogen preference towards the host where the pathogen has higher virulence, in contrast to the prediction in a temporally constant environment. Finally, in Scenario 3, we show that the evolution of virulence in response to the use of imperfect vaccines is affected by periodic fluctuations in vaccination coverage, which can select for lower virulence compared to constant environments. These three scenarios illustrate how the derivation of the selection gradient (using the general recipe detailed in Box 1) can yield important insights on the influence of temporal fluctuations on life‐history evolution.

From a methodological standpoint, our analysis extends previous adaptive dynamics studies, which use the Lyapunov (or Floquet, for periodic environments) exponent of a rare mutant as a measure of invasion fitness (Cornet et al., 2014; Donnelly et al., 2013; Ferris & Best, 2018; Geritz et al., 1998; Gandon, 2016; Klausmeier, 2008; Metz et al., 1992; Metz, 2008; Pigeault et al., 2018). Here, we circumvent this problem and derive an expression of the selection gradient on life‐history traits, which represents the first‐order approximation of invasion fitness for weak selection. We show how it can be used to better understand the impact of environmental fluctuations on the direction of selection and the potential evolutionary endpoints of life‐history evolution. In contrast with previous optimisation approaches (e.g. McNamara (1997)), our method does not rely on the assumption that there is no frequency‐ or density‐dependence in the population, and allows us to take into account a variety of ecological feedbacks.

Our approach is based on the analysis of the dynamics of a reproductive‐value‐weighted frequency of a mutant. This is a recurring idea in evolutionary biology (Fisher, 1930; Gardner, 2015; Grafen, 2015; Lion, 2018a; Lehmann & Rousset, 2014; Taylor & Frank, 1996), but, in contrast to previous approaches (see e.g. Gardner (2014)), the novelty here is that we use a dynamical definition of reproductive value to quantify the fluctuating quality of a class in a periodic environment (see Lion (2018a) for a general discussion on this topic, and Brommer et al. (2000); Bacaër & Abdurahman (2008); Caswell (2001) for other approaches). The resulting expression of the selection gradient can then be obtained by weighting the effect, at time t, of a mutation on the transition rates from class j to k by the frequency of class j at time t and by the individual reproductive value of class k at time t. This has two main implications. First, this allows us to quantify selection at time t in terms of the quantity and quality of each class. As in constant environments, the quantity of a class is given by the frequency of individuals in that class, and the quality of a class is measured by their reproductive value, which gives the relative share of future descendants left by individuals in that class. The only difference is that the quantity and quality of each class are allowed to fluctuate periodically over time. Second, the direction of selection can be obtained by computing the average, over one period of the resident attractor, of the instantaneous selection gradient, and we recover a periodic extension of the “invasion implies fixation” principle (Cai & Geritz, 2020; Geritz, 2005; Priklopil & Lehmann, 2020).

It could be argued that, with our approach, the problem of numerically computing a Floquet exponent is replaced by the problem of numerically computing the time‐dependent reproductive values and class frequencies. This is of course true if we are simply interested in the quantitative result, but the expression of the selection gradient allows for a qualitative discussion of the selective pressures. The different examples we examine above illustrate that this approach allows us to better understand the often counter‐intuitive effects of the periodic fluctuations of the environment on life‐history evolution. However, as our approach is currently limited to first‐order effects (e.g. convergence stability), the numerical calculation of invasion fitness using Floquet exponents is required to evaluate evolutionary stability and determine whether the predicted evolutionary singularities are ESSs or branching points (Appendix B).

Although periodic environmental fluctuations are important in nature, many organisms also experience stochastic environmental fluctuations. An interesting avenue for future work would be to extend our approach to stochastic stationary ecological dynamics. Since the “invasion implies fixation” principle has also been proven for stochastic fluctuations in the environment (Cai & Geritz, 2020), we think this extension is feasible and would help link our method with classical theory on the influence of environmental stochasticity on life‐history evolution (Frank & Slatkin, 1990; Lande et al., 2017; Sasaki & Ellner, 1995). However, this would require a careful definition of the concept of time‐dependent reproductive value in stochastic environments.

Our general recipe to study adaptation in periodically fluctuating environments relies on the classical Adaptive Dynamics assumption that the mutation rate is small. This analysis may predict the evolution of generalist strategies that balance the exploitation of the habitats that fluctuate periodically. Yet, as the period of the fluctuation between different environments increases, one may also expect that new mutations will introduce enough genetic variation to allow the population to adapt to this time‐varying environment through the recurrent selection of genotypes specialised to each habitat. This would be an interesting avenue for future work.

Although we assumed in our scenarios that all life‐history traits are constant, so that temporal fluctuations only come from periodic variations in birth rates which cause fluctuations in densities, it is straightforward to consider time‐dependent traits (such as a seasonal transmissibility ). This would give rise to additional terms capturing the temporal covariance between the trait and some life‐cycle‐specific ecological variable (Cornet et al., 2014; Kamo & Sasaki, 2005; Koelle et al., 2005; Pigeault et al., 2018).

More generally, fluctuating environments often select for plasticity in life‐history traits, which allows an organism to switch between phenotypes specialised to each habitat. It would be interesting to take this into account in our scenarios, especially since there are numerous examples of plastic life‐history strategies in pathogens. For instance, malaria parasites have evolved plastic transmission strategies to cope with the fluctuations of their within‐host environment as well as the fluctuations of the availability of their mosquito vectors (Birget et al., 2017; Cornet et al., 2014; Mideo & Reece, 2012). In addition, many viruses have the ability to change their host exploitation strategy when they perceive a change in the within‐cell environment or other cues from the environment (Gandon, 2016). Our approach can be a valuable tool to understand the evolution of fascinating strategies that enable viruses to coordinate the exploitation of their host population in fluctuating environments (see Bruce et al. (2021) for a recent application of our method).

DATA AND CODE ACCESSIBILITY

The Mathematica notebooks and R files used to run the numerical simulations and generate the figures are accessible from the following data repository, https://doi.org/10.5281/zenodo.6469843.

FUNDING INFORMATION

This work was funded by grants ANR‐16‐CE35‐0012 “STEEP” to S.L. and ANR‐17‐CE35‐0012 “EVOMALWILD” to S.G. from the Agence Nationale de la Recherche.

AUTHOR CONTRIBUTIONS

Both authors contributed equally to this work.

Associate Editor: M. Boots

Handling Editor: T. Chapman

Figure S.1: Scenario 1: The Curse of Pharaoh.

Figure S.2: Scenario 2: host preference.

Figure S.3: Scenario 3: imperfect vaccine and virulence.

Figure S.4: Scenario 3: imperfect vaccine and virulence.

Click here for additional data file.