Introducing evolution

This chapter could equally well belong in the population ecology section, or indeed a section of its own, but as I will focus on applications of evolution to disease systems we will examine it here.

In the models we have looked at so far, all individuals within a population (or compartment of a population) are identical – they have the same birth rates, death rates, risk of infection and so on. However, we often see mutations arise within populations, leading to offspring – and therefore a subpopulation – with slightly different rates. These two strains of the organism – the resident and the invader – will then compete with one another, with potentially one strain dominating (such a scenario could be applied to any sort of invading species, but our focus will be on mutants generated through reproduction). This gives rise to the potential for evolution.

The basic idea of natural selection is that strains of an organism that are the most successful – that is they make the best use of their resources, live the longest and, ultimately, produce the most offspring – are likely to be better represented in future generations. Evolution is the repeated process by which new mutants are generated and these most successful strains come to dominate. We often think of evolution as an immensely slow process – for example the evolution of humans from primates – but in large populations with short generation times, microbial populations for example, it can be witnessed in experiments lasting just a few days.

Like most biological processes, evolution is fundamentally complicated. An important element of mathematical modelling is deciding what assumptions we need to include, and what we will neglect. Broadly, evolutionary models take one of two approaches. One set of models focus on the genetics behind evolution. Mutation is ultimately a genetic process, and it is the genes that determine the behaviour of organisms. However, these models can have little to say about what it is that drives evolution. The other set of models focus on evolutionary ecology; that is, how the ecology (population dynamics) drives evolution (adaptive dynamics). This is because it is the ecological environment that creates selection that favours different strains. In these models we often assume very simple genetics, but are able to gain insight in to when we might expect organisms to evolve in a certain way. It is this second approach that we shall be introducing here, and in particular an approach developed through the 1990s called adaptive dynamics or more generally evolutionary invasion analysis.

Resident dynamics

To give us some focus, we will think about evolution in the context of infectious disease (this is my own area of research, so I had to really). In particular we will consider how a parasite might evolve to best exploit its host. We start by writing down the population dynamics of the resident strain in the absence of the mutant. We will take a variation of our wildlife disease model:

[latex]\begin{align} \frac{dS}{dt} &= (b-qN)S-dS-\beta SI\\ \frac{dI}{dt} &= \beta SI-(d+\alpha)I, \end{align}[/latex]

where [latex]N=S+I[/latex] is the total population size. Notice that here reproduction of hosts is reduced due to competition ([latex]q[/latex]), creating an emergent carrying capacity for the population in the absence of disease. We will not go through the full analysis of the resident dynamics here as it is very similar to the wildlife disease model, but you are encouraged to carry this out yourself.

The key points are:

• There is a trivial equilibrium at [latex]S = I = 0[/latex]. This equilibrium is never stable and biologically fairly irrelevant.
• There is a disease-free equilbirium at [latex]S_{df} = (b -d)/q, I_{df} = 0[/latex]. This is stable for [latex]R_0=\frac{\beta S_{df}}{d+\alpha}\lt 1[/latex].
• There is an endemic equilibrium at [latex]S_{e} = (d + \alpha)/\beta, I_{e} = (b -qS_e -d)/(\beta + q)[/latex]. This is stable for [latex]R_0=\frac{\beta S_{df}}{d+\alpha}\gt 1[/latex], which is true whenever [latex]I_e\gt0[/latex].

Fitness

We will assume that an organism, here a parasite, is able to evolve one or more traits. Pragmatically, this means that one or more parameters of the model will evolve. The key quantity for analysing evolution is fitness. This can be formally defined in slightly different ways, but it is ultimately a measure of how successful a given strain is in its environment. For our model, the fitness of a particular strain will be defined as its exponential growth rate.

Let us assume that initially all individuals in the population are the resident strain. We will then try and determine the fitness of a mutant strain, which has a small difference in one or more traits (we will usually assume two traits vary as we expect there to be a trade-off: if a strain is ‘stronger’ in one area we expect it to pay for this by being ‘weaker’ in another. However, these two traits will be linked by some defined trade-off function, meaning we can ultimately reduce it down to thinking about what happens to just one trait). We will make three key assumptions:

1. Mutants arise very rarely, meaning we can assume that the resident has reached an equilibrium of its population dynamics. This is called a separation of timescales assumption; the population dynamics are ‘fast’ and the adaptive dynamics are ‘slow’.
2. The mutant is initially rare, meaning it has a very small density (initially) relative to the resident.
3. The mutant is very similar to the resident, with just a small difference in the value of the parameter(s) that evolve.

The mutant’s growth rate depends both on its own strategy and on the ‘environment’, by which we mean the ecological conditions created by the current resident. Let us say that the current resident has some strategy [latex]x[/latex]. We can write the fitness of this resident – that is its growth rate – as [latex]r(x,E_x)[/latex], where [latex]E_x[/latex] represents the environment. In this case we must have [latex]r=0[/latex], because we have assumed the resident is at equilibrium, so it should be neither growing nor shrinking.

Now consider a mutant type with strategy [latex]x_m[/latex] invading a resident environment. The fitness of this mutant is [latex]r(x_m,E_x)[/latex]. In all likelihood [latex]r\neq0[/latex]. If [latex]r\gt0[/latex] then the mutant population will grow, meaning the mutant can invade the current resident. In that case it may coexist with the old resident, or more commonly will completely replace it. If [latex]r\lt 0[/latex] it will die out.

Remember we assumed mutants are very similar to the resident (so the mutation size is small). This means we can take a Taylor expansion of the mutant’s fitness up to the linear term as,
[latex]\begin{equation} r(x_m,E_x)=r(x,E_x)+\frac{\partial r}{\partial x_m}\Bigg|_{x_m=x}(x_m-x)+... \end{equation}[/latex]

Since [latex]r(x,E_x)=0[/latex], the direction of selection – whether it is positive or negative – is entirely dictated by the fitness gradient, [latex][\partial r/\partial x_m]_{x_m=x}[/latex]. If the mutant has a trait value [latex]x_m>x[/latex] then for that mutant to invade there must be a positive fitness gradient.

We can picture there being a sequence of these invasion events by mutants. We start with resident type A, and assume a mutant of type B arises. Where that mutant is successful it will usually go on to replace the type A population and now the resident population is type B. Then a new mutant of type C arises and attempts to invade, and so on. Given small mutations, we can now think of the dynamics of the evolving trait as being a dynamical system in its own right; the adaptive dynamics. If we define a (rather vague) evolutionary timescale, [latex]T\gt\gt t[/latex] (with [latex]t[/latex] the ecological timescale), then the change in the trait value [latex]x[/latex] will be given by,
[latex]\begin{equation} \frac{dx}{dT}=\mu \frac{\partial r}{\partial x_m}\Bigg|_{x_m=x}. \end{equation}[/latex]

The parameter [latex]\mu\gt0[/latex] contains information about how quickly new mutants appear (it is beyond the scope of this course to consider this in detail, so we shall just assume it is a constant). Since the fitness gradient will depend on the resident equilibrium, you might begin to see how the separation of timescales assumption works: we first have to solve the population dynamics to its equilibrium, then solve the adaptive dynamics. Note that [latex]\mu[/latex] only scales the ‘speed’ of evolution. Therefore if only one organism is evolving (so relative speeds do not matter) we will generally only need to focus on the fitness gradient to determine the long-term evolutionary behaviour.

Given this definition, we might then expect to look for ‘equilibria’ of these adaptive dynamics, like we have in all of our population models so far. That is, we would expect evolution of the trait [latex]x[/latex] to continue until a point is reached where the fitness gradient is zero. This is indeed what we do, but things are a little more complicated and there is some new terminology associated with it. Firstly, an ‘equilibrium’ of adaptive dynamics is known as a singular strategy,
[latex]\begin{equation} \frac{dx}{dT}\propto\frac{\partial r}{\partial x_m}\Big|_{x_m=x}= 0. \end{equation}[/latex]

What can happen at such a singular strategy? Perhaps confusingly, there are now two different types of stability we need to consider.

Evolutionary stability

Perhaps even more confusingly, our first term, evolutionary stability, is not really ‘stability’ in a traditional mathematical sense at all. Instead it relates to a long-standing ecological idea of the Evolutionarily Stable Strategy (ESS). In fact this term asks, if the strategy is adopted by the current resident, can it be invaded by any (nearby) mutants? Mathematically, this is given by,
[latex]\begin{equation} \frac{\partial^2 r}{\partial x_m^2}\Big|_{x_m=x}\lt0 \end{equation}[/latex]
You might spot that the ‘peak’ or ‘trough’ will always be at [latex]r=0[/latex], so all nearby mutants either have positive or negative fitness. In other words, it asks whether the strategy is a local fitness maximum or not.

Convergence stability

However, just because a strategy is uninvadable that does not mean that it will be reached through the evolutionary process of mutations and invasions. Instead we need to assess whether the strategy is locally attracting. This term is called convergence stability and is a more direct counterpart of classic stability.
[latex]\begin{equation} \frac{d}{dx}\left(\frac{\partial r}{\partial x_m}\Big|_{x_m=x}\right)=\frac{\partial^2 r}{\partial x_m^2}\Big|_{x_m=x}+\frac{\partial^2 r}{\partial x_m\partial x}\Big|_{x_m=x}\lt 0 \end{equation}[/latex]

Note that this is the sum of the evolutionary stability condition and an additional second derivative.

Evolutionary stability and convergence stability are independent properties. This therefore provides four potential outcomes at an evolutionary singular point:

• A Continuously Stable Strategy (CSS) is both locally attracting and uninvadable. It is therefore a long-term end-point of evolution.
• A repeller is neither attracting nor uninvadable. Therefore evolution will always take the system away from this strategy.
• A Garden of Eden is uninvadable and so would be a local fitness maximum were it ever reached. However, since it is not attracting, evolution will always take the system away from this strategy. In that sense it usually behaves roughly the same as a repeller.
• An evolutionary branching point is locally attracting but invadable, that is a fitness minimum. In this case evolution will drive the population to this strategy, but once there all nearby mutants can invade. What happens then is that the population divides in to two coexisting strategies either side of the singular strategy (subject to a couple of extra assumptions that we’ll ignore here), and evolution continues. This is an important outcome since it creates diversity directly from the evolutionary process.

The evolution of parasite virulence

Let us take these theoretical ideas and apply them to our infectious disease model to see what it all really means. We will assume here that it is the parasite that is able to evolve, and in particular that it can evolve to alter the amount of transmission, [latex]\beta[/latex]. We might well expect a disease to ‘want’ to infect and spread as quickly as possible (though we should do our best not to be quite so anthropomorphic). A classic assumption in the theoretical literature is that there is a transmission-virulence (mortality) trade-off for parasites. The argument is that for a parasite to increase its transmission rate it must grow more quickly, and this faster growth is likely to cause greater damage and hence greater virulence. As such we might define a trade-off function [latex]\alpha=\alpha(\beta)[/latex] which is an increasing function.

We shall assume that we have parameter values such that [latex]R_0\gt1[/latex] and the endemic equilibrium is stable. How do we define the fitness for the parasite? In this model we do not track explicitly the number of parasites, so we don’t have a growth rate for them as such. However, we do track the density of host individuals infected with parasites, and we can argue it is reasonable to say that this number of infected hosts is a good measure of the growth and success of the parasite, and therefore we will take the growth rate of infected hosts to be the fitness. We then need to write out the dynamics of the mutant parasite, [latex]I_m[/latex], interacting with a population of resident hosts and parasites. These are given by,

[latex]\begin{equation} \frac{dI_m}{dt}=\beta_m S_eI_m - (\alpha(\beta_m)+d)I_m. \end{equation}[/latex]

If we had any terms with mutant densities higher than first order we could set these to be 0, since the mutant is assumed to be rare. As we can see, that is not the case here (but it often will be). As such to find the fitness for the mutant parasite, i.e. its growth rate, we can write [latex]dI_m/dt=rI_m[/latex], where [latex]r[/latex] is the growth rate given by,

[latex]\begin{equation} r=\beta_m S_e - (\alpha(\beta_m)+d) \end{equation}[/latex]

The next step is to find the fitness gradient. We take the derivative of [latex]r[/latex] with respect to [latex]\beta_m[/latex] and then substitute in [latex]\beta_m=\beta[/latex] (this substitution is about having small mutation steps). This fitness gradient gives,

[latex]\begin{equation} \frac{\partial r}{\partial\beta_m}\Big|_{\beta_m=\beta}= S_e - \alpha'(\beta) \end{equation}[/latex]

If this fitness gradient is positive the parasite will increase the transmission rate, and hence also increase its virulence; if the fitness gradient is negative the parasite will reduce the transmission rate, and hence also decrease the virulence. Which occurs depends on the steepness of the transmission-virulence trade-off [latex]\alpha(\beta)[/latex]. Where the two terms are equal we will have a singular strategy.

Let us assume we are at a singular point, so that must mean [latex]S_e=\alpha'(\beta)[/latex]. What are the stability conditions at this singular point? Firstly we have evolutionary stability,

[latex]\begin{equation} \frac{\partial^2 r}{\partial\beta_m^2}\Big|_{\beta_m=\beta} = - \alpha''(\beta) \end{equation}[/latex]

So we have that the singular strategy will be evolutionarily stable whenever the curvature (second derivative) of the trade-off at the singular strategy is positive – this means that increased transmission is ‘increasingly costly’ in terms of virulence. For convergence stability we have the sum of the last term and,

[latex]\begin{equation} \frac{\partial^2 r}{\partial\beta_m\beta}\Big|_{\beta_m=\beta} = \frac{d{S_e}}{d\beta}. \end{equation}[/latex]

What does this term evaluate to? Using the expression we found for [latex]S_e[/latex], we have,

[latex]\begin{align} \frac{d S_e}{d\beta} &= \frac{\alpha'(\beta)}{\beta}-\frac{\alpha(\beta)+d}{\beta^2}\\ &=\frac{1}{\beta}\left(\alpha'(\beta)-\frac{\alpha(\beta)+d}{\beta}\right)\\ &=\frac{1}{\beta}\left(\alpha'(\beta)-S_e\right) \end{align}[/latex]

Now recall that we are evaluating this second derivative at the singular strategy, and we found above that at the singular strategy we have that [latex]S_e- \alpha'(\beta)=0[/latex]. Therefore the term above also evaluates to zero. The final convergence stability term, then, is identical to the evolutionary stability term,

[latex]\begin{equation} \frac{\partial^2 r}{\partial\beta_m^2}\Big|_{\beta_m=\beta}+\frac{\partial^2 r}{\partial\beta_m\partial\beta}\Big|_{\beta_m=\beta}=-\alpha''(\beta) \end{equation}[/latex]

So what does this tell us? That whenever the curvature of the trade-off is positive, then the singular strategy is both evolutionarily stable and convergence stable, and therefore a CSS. In contrast, whenever the curvature is negative, then the singular strategy is a repeller. Therefore we can never get evolutionary branching in this model (but will happen in many other models).

How does optimal parasite virulence vary?

Now let’s go to a bit more detail. Let us assume we have chosen parameters and a trade-off such that we predict a singular strategy that is a Continuously Stable Strategy (CSS). Now let us ask what would happen if we changed the system by increasing the host’s birth rate, [latex]b[/latex]. What would we expect to happen? Should the parasite increase or decrease its transmission (and therefore virulence)?

The answer, in fact, is neither. Changing the host’s birth rate will have no effect on parasite evolution. This may initially seem puzzling. After all, we have changed the environment (we would presume by increasing the size of the host population), and we know that the host birth rate appears in the population dynamics. Let’s look at the parasite’s fitness gradient, which is the driver of evolution,

[latex]\begin{equation} \frac{\partial r}{\partial\beta_m}\Big|_{\beta_m=\beta}= S_e - \alpha'(\beta)=\frac{\alpha(\beta)+d}{\beta}-\alpha'(\beta). \end{equation}[/latex]

As we can see, the host birth rate does not appear anywhere in the fitness gradient, and this is why it does not impact parasite evolution. Parasite evolution is simply determined by how infective the mutant is compared to the resident and how much more quickly it dies. We might have expected increasing the birth rate to lead to more infections since there are more hosts. However, it is only the equilibrium density of susceptible hosts that matters, and this value is actually not impacted by host births at all.

Now let us ask how parasite evolution might depend on the (natural) death rate, [latex]d[/latex]. As we can see, this parameter does appear in the parasite’s fitness since it partly controls the susceptible density. Let’s make our lives easier by choosing a functional form for the trade-off as,

[latex]\begin{equation} \alpha(\beta)=\alpha_{\text{min}}+\beta^2 \end{equation}[/latex]

(This is actually not a very good choice of trade-off form in general – think what happens to its gradient when [latex]\beta=0[/latex] – but it is good enough for now). Since [latex]\alpha'(\beta)\gt0[/latex] we know from above that this will produce a CSS (both evolutionarily stable and convergence stable). The fitness now becomes,

[latex]\begin{equation} \frac{\partial r}{\partial\beta_m}\Big|_{\beta_m=\beta}= S_e - \alpha'(\beta)=\frac{\alpha_{\text{min}}+\beta^2+d}{\beta}-2\beta. \end{equation}[/latex]

We can solve this to find the singular strategy,
[latex]\begin{align} \frac{\partial r}{\partial\beta_m}\Big|_{\beta_m=\beta} &= 0\\ \implies \frac{\alpha_{\text{min}}+\beta^2+d}{\beta}-2\beta &=0\\ \implies \beta&=\sqrt{\alpha_{\text{min}}+d} \end{align}[/latex]

Therefore the parasite’s evolutionarily stable transmission rate is an increasing function of the host’s death rate. In other words, we expect the parasite to evolve higher transmission – and therefore higher virulence – against hosts that die more quickly.

Why might this be? In general, what might we expect the parasite’s optimum strategy to be? Recall we defined the term [latex]R_0[/latex] to be a measure of how infective the parasite is. We would expect the parasite to evolve to be as infectious as possible, and therefore to maximise [latex]R_0[/latex]. A knock-on effect of this is that we should therefore expect the parasite to minimise the susceptible population, as this would mean it is exploiting the host to the best of its ability. Remember that [latex]S_e=S_{df}/R_0[/latex], clearly highlighting this link.

Looking at our equilibrium expression, we see that increasing the death rate, [latex]d[/latex], leads to a bigger susceptible density, [latex]S_e[/latex] (and so a lower [latex]R_0[/latex]). Since the parasite’s optimum strategy should be to minimise [latex]S_e[/latex] (maximise [latex]R_0[/latex]), we would therefore expect evolution to cause the parasite to exploit the host more. That means evolving higher transmission, and by the trade-off, higher virulence. Note that this means than when hosts die more quickly due to natural causes, the parasite is evolving to increase the disease-induced mortality as well.

Chapter references

• The theoretical underpinnings are developed from the paper by Geritz et al. (1998).